[cig-commits] r17148 - in short/3D/PyLith/branches/v1.5-stable: libsrc/meshio pylith/feassemble unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Mon Aug 30 12:00:31 PDT 2010
Author: brad
Date: 2010-08-30 12:00:30 -0700 (Mon, 30 Aug 2010)
New Revision: 17148
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
Log:
Fixed discrepancy in vertex order between CUBIT, PyLith, FIAT, and VTK.
Modified: short/3D/PyLith/branches/v1.5-stable/libsrc/meshio/MeshIOCubit.cc
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/libsrc/meshio/MeshIOCubit.cc 2010-08-30 15:49:48 UTC (rev 17147)
+++ short/3D/PyLith/branches/v1.5-stable/libsrc/meshio/MeshIOCubit.cc 2010-08-30 19:00:30 UTC (rev 17148)
@@ -396,6 +396,8 @@
; // do nothing
} else if (3 == meshDim && 8 == numCorners) { // HEX8
; // do nothing
+ } else if (2 == meshDim && 6 == numCorners) { // TRI6
+ ; // do nothing
} 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-08-30 15:49:48 UTC (rev 17147)
+++ short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATSimplex.py 2010-08-30 19:00:30 UTC (rev 17148)
@@ -93,42 +93,64 @@
"""
self._setupGeometry(spaceDim)
- if "point" != self.shape.lower():
+ if "point" == self.shape.lower():
+ # Need 0-D quadrature for boundary conditions of 1-D meshes
+ self.cellDim = 0
+ self.numCorners = 1
+ self.numQuadPts = 1
+ self.basis = numpy.array([[1.0]])
+ self.basisDeriv = numpy.array([[[1.0]]])
+ self.quadPts = numpy.array([[0.0]])
+ self.quadWts = numpy.array([1.0])
+ self.vertices = numpy.array([[0.0]])
+ else:
quadrature = self._setupQuadrature()
basisFns = self._setupBasisFns()
# Get coordinates of vertices (dual basis)
- self.vertices = numpy.array(self._setupVertices(), dtype=numpy.float64)
+ vertices = numpy.array(self._setupVertices(), dtype=numpy.float64)
# Evaluate basis functions at quadrature points
quadpts = quadrature.get_points()
basis = numpy.array(basisFns.tabulate(quadpts)).transpose()
- self.basis = numpy.reshape(basis.flatten(), basis.shape)
# Evaluate derivatives of basis functions at quadrature points
import FIAT.shapes
dim = FIAT.shapes.dimension(basisFns.base.shape)
basisDeriv = numpy.array([basisFns.deriv_all(d).tabulate(quadpts) \
for d in range(dim)]).transpose()
- self.basisDeriv = numpy.reshape(basisDeriv.flatten(), basisDeriv.shape)
- self.quadPts = numpy.array(quadrature.get_points())
- self.quadWts = numpy.array(quadrature.get_weights())
-
self.cellDim = dim
self.numCorners = len(basisFns)
self.numQuadPts = len(quadrature.get_weights())
- else:
- # Need 0-D quadrature for boundary conditions of 1-D meshes
- self.cellDim = 0
- self.numCorners = 1
- self.numQuadPts = 1
- self.basis = numpy.array([[1.0]])
- self.basisDeriv = numpy.array([[[1.0]]])
- self.quadPts = numpy.array([[0.0]])
- self.quadWts = numpy.array([1.0])
- self.vertices = numpy.array([[0.0]])
+ 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)
+
+ 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: ")
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-08-30 15:49:48 UTC (rev 17147)
+++ short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATSimplex.py 2010-08-30 19:00:30 UTC (rev 17148)
@@ -198,9 +198,9 @@
vertices = numpy.array([[-1.0, -1.0],
[+1.0, -1.0],
[-1.0, +1.0],
+ [ 0.0, -1.0],
[ 0.0, 0.0],
- [-1.0, 0.0],
- [ 0.0, -1.0]])
+ [-1.0, 0.0]])
quadPts = numpy.array([ [-0.64288254, -0.68989795],
[-0.84993778, 0.28989795],
[ 0.33278049, -0.68989795],
@@ -267,33 +267,33 @@
def N3(self, p):
- return (1.0+p[0])*(1+p[1])
+ return (-p[0]-p[1])*(1+p[0])
def N3p(self, p):
- return (1+p[1])
+ return -1.0-2*p[0]-p[1]
def N3q(self, p):
- return (1.0+p[0])
+ return -(1+p[0])
def N4(self, p):
- return (-p[0]-p[1])*(1+p[1])
+ return (1.0+p[0])*(1+p[1])
def N4p(self, p):
- return -(1+p[1])
+ return (1+p[1])
def N4q(self, p):
- return -1.0-p[0]-2*p[1]
+ return (1.0+p[0])
def N5(self, p):
- return (-p[0]-p[1])*(1+p[0])
+ return (-p[0]-p[1])*(1+p[1])
def N5p(self, p):
- return -1.0-2*p[0]-p[1]
+ return -(1+p[1])
def N5q(self, p):
- return -(1+p[0])
+ return -1.0-p[0]-2*p[1]
# ----------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list