[cig-commits] r17166 - in short/3D/PyLith/branches/v1.5-stable: libsrc/meshio pylith/feassemble unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu Sep 2 17:23:39 PDT 2010
Author: brad
Date: 2010-09-02 17:23:39 -0700 (Thu, 02 Sep 2010)
New Revision: 17166
Modified:
short/3D/PyLith/branches/v1.5-stable/libsrc/meshio/MeshIOCubit.cc
short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATSimplex.py
short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATSimplex.py
short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestMeshQuadrature.py
short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestSubMeshQuadrature.py
Log:
Updated to Sieve BFS ordering.
Modified: short/3D/PyLith/branches/v1.5-stable/libsrc/meshio/MeshIOCubit.cc
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/libsrc/meshio/MeshIOCubit.cc 2010-09-02 23:05:04 UTC (rev 17165)
+++ short/3D/PyLith/branches/v1.5-stable/libsrc/meshio/MeshIOCubit.cc 2010-09-03 00:23:39 UTC (rev 17166)
@@ -394,10 +394,35 @@
if (2 == meshDim && 4 == numCorners) { // QUAD4
; // do nothing
+
} else if (3 == meshDim && 8 == numCorners) { // HEX8
; // do nothing
+
} else if (2 == meshDim && 6 == numCorners) { // TRI6
- ; // do nothing
+ // CUBIT
+ // corners,
+ // bottom edges, middle edges, top edges
+
+ // PyLith (Sieve)
+ // bottom edge, right edge, left edge, corners
+
+ // Permutation: 3, 4, 5, 0, 1, 2
+ int tmp = 0;
+ for (int iCell=0; iCell < numCells; ++iCell) {
+ const int ii = iCell*numCorners;
+ tmp = (*cells)[ii+0];
+ (*cells)[ii+0] = (*cells)[ii+3];
+ (*cells)[ii+3] = tmp;
+
+ tmp = (*cells)[ii+1];
+ (*cells)[ii+1] = (*cells)[ii+4];
+ (*cells)[ii+4] = tmp;
+
+ tmp = (*cells)[ii+2];
+ (*cells)[ii+2] = (*cells)[ii+5];
+ (*cells)[ii+5] = tmp;
+ } // for
+
} else if (3 == meshDim && 27 == numCorners) { // HEX27
// CUBIT
// corners,
Modified: short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATSimplex.py
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATSimplex.py 2010-09-02 23:05:04 UTC (rev 17165)
+++ short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATSimplex.py 2010-09-03 00:23:39 UTC (rev 17166)
@@ -124,33 +124,18 @@
self.numCorners = len(basisFns)
self.numQuadPts = len(quadrature.get_weights())
- if 1 == self.degree or "line" == self.shape.lower():
- self.vertices = vertices
- self.basis = numpy.reshape(basis.flatten(), basis.shape)
- self.basisDeriv = numpy.reshape(basisDeriv.flatten(), basisDeriv.shape)
+ # Permute from FIAT order to Sieve order
+ p = self._permutationFIATToSieve()
+ self.vertices = vertices[p,:]
+ self.basis = numpy.reshape(basis[:,p].flatten(), basis.shape)
+ self.basisDeriv = numpy.reshape(basisDeriv[:,p,:].flatten(),
+ basisDeriv.shape)
- self.quadPts = numpy.array(quadrature.get_points())
- self.quadWts = numpy.array(quadrature.get_weights())
+ # No permutation in order of quadrature points
+ self.quadPts = numpy.array(quadrature.get_points())
+ self.quadWts = numpy.array(quadrature.get_weights())
- elif 2 == self.degree:
- # Set order of vertices and basis functions.
- if "triangle" == self.shape.lower():
- v = [0, 1, 2, 5, 3, 4]
- elif "tetrahedron" == self.shape.lower():
- v = [0, 1, 2, 3, 6, 4, 5, 7, 8, 9]
- else:
- raise ValueError("Unknown cell type '"+self.shape.lower()+"'.")
- self.vertices = vertices[v,:]
- self.basis = numpy.reshape(basis[:,v].flatten(), basis.shape)
- self.basisDeriv = numpy.reshape(basisDeriv[:,v,:].flatten(),
- basisDeriv.shape)
-
- self.quadPts = numpy.array(quadrature.get_points())
- self.quadWts = numpy.array(quadrature.get_weights())
- else:
- raise NotImplemented("Degree "+str(self.degree)+" not implemented.")
-
self._info.line("Cell geometry: ")
self._info.line(self.geometry)
self._info.line("Vertices: ")
@@ -181,7 +166,50 @@
self.order = self.degree
return
+
+ def _permutationFIATToSieve(self):
+ """
+ Permute from FIAT basis order to Sieve basis order.
+ FIAT: corners, edges, faces
+
+ Sieve: breadth search first (faces, edges, corners)
+ """
+ import FIAT.shapes
+
+ basis = self.cell.function_space()
+ dim = FIAT.shapes.dimension(self._getShape())
+ ids = self.cell.Udual.entity_ids
+ permutation = []
+ if dim == 1:
+ for edge in ids[1]:
+ permutation.extend(ids[1][edge])
+ for vertex in ids[0]:
+ permutation.extend(ids[0][vertex])
+ elif dim == 2:
+ for face in ids[2]:
+ permutation.extend(ids[2][face])
+ for edge in ids[1]:
+ permutation.extend(ids[1][(edge+2)%3])
+ for vertex in ids[0]:
+ permutation.extend(ids[0][vertex])
+ elif dim == 3:
+ for volume in ids[3]:
+ permutation.extend(ids[3][volume])
+ for face in [3, 2, 0, 1]:
+ permutation.extend(ids[2][face])
+ for edge in [2, 0, 1, 3]:
+ permutation.extend(ids[1][edge])
+ for edge in [4, 5]:
+ if len(ids[1][edge]) > 0:
+ permutation.extend(ids[1][edge][::-1])
+ for vertex in ids[0]:
+ permutation.extend(ids[0][vertex])
+ else:
+ raise ValueError("Unknown dimension '%d'." % dim)
+ return permutation
+
+
def _setupGeometry(self, spaceDim):
"""
Setup reference cell geometry object.
@@ -215,6 +243,7 @@
if None == self.geometry:
raise ValueError("Could not set shape '%s' of cell for '%s' in spatial " \
"dimension '%s'." % (name, self.name, spaceDim))
+
return
@@ -233,15 +262,15 @@
Setup basis functions for reference cell.
"""
from FIAT.Lagrange import Lagrange
- return Lagrange(self._getShape(), self.degree).function_space()
+ self.cell = Lagrange(self._getShape(), self.degree)
+ return self.cell.function_space()
def _setupVertices(self):
"""
Setup evaluation functional points for reference cell.
"""
- from FIAT.Lagrange import Lagrange
- return Lagrange(self._getShape(), self.degree).Udual.pts
+ return self.cell.Udual.pts
def _getShape(self):
Modified: short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATSimplex.py
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATSimplex.py 2010-09-02 23:05:04 UTC (rev 17165)
+++ short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATSimplex.py 2010-09-03 00:23:39 UTC (rev 17166)
@@ -79,7 +79,7 @@
"""
Setup line3 cell.
"""
- vertices = numpy.array([[-1.0], [1.0], [0.0]])
+ vertices = numpy.array([[0.0], [-1.0], [1.0]])
quadPts = numpy.array([ [-1.0/3**0.5],
[+1.0/3**0.5] ])
quadWts = numpy.array( [1.0, 1.0])
@@ -91,7 +91,9 @@
for q in quadPts:
basis[iQuad] = numpy.array([self.N0(q), self.N1(q), self.N2(q)],
dtype=numpy.float64).reshape( (3,) )
- deriv = numpy.array([[self.N0p(q)], [self.N1p(q)], [self.N2p(q)]])
+ deriv = numpy.array([[self.N0p(q)],
+ [self.N1p(q)],
+ [self.N2p(q)]])
basisDeriv[iQuad] = deriv.reshape((3, 1))
iQuad += 1
@@ -107,23 +109,26 @@
def N0(self, p):
- return -0.5*p*(1.0-p)
+ return (1.0-p*p)
def N0p(self, p):
- return 1.0*p - 0.5
+ return -2.0*p
+
def N1(self, p):
- return 0.5*p*(1.0+p)
+ return -0.5*p*(1.0-p)
def N1p(self, p):
- return 1.0*p + 0.5
+ return 1.0*p - 0.5
+
def N2(self, p):
- return (1.0-p*p)
+ return 0.5*p*(1.0+p)
def N2p(self, p):
- return -2.0*p
+ return 1.0*p + 0.5
+
# ----------------------------------------------------------------------
class Tri3(object):
@@ -170,6 +175,7 @@
def N0q(self, p):
return -0.5
+
def N1(self, p):
return 0.5*(1.0+p[0])
@@ -179,6 +185,7 @@
def N1q(self, p):
return 0.0
+
def N2(self, p):
return 0.5*(1.0+p[1])
@@ -188,6 +195,7 @@
def N2q(self, p):
return 0.5
+
# ----------------------------------------------------------------------
class Tri6(object):
@@ -195,12 +203,12 @@
"""
Setup tri33 cell.
"""
- vertices = numpy.array([[-1.0, -1.0],
+ vertices = numpy.array([[ 0.0, -1.0],
+ [ 0.0, 0.0],
+ [-1.0, 0.0],
+ [-1.0, -1.0],
[+1.0, -1.0],
- [-1.0, +1.0],
- [ 0.0, -1.0],
- [ 0.0, 0.0],
- [-1.0, 0.0]])
+ [-1.0, +1.0]])
quadPts = numpy.array([ [-0.64288254, -0.68989795],
[-0.84993778, 0.28989795],
[ 0.33278049, -0.68989795],
@@ -237,63 +245,63 @@
def N0(self, p):
- return 0.5*(-p[0]-p[1])*(-1.0-p[0]-p[1])
+ return (-p[0]-p[1])*(1+p[0])
def N0p(self, p):
- return 0.5+p[0]+p[1]
+ return -1.0-2*p[0]-p[1]
def N0q(self, p):
- return 0.5+p[0]+p[1]
+ return -(1+p[0])
def N1(self, p):
- return 0.5*(1.0+p[0])*(p[0])
+ return (1.0+p[0])*(1+p[1])
def N1p(self, p):
- return 0.5+p[0]
+ return (1+p[1])
def N1q(self, p):
- return 0
+ return (1.0+p[0])
def N2(self, p):
- return 0.5*(1.0+p[1])*(p[1])
+ return (-p[0]-p[1])*(1+p[1])
def N2p(self, p):
- return 0
+ return -(1+p[1])
def N2q(self, p):
- return 0.5+p[1]
+ return -1.0-p[0]-2*p[1]
def N3(self, p):
- return (-p[0]-p[1])*(1+p[0])
+ return 0.5*(-p[0]-p[1])*(-1.0-p[0]-p[1])
def N3p(self, p):
- return -1.0-2*p[0]-p[1]
+ return 0.5+p[0]+p[1]
def N3q(self, p):
- return -(1+p[0])
+ return 0.5+p[0]+p[1]
def N4(self, p):
- return (1.0+p[0])*(1+p[1])
+ return 0.5*(1.0+p[0])*(p[0])
def N4p(self, p):
- return (1+p[1])
+ return 0.5+p[0]
def N4q(self, p):
- return (1.0+p[0])
+ return 0
def N5(self, p):
- return (-p[0]-p[1])*(1+p[1])
+ return 0.5*(1.0+p[1])*(p[1])
def N5p(self, p):
- return -(1+p[1])
+ return 0
def N5q(self, p):
- return -1.0-p[0]-2*p[1]
+ return 0.5+p[1]
# ----------------------------------------------------------------------
@@ -384,6 +392,7 @@
def N3r(self, p):
return 0.5
+
# ----------------------------------------------------------------------
class TestFIATSimplex(unittest.TestCase):
"""
Modified: short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestMeshQuadrature.py
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestMeshQuadrature.py 2010-09-02 23:05:04 UTC (rev 17165)
+++ short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestMeshQuadrature.py 2010-09-03 00:23:39 UTC (rev 17166)
@@ -29,22 +29,22 @@
# ----------------------------------------------------------------------
def N0(p):
- return -0.5*p*(1.0-p)
+ return (1.0-p**2)
def N0p(p):
- return -0.5*(1.0-p) + 0.5*p
+ return -2.0*p
def N1(p):
- return 0.5*p*(1.0+p)
+ return -0.5*p*(1.0-p)
def N1p(p):
- return +0.5*(1.0+p) + 0.5*p
+ return -0.5*(1.0-p) + 0.5*p
def N2(p):
- return (1.0-p**2)
+ return 0.5*p*(1.0+p)
def N2p(p):
- return -2.0*p
+ return +0.5*(1.0+p) + 0.5*p
# ----------------------------------------------------------------------
class TestMeshQuadrature(unittest.TestCase):
Modified: short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestSubMeshQuadrature.py
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestSubMeshQuadrature.py 2010-09-02 23:05:04 UTC (rev 17165)
+++ short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestSubMeshQuadrature.py 2010-09-03 00:23:39 UTC (rev 17166)
@@ -29,22 +29,22 @@
# ----------------------------------------------------------------------
def N0(p):
- return -0.5*p*(1.0-p)
+ return (1.0-p**2)
def N0p(p):
- return -0.5*(1.0-p) + 0.5*p
+ return -2.0*p
def N1(p):
- return 0.5*p*(1.0+p)
+ return -0.5*p*(1.0-p)
def N1p(p):
- return +0.5*(1.0+p) + 0.5*p
+ return -0.5*(1.0-p) + 0.5*p
def N2(p):
- return (1.0-p**2)
+ return 0.5*p*(1.0+p)
def N2p(p):
- return -2.0*p
+ return +0.5*(1.0+p) + 0.5*p
# ----------------------------------------------------------------------
class TestSubMeshQuadrature(unittest.TestCase):
More information about the CIG-COMMITS
mailing list