[cig-commits] r17178 - in short/3D/PyLith/trunk: . libsrc/faults libsrc/feassemble libsrc/materials libsrc/meshio libsrc/topology playpen pylith/feassemble pylith/perf pylith/utils tests/3d tests_auto/2d tests_auto/3dnew tests_auto/3dnew/hex8 unittests/libtests/bc unittests/libtests/bc/data unittests/libtests/faults unittests/libtests/faults/data unittests/libtests/meshio/data unittests/libtests/topology unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Tue Sep 7 15:59:06 PDT 2010
Author: brad
Date: 2010-09-07 15:59:05 -0700 (Tue, 07 Sep 2010)
New Revision: 17178
Added:
short/3D/PyLith/trunk/playpen/quadratic/
short/3D/PyLith/trunk/tests/3d/plasticity/
short/3D/PyLith/trunk/tests_auto/2d/quad9/
short/3D/PyLith/trunk/tests_auto/2d/tri6/
short/3D/PyLith/trunk/tests_auto/3dnew/hex27/
short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axial_disp.spatialdb
short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp.cfg
short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp_gendb.py
short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp_soln.py
short/3D/PyLith/trunk/tests_auto/3dnew/tet10/
short/3D/PyLith/trunk/tests_auto/3dnew/tet4/
Modified:
short/3D/PyLith/trunk/configure.ac
short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/faults/TopologyOps.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshIOCubit.cc
short/3D/PyLith/trunk/libsrc/topology/SubMesh.cc
short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py
short/3D/PyLith/trunk/pylith/feassemble/FIATSimplex.py
short/3D/PyLith/trunk/pylith/perf/Mesh.py
short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py
short/3D/PyLith/trunk/tests/3d/Makefile.am
short/3D/PyLith/trunk/tests_auto/3dnew/hex8/mesh.exo
short/3D/PyLith/trunk/tests_auto/3dnew/hex8/mesh.jou
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependent.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/BoundaryMeshDataHex8.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.hh
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.hh
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataLine2.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataLine2.hh
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataQuad4.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataQuad4.hh
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTet4.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTet4.hh
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTri3.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTri3.hh
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8b.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8d.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8e.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8f.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8g.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8h.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8i.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4e.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinSrcsDataTet4.cc
short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_cell_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_vertex_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_cell_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_vertex_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_cell_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_vertex_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_cell_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_vertex_t10.vtk
short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldsNewMesh.hh
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATSimplex.py
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestMeshQuadrature.py
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestSubMeshQuadrature.py
Log:
Merge from trunk.
Modified: short/3D/PyLith/trunk/configure.ac
===================================================================
--- short/3D/PyLith/trunk/configure.ac 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/configure.ac 2010-09-07 22:59:05 UTC (rev 17178)
@@ -304,6 +304,17 @@
tests/3d/Makefile
tests/3d/matprops/Makefile
tests/3d/slipdir/Makefile
+ tests/3d/plasticity/Makefile
+ tests/3d/plasticity/threehex27/Makefile
+ tests/3d/plasticity/threehex27/config/Makefile
+ tests/3d/plasticity/threehex27/mesh/Makefile
+ tests/3d/plasticity/threehex27/results/Makefile
+ tests/3d/plasticity/threehex27/spatialdb/Makefile
+ tests/3d/plasticity/threehex8/Makefile
+ tests/3d/plasticity/threehex8/config/Makefile
+ tests/3d/plasticity/threehex8/mesh/Makefile
+ tests/3d/plasticity/threehex8/results/Makefile
+ tests/3d/plasticity/threehex8/spatialdb/Makefile
doc/Makefile
doc/developer/Makefile
doc/install/Makefile
Modified: short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -545,7 +545,7 @@
<< cellNeighbors.size() << std::endl;
throw ALE::Exception("Invalid number of neighbors");
} // if
- } else if (numCorners == 4) {
+ } else if (numCorners == 4 || numCorners == 9) {
if (cellNeighbors.size() > 4) {
std::cout << "Cell " << *c_iter
<< " has an invalid number of neighbors "
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -1374,7 +1374,7 @@
} // if
// Check quadrature against mesh
- const int numCorners = _quadrature->refGeometry().numCorners();
+ const int numCorners = _quadrature->numBasis();
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->getLabelStratum("material-id", id());
assert(!cells.isNull());
Modified: short/3D/PyLith/trunk/libsrc/faults/TopologyOps.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/TopologyOps.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/faults/TopologyOps.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -155,6 +155,7 @@
PointArray origVertices;
PointArray faceVertices;
+ faultSieve->setDebug(2);
if (!faultSieve->commRank()) {
numCorners = mesh->getNumCellCorners();
faceSize = selection::numFaceVertices(mesh);
@@ -226,12 +227,33 @@
if (dim == 0) {
f = *faceVertices.begin();
}
- if (faceSize != dim+1) {
+
+ std::cout << "dim: " << dim << ", faceSize: " << faceSize << ", numCorners: " << numCorners << std::endl;
+
+ if (2 == dim && 4 == faceSize){
if (debug) std::cout << " Adding hex face " << f << std::endl;
- ALE::SieveBuilder<ALE::Mesh>::buildHexFaces(faultSieve, orientation, dim, curElement, bdVertices, oFaultFaces, f, o);
+ ALE::SieveBuilder<ALE::Mesh>::buildHexFaces(
+ faultSieve, orientation, dim, curElement,
+ bdVertices, oFaultFaces, f, o);
+ } else if ((1 == dim && 3 == faceSize) ||
+ (2 == dim && 9 == faceSize)){
+ if (debug) std::cout << " Adding quadratic hex face " << f
+ << std::endl;
+ ALE::SieveBuilder<ALE::Mesh>::buildQuadraticHexFaces(
+ faultSieve, orientation, dim, curElement,
+ bdVertices, oFaultFaces, f, o);
+ } else if ((1 == dim && 3 == faceSize) ||
+ (2 == dim && 6 == faceSize)){
+ if (debug) std::cout << " Adding quadratic tri face " << f
+ << std::endl;
+ ALE::SieveBuilder<ALE::Mesh>::buildQuadraticTetFaces(
+ faultSieve, orientation, dim, curElement,
+ bdVertices, oFaultFaces, f, o);
} else {
if (debug) std::cout << " Adding simplicial face " << f << std::endl;
- ALE::SieveBuilder<ALE::Mesh>::buildFaces(faultSieve, orientation, dim, curElement, bdVertices, oFaultFaces, f, o);
+ ALE::SieveBuilder<ALE::Mesh>::buildFaces(
+ faultSieve, orientation, dim, curElement,
+ bdVertices, oFaultFaces, f, o);
}
faultSieve->addArrow(f, support[s]);
//faultSieve->view("");
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -172,7 +172,8 @@
assert(0 != jacobian);
assert(0 != det);
- assert(numCorners()*spaceDim() == vertices.size());
+ assert( (numCorners()*spaceDim() == vertices.size()) || // linear hex
+ ((numCorners()+19)*spaceDim() == vertices.size()) ); // quadratic hex
assert(cellDim() == location.size());
assert(spaceDim()*cellDim() == jacobian->size());
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -99,7 +99,8 @@
assert(0 != jacobian);
assert(0 != det);
- assert(numCorners()*spaceDim() == vertices.size());
+ assert( (numCorners()*spaceDim() == vertices.size()) || // linear edge
+ ((numCorners()+1)*spaceDim() == vertices.size()) ); // quadratic edge
assert(spaceDim()*cellDim() == jacobian->size());
const double x0 = vertices[0];
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -106,7 +106,8 @@
assert(0 != jacobian);
assert(0 != det);
- assert(numCorners()*spaceDim() == vertices.size());
+ assert( (numCorners()*spaceDim() == vertices.size()) || // linear
+ ((numCorners()+1)*spaceDim() == vertices.size()) ); // quadratic
assert(spaceDim()*cellDim() == jacobian->size());
const double x0 = vertices[0];
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -110,7 +110,8 @@
assert(0 != jacobian);
assert(0 != det);
- assert(numCorners()*spaceDim() == vertices.size());
+ assert( (numCorners()*spaceDim() == vertices.size()) || // linear edge
+ ((numCorners()+1)*spaceDim() == vertices.size()) ); // quadratic edge
assert(spaceDim()*cellDim() == jacobian->size());
const double x0 = vertices[0];
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -121,7 +121,8 @@
assert(0 != jacobian);
assert(0 != det);
- assert(numCorners()*spaceDim() == vertices.size());
+ assert( (numCorners()*spaceDim() == vertices.size()) || // linear quad
+ ((numCorners()+5)*spaceDim() == vertices.size()) ); // quadratic quad
assert(cellDim() == location.size());
assert(spaceDim()*cellDim() == jacobian->size());
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -131,7 +131,8 @@
assert(0 != jacobian);
assert(0 != det);
- assert(numCorners()*spaceDim() == vertices.size());
+ assert( (numCorners()*spaceDim() == vertices.size()) || // linear quad
+ ((numCorners()+5)*spaceDim() == vertices.size()) ); // quadratic quad
assert(cellDim() == location.size());
assert(spaceDim()*cellDim() == jacobian->size());
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -27,7 +27,7 @@
#include <cassert> // USES assert()
-//#define ISOPARAMETRIC
+#define ISOPARAMETRIC
// ----------------------------------------------------------------------
// Constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -27,7 +27,7 @@
#include <cassert> // USES assert()
-//#define ISOPARAMETRIC
+#define ISOPARAMETRIC
// ----------------------------------------------------------------------
// Constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -27,7 +27,7 @@
#include <cassert> // USES assert()
-//#define ISOPARAMETRIC
+#define ISOPARAMETRIC
// ----------------------------------------------------------------------
// Constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -27,7 +27,7 @@
#include <cassert> // USES assert()
-//#define ISOPARAMETRIC
+#define ISOPARAMETRIC
// ----------------------------------------------------------------------
// Constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -29,7 +29,7 @@
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error()
-//#define ISOPARAMETRIC
+#define ISOPARAMETRIC
// ----------------------------------------------------------------------
// Constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -27,7 +27,7 @@
#include <cassert> // USES assert()
-//#define ISOPARAMETRIC
+#define ISOPARAMETRIC
// ----------------------------------------------------------------------
// Constructor
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -408,6 +408,9 @@
// Initial stress and strain values
const double meanStrainInitial = (initialStrain[0] + initialStrain[1]) / 3.0;
+ // :BUG:? CHARLES - Doesn't sigma_zz contribute to
+ // meanStressInitial? Isn't sigma_zz nonzero and found from sigma_xx
+ // and sigma_yy?
const double meanStressInitial = (initialStress[0] + initialStress[1]) / 3.0;
const double devStrainInitial[] = {initialStrain[0] - meanStrainInitial,
initialStrain[1] - meanStrainInitial,
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIOCubit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIOCubit.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIOCubit.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -392,38 +392,98 @@
assert(0 != cells);
assert(cells->size() == numCells*numCorners);
- if (2 == meshDim && 4 == numCorners) // QUAD
-#if 0 // OLD
- // 0 1 2 3 -> 0 1 3 2
+ if (2 == meshDim && 4 == numCorners) { // QUAD4
+ ; // do nothing
+
+ } else if (3 == meshDim && 8 == numCorners) { // HEX8
+ ; // do nothing
+
+ } else if (2 == meshDim && 6 == numCorners) { // TRI6
+ // 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 i2 = iCell*numCorners+2;
- const int i3 = iCell*numCorners+3;
- const int tmp = (*cells)[i2];
- (*cells)[i2] = (*cells)[i3];
- (*cells)[i3] = tmp;
+ 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
- ; // do nothing
-#endif
- else if (3 == meshDim && 8 == numCorners) // HEX
-#if 0 // OLD
- // 0 1 2 3 4 5 6 7 -> 0 1 3 2 4 5 7 6
+
+ } else if (3 == meshDim && 27 == numCorners) { // HEX27
+ // CUBIT
+ // corners,
+ // bottom edges, middle edges, top edges
+ // interior
+ // bottom/top, left/right, front/back
+
+ // PyLith
+ // corners,
+ // bottom edges, top edges, middle edges
+ // left/right, front/back, bottom/top
+ // interior
+ int tmp = 0;
for (int iCell=0; iCell < numCells; ++iCell) {
- const int i2 = iCell*numCorners+2;
- const int i3 = iCell*numCorners+3;
- int tmp = (*cells)[i2];
- (*cells)[i2] = (*cells)[i3];
- (*cells)[i3] = tmp;
+ const int i12 = iCell*numCorners+12;
+ const int i13 = iCell*numCorners+13;
+ const int i14 = iCell*numCorners+14;
+ const int i15 = iCell*numCorners+15;
+ const int i16 = iCell*numCorners+16;
+ const int i17 = iCell*numCorners+17;
+ const int i18 = iCell*numCorners+18;
+ const int i19 = iCell*numCorners+19;
+ const int i20 = iCell*numCorners+20;
+ const int i21 = iCell*numCorners+21;
+ const int i22 = iCell*numCorners+22;
+ const int i23 = iCell*numCorners+23;
+ const int i24 = iCell*numCorners+24;
+ const int i25 = iCell*numCorners+25;
+ const int i26 = iCell*numCorners+26;
- const int i6 = iCell*numCorners+6;
- const int i7 = iCell*numCorners+7;
- tmp = (*cells)[i6];
- (*cells)[i6] = (*cells)[i7];
- (*cells)[i7] = tmp;
+ tmp = (*cells)[i12];
+ (*cells)[i12] = (*cells)[i16];
+ (*cells)[i16] = tmp;
+
+ tmp = (*cells)[i13];
+ (*cells)[i13] = (*cells)[i17];
+ (*cells)[i17] = tmp;
+
+ tmp = (*cells)[i14];
+ (*cells)[i14] = (*cells)[i18];
+ (*cells)[i18] = tmp;
+
+ tmp = (*cells)[i15];
+ (*cells)[i15] = (*cells)[i19];
+ (*cells)[i19] = tmp;
+
+ tmp = (*cells)[i20];
+ (*cells)[i20] = (*cells)[i23];
+ (*cells)[i23] = (*cells)[i26];
+ (*cells)[i26] = tmp;
+
+ tmp = (*cells)[i21];
+ (*cells)[i21] = (*cells)[i24];
+ (*cells)[i24] = tmp;
+
+ tmp = (*cells)[i22];
+ (*cells)[i22] = (*cells)[i25];
+ (*cells)[i25] = tmp;
} // for
-#else
- ; // do nothing
-#endif
+ } // if/else
+
} // _orientCells
Modified: short/3D/PyLith/trunk/libsrc/topology/SubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/SubMesh.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/libsrc/topology/SubMesh.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -93,8 +93,8 @@
_mesh->setRealSection("coordinates",
meshSieveMesh->getRealSection("coordinates"));
if (meshSieveMesh->hasRealSection("coordinates_dimensioned"))
- _mesh->setRealSection("coordinates_dimensioned",
- meshSieveMesh->getRealSection("coordinates_dimensioned"));
+ _mesh->setRealSection("coordinates_dimensioned",
+ meshSieveMesh->getRealSection("coordinates_dimensioned"));
// Create the parallel overlap
const ALE::Obj<SieveMesh::sieve_type>& sieve = _mesh->getSieve();
Copied: short/3D/PyLith/trunk/playpen/quadratic (from rev 17177, short/3D/PyLith/branches/v1.5-stable/playpen/quadratic)
Modified: short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py 2010-09-07 22:59:05 UTC (rev 17178)
@@ -103,14 +103,14 @@
quadwts = numpy.array(quadrature.get_weights())
numQuadPts = len(quadpts)
basis = numpy.array(element.function_space().tabulate(quadrature.get_points())).transpose()
- numBasisFns = len(element.function_space())
+ numBasis = len(element.function_space())
# Evaluate derivatives of basis functions at quadrature points
basisDeriv = numpy.array([element.function_space().deriv_all(d).tabulate(quadrature.get_points()) \
for d in range(1)]).transpose()
self.numQuadPts = numQuadPts**dim
- self.numCorners = numBasisFns**dim
+ self.numCorners = numBasis**dim
if dim == 1:
self.vertices = numpy.array(vertices)
@@ -120,520 +120,219 @@
self.basisDeriv = numpy.reshape(basisDeriv.flatten(), basisDeriv.shape)
else:
if dim == 2:
+ # Set order of vertices and basis functions.
+ # Corners
+ vertexOrder = [(0,0), (1,0), (1,1), (0,1)]
+ # Edges
+ # Bottom
+ p = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q)
+ # Right
+ p = numpy.ones(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q)
+ # Top
+ p = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ q = numpy.ones(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q)
+ # Left
+ p = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ vertexOrder += zip(p,q)
+ # Face
+ p = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q)
+
self.vertices = numpy.zeros((self.numCorners, dim))
n = 0
- # Bottom
- for r in range(0, numBasisFns-1):
- self.vertices[n][0] = vertices[r]
- self.vertices[n][1] = vertices[0]
- n += 1
- # Right
- for q in range(0, numBasisFns-1):
- self.vertices[n][0] = vertices[numBasisFns-1]
+ for (p,q) in vertexOrder:
+ self.vertices[n][0] = vertices[p]
self.vertices[n][1] = vertices[q]
n += 1
- # Top
- for r in range(numBasisFns-1, 0, -1):
- self.vertices[n][0] = vertices[r]
- self.vertices[n][1] = vertices[numBasisFns-1]
- n += 1
- # Left
- for q in range(numBasisFns-1, 0, -1):
- self.vertices[n][0] = vertices[0]
- self.vertices[n][1] = vertices[q]
- n += 1
- # Interior
- for q in range(1, numBasisFns-1):
- for r in range(1, numBasisFns-1):
- self.vertices[n][0] = vertices[r]
- self.vertices[n][1] = vertices[q]
- n += 1
- if not n == self.numCorners: raise RuntimeError('Invalid 2D function tabulation')
+ if not n == self.numCorners:
+ raise RuntimeError('Invalid 2-D vertex ordering: '+str(n)+ \
+ ' should be '+str(self.numCorners))
self.quadPts = numpy.zeros((numQuadPts*numQuadPts, dim))
self.quadWts = numpy.zeros((numQuadPts*numQuadPts,))
self.basis = numpy.zeros((numQuadPts*numQuadPts,
- numBasisFns*numBasisFns))
+ numBasis*numBasis))
self.basisDeriv = numpy.zeros((numQuadPts*numQuadPts,
- numBasisFns*numBasisFns, dim))
+ numBasis*numBasis, dim))
+
+ # Order of quadrature points doesn't matter
+ # Order of basis functions should match vertices for isoparametric
n = 0
- # Bottom
- for r in range(0, numQuadPts-1):
- self.quadPts[n][0] = quadpts[r]
- self.quadPts[n][1] = quadpts[0]
- self.quadWts[n] = quadwts[r]*quadwts[0]
- m = 0
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[0][0]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[0][0]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[0][0][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[0][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[0][f]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[0][f][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[0][numBasisFns-1]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[0][numBasisFns-1]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[0][numBasisFns-1][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[0][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[0][f]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[0][f][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[0][f]
- self.basisDeriv[0][r][f][g][0] = basisDeriv[r][g][0]*basis[0][f]
- self.basisDeriv[0][r][f][g][1] = basis[r][g]*basisDeriv[0][f][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 2D function tabulation')
- n += 1
- # Right
- for q in range(0, numQuadPts-1):
- self.quadPts[n][0] = quadpts[numQuadPts-1]
- self.quadPts[n][1] = quadpts[q]
- self.quadWts[n] = quadwts[numQuadPts-1]*quadwts[q]
- m = 0
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[numQuadPts-1][g]*basis[q][0]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][0]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][g]*basisDeriv[q][0][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[numQuadPts-1][numBasisFns-1]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][numBasisFns-1][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][numBasisFns-1]*basisDeriv[q][f][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[numQuadPts-1][g]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][g]*basisDeriv[q][numBasisFns-1][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[numQuadPts-1][0]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][0][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][0]*basisDeriv[q][f][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[numQuadPts-1][g]*basis[0][f]
- self.basisDeriv[q][numQuadPts-1][f][g][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][f]
- self.basisDeriv[q][numQuadPts-1][f][g][1] = basis[numQuadPts-1][g]*basisDeriv[q][f][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 2D function tabulation')
- n += 1
- # Top
- for r in range(numQuadPts-1, 0, -1):
- self.quadPts[n][0] = quadpts[r]
- self.quadPts[n][1] = quadpts[numQuadPts-1]
- self.quadWts[n] = quadwts[r]*quadwts[numQuadPts-1]
- m = 0
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[numQuadPts-1][0]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][0]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[numQuadPts-1][0][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[numQuadPts-1][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[numQuadPts-1][f]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[numQuadPts-1][f][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[numQuadPts-1][numBasisFns-1]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][numBasisFns-1]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[numQuadPts-1][numBasisFns-1][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[numQuadPts-1][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[numQuadPts-1][f]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[numQuadPts-1][f][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[numQuadPts-1][f]
- self.basisDeriv[numQuadPts-1][r][f][g][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][f]
- self.basisDeriv[numQuadPts-1][r][f][g][1] = basis[r][g]*basisDeriv[numQuadPts-1][f][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 2D function tabulation')
- n += 1
- # Left
- for q in range(numQuadPts-1, 0, -1):
- self.quadPts[n][0] = quadpts[0]
- self.quadPts[n][1] = quadpts[q]
- self.quadWts[n] = quadwts[0]*quadwts[q]
- m = 0
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[0][g]*basis[q][0]
- self.basisDeriv[n][m][0] = basisDeriv[0][g][0]*basis[q][0]
- self.basisDeriv[n][m][1] = basis[0][g]*basisDeriv[q][0][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[0][numBasisFns-1]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[0][numBasisFns-1][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[0][numBasisFns-1]*basisDeriv[q][f][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[0][g]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][0] = basisDeriv[0][g][0]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][1] = basis[0][g]*basisDeriv[q][numBasisFns-1][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[0][0]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[0][0][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[0][0]*basisDeriv[q][f][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[0][g]*basis[0][f]
- self.basisDeriv[q][0][f][g][0] = basisDeriv[0][g][0]*basis[q][f]
- self.basisDeriv[q][0][f][g][1] = basis[0][g]*basisDeriv[q][f][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 2D function tabulation')
- n += 1
- # Interior
- for q in range(1, numQuadPts-1):
- for r in range(1, numQuadPts-1):
- self.quadPts[n][0] = quadpts[r]
+ for q in range(0, numQuadPts):
+ for p in range(0, numQuadPts):
+ self.quadPts[n][0] = quadpts[p]
self.quadPts[n][1] = quadpts[q]
- self.quadWts[n] = quadwts[r]*quadwts[q]
+ self.quadWts[n] = quadwts[p]*quadwts[q]
+
m = 0
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[q][0]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[q][0]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[q][0][0]
+ for (bp,bq) in vertexOrder:
+ self.basis[n][m] = basis[p][bp]*basis[q][bq]
+ self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]
+ self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]
m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[q][f][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[q][numBasisFns-1][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[q][f][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[q][f]
- self.basisDeriv[q][r][f][g][0] = basisDeriv[r][g][0]*basis[q][f]
- self.basisDeriv[q][r][f][g][1] = basis[r][g]*basisDeriv[q][f][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 2D function tabulation')
+ if not m == self.numCorners:
+ raise RuntimeError('Invalid 2-D basis tabulation: '+str(m)+ \
+ ' should be '+str(self.numCorners))
n += 1
- if not n == self.numQuadPts: raise RuntimeError('Invalid 2D quadrature')
+ if not n == self.numQuadPts:
+ raise RuntimeError('Invalid 2-D quadrature: '+str(n)+ \
+ ' should be '+str(self.numQuadPts))
+
elif dim == 3:
+ # Set order of vertices and basis functions.
+ # Corners
+ vertexOrder = [(0,0,0), (1,0,0), (1,1,0), (0,1,0),
+ (0,0,1), (1,0,1), (1,1,1), (0,1,1)]
+ # Edges
+ # Bottom front
+ p = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ r = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+ # Bottom right
+ p = numpy.ones(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(2, numBasis, dtype=numpy.int32)
+ r = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+ # Bottom back
+ p = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ q = numpy.ones(numBasis-2, dtype=numpy.int32)
+ r = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+ # Bottom left
+ p = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ r = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+ # Top front
+ p = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ r = numpy.ones(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+ # Top right
+ p = numpy.ones(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(2, numBasis, dtype=numpy.int32)
+ r = numpy.ones(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+ # Top back
+ p = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ q = numpy.ones(numBasis-2, dtype=numpy.int32)
+ r = numpy.ones(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+ # Top left
+ p = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ r = numpy.ones(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+ # Middle left front
+ p = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ q = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ r = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+ # Middle right front
+ p = numpy.ones(numBasis-2, dtype=numpy.int32)
+ q = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ r = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+ # Middle right back
+ p = numpy.ones(numBasis-2, dtype=numpy.int32)
+ q = numpy.ones(numBasis-2, dtype=numpy.int32)
+ r = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+ # Middle left back
+ p = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ q = numpy.ones(numBasis-2, dtype=numpy.int32)
+ r = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
+
+ # Face
+ # Left / Right
+ ip = numpy.arange(0, 2, dtype=numpy.int32)
+ p = numpy.tile(ip, ((numBasis-2)*(numBasis-2), 1)).transpose()
+ iq = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.tile(iq, (1, 2*(numBasis-2)))
+ ir = numpy.arange(2, numBasis, dtype=numpy.int32)
+ r = numpy.tile(ir, (2, numBasis-2)).transpose()
+ vertexOrder += zip(p.ravel(),q.ravel(),r.ravel())
+ # Front / Back
+ ip = numpy.arange(2, numBasis, dtype=numpy.int32)
+ p = numpy.tile(ip, (1, 2*(numBasis-2)))
+ iq = numpy.arange(0, 2, dtype=numpy.int32)
+ q = numpy.tile(iq, ((numBasis-2)*(numBasis-2), 1)).transpose()
+ ir = numpy.arange(2, numBasis, dtype=numpy.int32)
+ r = numpy.tile(ir, (2, numBasis-2)).transpose()
+ vertexOrder += zip(p.ravel(),q.ravel(),r.ravel())
+ # Bottom / Top
+ ip = numpy.arange(2, numBasis, dtype=numpy.int32)
+ p = numpy.tile(ip, (1, 2*(numBasis-2)))
+ iq = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.tile(iq, (2, numBasis-2)).transpose()
+ ir = numpy.arange(0, 2, dtype=numpy.int32)
+ r = numpy.tile(ir, ((numBasis-2)*(numBasis-2), 1)).transpose()
+ vertexOrder += zip(p.ravel(),q.ravel(),r.ravel())
+
+ # Interior
+ ip = numpy.arange(2, numBasis, dtype=numpy.int32)
+ p = numpy.tile(ip, (1, (numBasis-2)*(numBasis-2)))
+ iq = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.tile(iq, ((numBasis-2), numBasis-2)).transpose()
+ ir = numpy.arange(2, numBasis, dtype=numpy.int32)
+ r = numpy.tile(ir, ((numBasis-2)*(numBasis-2), 1)).transpose()
+ vertexOrder += zip(p.ravel(),q.ravel(),r.ravel())
+
self.vertices = numpy.zeros((self.numCorners, dim))
n = 0
- # Depth
- for s in range(numBasisFns):
- # Bottom
- for r in range(0, numBasisFns-1):
- self.vertices[n][0] = vertices[r]
- self.vertices[n][1] = vertices[0]
- self.vertices[n][2] = vertices[s]
- n += 1
- # Right
- for q in range(0, numBasisFns-1):
- self.vertices[n][0] = vertices[numBasisFns-1]
- self.vertices[n][1] = vertices[q]
- self.vertices[n][2] = vertices[s]
- n += 1
- # Top
- for r in range(numBasisFns-1, 0, -1):
- self.vertices[n][0] = vertices[r]
- self.vertices[n][1] = vertices[numBasisFns-1]
- self.vertices[n][2] = vertices[s]
- n += 1
- # Left
- for q in range(numBasisFns-1, 0, -1):
- self.vertices[n][0] = vertices[0]
- self.vertices[n][1] = vertices[q]
- self.vertices[n][2] = vertices[s]
- n += 1
- # Interior
- for q in range(1, numBasisFns-1):
- for r in range(1, numBasisFns-1):
- self.vertices[n][0] = vertices[r]
- self.vertices[n][1] = vertices[q]
- self.vertices[n][2] = vertices[s]
- n += 1
- if not n == self.numCorners: raise RuntimeError('Invalid 3D function tabulation: '+str(n)+' should be '+str(self.numCorners))
+ for (p,q,r) in vertexOrder:
+ self.vertices[n][0] = vertices[p]
+ self.vertices[n][1] = vertices[q]
+ self.vertices[n][2] = vertices[r]
+ n += 1
+ if not n == self.numCorners:
+ raise RuntimeError('Invalid 3-D vertex ordering: '+str(n)+ \
+ ' should be '+str(self.numCorners))
self.quadPts = numpy.zeros((numQuadPts*numQuadPts*numQuadPts, dim))
self.quadWts = numpy.zeros((numQuadPts*numQuadPts*numQuadPts,))
self.basis = numpy.zeros((numQuadPts*numQuadPts*numQuadPts,
- numBasisFns*numBasisFns*numBasisFns))
+ numBasis*numBasis*numBasis))
self.basisDeriv = numpy.zeros((numQuadPts*numQuadPts*numQuadPts,
- numBasisFns*numBasisFns*numBasisFns,
+ numBasis*numBasis*numBasis,
dim))
+
+ # Order of quadrature points doesn't matter
+ # Order of basis functions should match vertices for isoparametric
n = 0
- # Depth
- for s in range(numQuadPts):
- # Bottom
- for r in range(0, numQuadPts-1):
- self.quadPts[n][0] = quadpts[r]
- self.quadPts[n][1] = quadpts[0]
- self.quadPts[n][2] = quadpts[s]
- self.quadWts[n] = quadwts[r]*quadwts[0]*quadwts[s]
- m = 0
- for h in range(numBasisFns):
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[0][0]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[0][0]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[0][0][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[0][0]*basisDeriv[s][h][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[0][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[0][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[0][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][numBasisFns-1]*basis[0][f]*basisDeriv[s][h][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[0][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[0][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[0][numBasisFns-1][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[0][numBasisFns-1]*basisDeriv[s][h][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[0][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[0][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[0][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][0]*basis[0][f]*basisDeriv[s][h][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[0][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[0][f]*basis[s][h]
- self.basisDeriv[m][m][1] = basis[r][g]*basisDeriv[0][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[0][f]*basisDeriv[s][h][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 3D function tabulation')
- n += 1
- # Right
- for q in range(0, numQuadPts-1):
- self.quadPts[n][0] = quadpts[numQuadPts-1]
- self.quadPts[n][1] = quadpts[q]
- self.quadPts[n][2] = quadpts[s]
- self.quadWts[n] = quadwts[numQuadPts-1]*quadwts[q]*quadwts[s]
- m = 0
- for h in range(numBasisFns):
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[numQuadPts-1][g]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][g]*basisDeriv[q][0][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[numQuadPts-1][g]*basis[q][0]*basisDeriv[s][h][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[numQuadPts-1][numBasisFns-1]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][numBasisFns-1][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][numBasisFns-1]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[numQuadPts-1][numBasisFns-1]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[numQuadPts-1][g]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][g]*basisDeriv[q][numBasisFns-1][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[numQuadPts-1][g]*basis[q][numBasisFns-1]*basisDeriv[s][h][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[numQuadPts-1][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][0][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][0]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[numQuadPts-1][0]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[numQuadPts-1][g]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[m][m][1] = basis[numQuadPts-1][g]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[numQuadPts-1][g]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 3D function tabulation')
- n += 1
- # Top
- for r in range(numQuadPts-1, 0, -1):
- self.quadPts[n][0] = quadpts[r]
- self.quadPts[n][1] = quadpts[numQuadPts-1]
- self.quadPts[n][2] = quadpts[s]
- self.quadWts[n] = quadwts[r]*quadwts[numQuadPts-1]*quadwts[s]
- m = 0
- for h in range(numBasisFns):
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[numQuadPts-1][0]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][0]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[numQuadPts-1][0][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[numQuadPts-1][0]*basisDeriv[s][h][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[numQuadPts-1][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][numBasisFns-1]*basis[numQuadPts-1][f]*basisDeriv[s][h][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[numQuadPts-1][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[numQuadPts-1][numBasisFns-1][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[numQuadPts-1][numBasisFns-1]*basisDeriv[s][h][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[numQuadPts-1][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][0]*basis[numQuadPts-1][f]*basisDeriv[s][h][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[m][m][1] = basis[r][g]*basisDeriv[numQuadPts-1][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[numQuadPts-1][f]*basisDeriv[s][h][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 3D function tabulation')
- n += 1
- # Left
- for q in range(numQuadPts-1, 0, -1):
- self.quadPts[n][0] = quadpts[0]
- self.quadPts[n][1] = quadpts[q]
- self.quadPts[n][2] = quadpts[s]
- self.quadWts[n] = quadwts[0]*quadwts[q]*quadwts[s]
- m = 0
- for h in range(numBasisFns):
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[0][g]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[0][g][0]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[0][g]*basisDeriv[q][0][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[0][g]*basis[q][0]*basisDeriv[s][h][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[0][numBasisFns-1]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[0][numBasisFns-1][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[0][numBasisFns-1]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[0][numBasisFns-1]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[0][g]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[0][g][0]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[0][g]*basisDeriv[q][numBasisFns-1][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[0][g]*basis[q][numBasisFns-1]*basisDeriv[s][h][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[0][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[0][0][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[0][0]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[0][0]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[0][g]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[0][g][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[m][m][1] = basis[0][g]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[0][g]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 3D function tabulation')
- n += 1
- # Interior
- for q in range(1, numQuadPts-1):
- for r in range(1, numQuadPts-1):
- self.quadPts[n][0] = quadpts[r]
+ for r in range(0, numQuadPts):
+ for q in range(0, numQuadPts):
+ for p in range(0, numQuadPts):
+ self.quadPts[n][0] = quadpts[p]
self.quadPts[n][1] = quadpts[q]
- self.quadPts[n][2] = quadpts[s]
- self.quadWts[n] = quadwts[r]*quadwts[q]*quadwts[s]
+ self.quadPts[n][2] = quadpts[r]
+ self.quadWts[n] = quadwts[p]*quadwts[q]*quadwts[r]
+
m = 0
- for h in range(numBasisFns):
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[q][0][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[q][0]*basisDeriv[s][h][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][numBasisFns-1]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[q][numBasisFns-1][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[q][numBasisFns-1]*basisDeriv[s][h][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][0]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[m][m][1] = basis[r][g]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 3D function tabulation')
+ for (bp,bq,br) in vertexOrder:
+ self.basis[n][m] = basis[p][bp]*basis[q][bq]*basis[r][br]
+ self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]*basis[r][br]
+ self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]*basis[r][br]
+ self.basisDeriv[n][m][2] = basis[p][bp]*basis[q][bq]*basisDeriv[r][br][0]
+ m += 1
+
+ if not m == self.numCorners:
+ raise RuntimeError('Invalid 3-D basis tabulation: '+str(m)+ \
+ ' should be '+str(self.numCorners))
n += 1
- if not n == self.numQuadPts: raise RuntimeError('Invalid 3D quadrature')
+ if not n == self.numQuadPts:
+ raise RuntimeError('Invalid 3-D quadrature: '+str(n)+ \
+ ' should be '+str(self.numQuadPts))
+
self.vertices = numpy.reshape(self.vertices, (self.numCorners, dim))
self.quadPts = numpy.reshape(self.quadPts, (self.numQuadPts, dim))
self.quadWts = numpy.reshape(self.quadWts, (self.numQuadPts))
Modified: short/3D/PyLith/trunk/pylith/feassemble/FIATSimplex.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/FIATSimplex.py 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/pylith/feassemble/FIATSimplex.py 2010-09-07 22:59:05 UTC (rev 17178)
@@ -93,42 +93,49 @@
"""
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]])
+ # 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)
+
+ # No permutation in order of quadrature points
+ self.quadPts = numpy.array(quadrature.get_points())
+ self.quadWts = numpy.array(quadrature.get_weights())
+
+
self._info.line("Cell geometry: ")
self._info.line(self.geometry)
self._info.line("Vertices: ")
@@ -159,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.
@@ -193,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
@@ -211,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/trunk/pylith/perf/Mesh.py
===================================================================
--- short/3D/PyLith/trunk/pylith/perf/Mesh.py 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/pylith/perf/Mesh.py 2010-09-07 22:59:05 UTC (rev 17178)
@@ -26,9 +26,13 @@
cellTypes = {(1,2): 'line2',
(1,3): 'line3',
(2,3): 'tri3',
+ (2,6): 'tri6',
(2,4): 'quad4',
+ (2,9): 'quad9',
(3,4): 'tet4',
- (3,8): 'hex8'
+ (3,8): 'hex8',
+ (3,10): 'tet10',
+ (3,27): 'hex27'
}
"""
Modified: short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py
===================================================================
--- short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py 2010-09-07 22:59:05 UTC (rev 17178)
@@ -81,14 +81,28 @@
cellId = cellTypes[0]
if numpy.sum(cellTypes-cellId) != 0:
raise ValueError("Expecting cells to all be the same type.")
- if cellId == 5: # tri3
- ncorners = 3
- elif cellId == 12: # tri3
- ncorners = 8
+ if cellId == 1: # vertex
+ ncorners = 1
elif cellId == 3: # line2
ncorners = 2
+ elif cellId == 5: # tri3
+ ncorners = 3
elif cellId == 9: # quad4
ncorners = 4
+ elif cellId == 10: # tet4
+ ncorners = 4
+ elif cellId == 12: # hex8
+ ncorners = 8
+ elif cellId == 21: # line3
+ ncorners = 3
+ elif cellId == 22: # tri6
+ ncorners = 6
+ elif cellId == 22: # quad9
+ ncorners = 9
+ elif cellId == 24: # tet10
+ ncorners = 10
+ elif cellId == 29: # hex27
+ ncorners = 27
elif cellId == 255: # unknown?
ncorners = 1
else:
Modified: short/3D/PyLith/trunk/tests/3d/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/3d/Makefile.am 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/tests/3d/Makefile.am 2010-09-07 22:59:05 UTC (rev 17178)
@@ -18,6 +18,7 @@
SUBDIRS = \
matprops \
+ plasticity \
slipdir
# End of file
Copied: short/3D/PyLith/trunk/tests/3d/plasticity (from rev 17177, short/3D/PyLith/branches/v1.5-stable/tests/3d/plasticity)
Copied: short/3D/PyLith/trunk/tests_auto/2d/quad9 (from rev 17177, short/3D/PyLith/branches/v1.5-stable/tests_auto/2d/quad9)
Copied: short/3D/PyLith/trunk/tests_auto/2d/tri6 (from rev 17177, short/3D/PyLith/branches/v1.5-stable/tests_auto/2d/tri6)
Copied: short/3D/PyLith/trunk/tests_auto/3dnew/hex27 (from rev 17177, short/3D/PyLith/branches/v1.5-stable/tests_auto/3dnew/hex27)
Copied: short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axial_disp.spatialdb (from rev 17177, short/3D/PyLith/branches/v1.5-stable/tests_auto/3dnew/hex8/axial_disp.spatialdb)
===================================================================
--- short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axial_disp.spatialdb (rev 0)
+++ short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axial_disp.spatialdb 2010-09-07 22:59:05 UTC (rev 17178)
@@ -0,0 +1,113 @@
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 3
+ value-names = displacement-x displacement-y displacement-z
+ value-units = m m m
+ num-locs = 100
+ data-dim = 3
+ space-dim = 3
+ cs-data = cartesian {
+ to-meters = 1
+ space-dim = 3
+}
+}
+ -4.000000e+03 -4.000000e+03 -6.000000e+03 -7.037037e-01 1.851852e-01 2.777778e-01
+ -2.000000e+03 -4.000000e+03 -6.000000e+03 -3.518519e-01 1.851852e-01 2.777778e-01
+ 0.000000e+00 -4.000000e+03 -6.000000e+03 0.000000e+00 1.851852e-01 2.777778e-01
+ 2.000000e+03 -4.000000e+03 -6.000000e+03 3.518519e-01 1.851852e-01 2.777778e-01
+ 4.000000e+03 -4.000000e+03 -6.000000e+03 7.037037e-01 1.851852e-01 2.777778e-01
+ -4.000000e+03 -2.000000e+03 -6.000000e+03 -7.037037e-01 9.259259e-02 2.777778e-01
+ -2.000000e+03 -2.000000e+03 -6.000000e+03 -3.518519e-01 9.259259e-02 2.777778e-01
+ 0.000000e+00 -2.000000e+03 -6.000000e+03 0.000000e+00 9.259259e-02 2.777778e-01
+ 2.000000e+03 -2.000000e+03 -6.000000e+03 3.518519e-01 9.259259e-02 2.777778e-01
+ 4.000000e+03 -2.000000e+03 -6.000000e+03 7.037037e-01 9.259259e-02 2.777778e-01
+ -4.000000e+03 0.000000e+00 -6.000000e+03 -7.037037e-01 -0.000000e+00 2.777778e-01
+ -2.000000e+03 0.000000e+00 -6.000000e+03 -3.518519e-01 -0.000000e+00 2.777778e-01
+ 0.000000e+00 0.000000e+00 -6.000000e+03 0.000000e+00 0.000000e+00 2.777778e-01
+ 2.000000e+03 0.000000e+00 -6.000000e+03 3.518519e-01 0.000000e+00 2.777778e-01
+ 4.000000e+03 0.000000e+00 -6.000000e+03 7.037037e-01 0.000000e+00 2.777778e-01
+ -4.000000e+03 2.000000e+03 -6.000000e+03 -7.037037e-01 -9.259259e-02 2.777778e-01
+ -2.000000e+03 2.000000e+03 -6.000000e+03 -3.518519e-01 -9.259259e-02 2.777778e-01
+ 0.000000e+00 2.000000e+03 -6.000000e+03 0.000000e+00 -9.259259e-02 2.777778e-01
+ 2.000000e+03 2.000000e+03 -6.000000e+03 3.518519e-01 -9.259259e-02 2.777778e-01
+ 4.000000e+03 2.000000e+03 -6.000000e+03 7.037037e-01 -9.259259e-02 2.777778e-01
+ -4.000000e+03 4.000000e+03 -6.000000e+03 -7.037037e-01 -1.851852e-01 2.777778e-01
+ -2.000000e+03 4.000000e+03 -6.000000e+03 -3.518519e-01 -1.851852e-01 2.777778e-01
+ 0.000000e+00 4.000000e+03 -6.000000e+03 0.000000e+00 -1.851852e-01 2.777778e-01
+ 2.000000e+03 4.000000e+03 -6.000000e+03 3.518519e-01 -1.851852e-01 2.777778e-01
+ 4.000000e+03 4.000000e+03 -6.000000e+03 7.037037e-01 -1.851852e-01 2.777778e-01
+ -4.000000e+03 -4.000000e+03 -4.000000e+03 -7.037037e-01 1.851852e-01 1.851852e-01
+ -2.000000e+03 -4.000000e+03 -4.000000e+03 -3.518519e-01 1.851852e-01 1.851852e-01
+ 0.000000e+00 -4.000000e+03 -4.000000e+03 0.000000e+00 1.851852e-01 1.851852e-01
+ 2.000000e+03 -4.000000e+03 -4.000000e+03 3.518519e-01 1.851852e-01 1.851852e-01
+ 4.000000e+03 -4.000000e+03 -4.000000e+03 7.037037e-01 1.851852e-01 1.851852e-01
+ -4.000000e+03 -2.000000e+03 -4.000000e+03 -7.037037e-01 9.259259e-02 1.851852e-01
+ -2.000000e+03 -2.000000e+03 -4.000000e+03 -3.518519e-01 9.259259e-02 1.851852e-01
+ 0.000000e+00 -2.000000e+03 -4.000000e+03 0.000000e+00 9.259259e-02 1.851852e-01
+ 2.000000e+03 -2.000000e+03 -4.000000e+03 3.518519e-01 9.259259e-02 1.851852e-01
+ 4.000000e+03 -2.000000e+03 -4.000000e+03 7.037037e-01 9.259259e-02 1.851852e-01
+ -4.000000e+03 0.000000e+00 -4.000000e+03 -7.037037e-01 -0.000000e+00 1.851852e-01
+ -2.000000e+03 0.000000e+00 -4.000000e+03 -3.518519e-01 -0.000000e+00 1.851852e-01
+ 0.000000e+00 0.000000e+00 -4.000000e+03 0.000000e+00 0.000000e+00 1.851852e-01
+ 2.000000e+03 0.000000e+00 -4.000000e+03 3.518519e-01 0.000000e+00 1.851852e-01
+ 4.000000e+03 0.000000e+00 -4.000000e+03 7.037037e-01 0.000000e+00 1.851852e-01
+ -4.000000e+03 2.000000e+03 -4.000000e+03 -7.037037e-01 -9.259259e-02 1.851852e-01
+ -2.000000e+03 2.000000e+03 -4.000000e+03 -3.518519e-01 -9.259259e-02 1.851852e-01
+ 0.000000e+00 2.000000e+03 -4.000000e+03 0.000000e+00 -9.259259e-02 1.851852e-01
+ 2.000000e+03 2.000000e+03 -4.000000e+03 3.518519e-01 -9.259259e-02 1.851852e-01
+ 4.000000e+03 2.000000e+03 -4.000000e+03 7.037037e-01 -9.259259e-02 1.851852e-01
+ -4.000000e+03 4.000000e+03 -4.000000e+03 -7.037037e-01 -1.851852e-01 1.851852e-01
+ -2.000000e+03 4.000000e+03 -4.000000e+03 -3.518519e-01 -1.851852e-01 1.851852e-01
+ 0.000000e+00 4.000000e+03 -4.000000e+03 0.000000e+00 -1.851852e-01 1.851852e-01
+ 2.000000e+03 4.000000e+03 -4.000000e+03 3.518519e-01 -1.851852e-01 1.851852e-01
+ 4.000000e+03 4.000000e+03 -4.000000e+03 7.037037e-01 -1.851852e-01 1.851852e-01
+ -4.000000e+03 -4.000000e+03 -2.000000e+03 -7.037037e-01 1.851852e-01 9.259259e-02
+ -2.000000e+03 -4.000000e+03 -2.000000e+03 -3.518519e-01 1.851852e-01 9.259259e-02
+ 0.000000e+00 -4.000000e+03 -2.000000e+03 0.000000e+00 1.851852e-01 9.259259e-02
+ 2.000000e+03 -4.000000e+03 -2.000000e+03 3.518519e-01 1.851852e-01 9.259259e-02
+ 4.000000e+03 -4.000000e+03 -2.000000e+03 7.037037e-01 1.851852e-01 9.259259e-02
+ -4.000000e+03 -2.000000e+03 -2.000000e+03 -7.037037e-01 9.259259e-02 9.259259e-02
+ -2.000000e+03 -2.000000e+03 -2.000000e+03 -3.518519e-01 9.259259e-02 9.259259e-02
+ 0.000000e+00 -2.000000e+03 -2.000000e+03 0.000000e+00 9.259259e-02 9.259259e-02
+ 2.000000e+03 -2.000000e+03 -2.000000e+03 3.518519e-01 9.259259e-02 9.259259e-02
+ 4.000000e+03 -2.000000e+03 -2.000000e+03 7.037037e-01 9.259259e-02 9.259259e-02
+ -4.000000e+03 0.000000e+00 -2.000000e+03 -7.037037e-01 -0.000000e+00 9.259259e-02
+ -2.000000e+03 0.000000e+00 -2.000000e+03 -3.518519e-01 -0.000000e+00 9.259259e-02
+ 0.000000e+00 0.000000e+00 -2.000000e+03 0.000000e+00 0.000000e+00 9.259259e-02
+ 2.000000e+03 0.000000e+00 -2.000000e+03 3.518519e-01 0.000000e+00 9.259259e-02
+ 4.000000e+03 0.000000e+00 -2.000000e+03 7.037037e-01 0.000000e+00 9.259259e-02
+ -4.000000e+03 2.000000e+03 -2.000000e+03 -7.037037e-01 -9.259259e-02 9.259259e-02
+ -2.000000e+03 2.000000e+03 -2.000000e+03 -3.518519e-01 -9.259259e-02 9.259259e-02
+ 0.000000e+00 2.000000e+03 -2.000000e+03 0.000000e+00 -9.259259e-02 9.259259e-02
+ 2.000000e+03 2.000000e+03 -2.000000e+03 3.518519e-01 -9.259259e-02 9.259259e-02
+ 4.000000e+03 2.000000e+03 -2.000000e+03 7.037037e-01 -9.259259e-02 9.259259e-02
+ -4.000000e+03 4.000000e+03 -2.000000e+03 -7.037037e-01 -1.851852e-01 9.259259e-02
+ -2.000000e+03 4.000000e+03 -2.000000e+03 -3.518519e-01 -1.851852e-01 9.259259e-02
+ 0.000000e+00 4.000000e+03 -2.000000e+03 0.000000e+00 -1.851852e-01 9.259259e-02
+ 2.000000e+03 4.000000e+03 -2.000000e+03 3.518519e-01 -1.851852e-01 9.259259e-02
+ 4.000000e+03 4.000000e+03 -2.000000e+03 7.037037e-01 -1.851852e-01 9.259259e-02
+ -4.000000e+03 -4.000000e+03 0.000000e+00 -7.037037e-01 1.851852e-01 -0.000000e+00
+ -2.000000e+03 -4.000000e+03 0.000000e+00 -3.518519e-01 1.851852e-01 -0.000000e+00
+ 0.000000e+00 -4.000000e+03 0.000000e+00 0.000000e+00 1.851852e-01 0.000000e+00
+ 2.000000e+03 -4.000000e+03 0.000000e+00 3.518519e-01 1.851852e-01 0.000000e+00
+ 4.000000e+03 -4.000000e+03 0.000000e+00 7.037037e-01 1.851852e-01 0.000000e+00
+ -4.000000e+03 -2.000000e+03 0.000000e+00 -7.037037e-01 9.259259e-02 -0.000000e+00
+ -2.000000e+03 -2.000000e+03 0.000000e+00 -3.518519e-01 9.259259e-02 -0.000000e+00
+ 0.000000e+00 -2.000000e+03 0.000000e+00 0.000000e+00 9.259259e-02 0.000000e+00
+ 2.000000e+03 -2.000000e+03 0.000000e+00 3.518519e-01 9.259259e-02 0.000000e+00
+ 4.000000e+03 -2.000000e+03 0.000000e+00 7.037037e-01 9.259259e-02 0.000000e+00
+ -4.000000e+03 0.000000e+00 0.000000e+00 -7.037037e-01 0.000000e+00 0.000000e+00
+ -2.000000e+03 0.000000e+00 0.000000e+00 -3.518519e-01 0.000000e+00 0.000000e+00
+ 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
+ 2.000000e+03 0.000000e+00 0.000000e+00 3.518519e-01 0.000000e+00 0.000000e+00
+ 4.000000e+03 0.000000e+00 0.000000e+00 7.037037e-01 0.000000e+00 0.000000e+00
+ -4.000000e+03 2.000000e+03 0.000000e+00 -7.037037e-01 -9.259259e-02 0.000000e+00
+ -2.000000e+03 2.000000e+03 0.000000e+00 -3.518519e-01 -9.259259e-02 0.000000e+00
+ 0.000000e+00 2.000000e+03 0.000000e+00 0.000000e+00 -9.259259e-02 0.000000e+00
+ 2.000000e+03 2.000000e+03 0.000000e+00 3.518519e-01 -9.259259e-02 0.000000e+00
+ 4.000000e+03 2.000000e+03 0.000000e+00 7.037037e-01 -9.259259e-02 0.000000e+00
+ -4.000000e+03 4.000000e+03 0.000000e+00 -7.037037e-01 -1.851852e-01 0.000000e+00
+ -2.000000e+03 4.000000e+03 0.000000e+00 -3.518519e-01 -1.851852e-01 0.000000e+00
+ 0.000000e+00 4.000000e+03 0.000000e+00 0.000000e+00 -1.851852e-01 0.000000e+00
+ 2.000000e+03 4.000000e+03 0.000000e+00 3.518519e-01 -1.851852e-01 0.000000e+00
+ 4.000000e+03 4.000000e+03 0.000000e+00 7.037037e-01 -1.851852e-01 0.000000e+00
Copied: short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp.cfg (from rev 17177, short/3D/PyLith/branches/v1.5-stable/tests_auto/3dnew/hex8/axialdisp.cfg)
===================================================================
--- short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp.cfg (rev 0)
+++ short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp.cfg 2010-09-07 22:59:05 UTC (rev 17178)
@@ -0,0 +1,134 @@
+# -*- Python -*-
+[pylithapp]
+
+[pylithapp.launcher] # WARNING: THIS IS NOT PORTABLE
+command = mpirun -np ${nodes}
+
+# ----------------------------------------------------------------------
+# journal
+# ----------------------------------------------------------------------
+[pylithapp.journal.info]
+#timedependent = 1
+#implicit = 1
+#petsc = 1
+#solverlinear = 1
+#meshiocubit = 1
+#implicitelasticity = 1
+#quadrature3d = 1
+#fiatlagrange = 1
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+[pylithapp.mesh_generator]
+#debug = 1
+reader = pylith.meshio.MeshIOCubit
+
+[pylithapp.mesh_generator.reader]
+filename = mesh.exo
+use_nodeset_names = True
+coordsys.space_dim = 3
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[pylithapp.timedependent]
+dimension = 3
+bc = [x_neg,x_pos,y_neg,z_neg]
+
+[pylithapp.timedependent.formulation]
+output = [domain,subdomain]
+output.subdomain = pylith.meshio.OutputSolnSubset
+
+[pylithapp.timedependent.formulation.time_step]
+total_time = 0.0*s
+
+# ----------------------------------------------------------------------
+# materials
+# ----------------------------------------------------------------------
+[pylithapp.timedependent]
+materials = [elastic,viscoelastic]
+materials.elastic = pylith.materials.ElasticIsotropic3D
+materials.viscoelastic = pylith.materials.ElasticIsotropic3D
+
+[pylithapp.timedependent.materials.elastic]
+label = Elastic material
+id = 1
+db_properties.iohandler.filename = matprops.spatialdb
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 3
+
+[pylithapp.timedependent.materials.viscoelastic]
+label = Elastic material 2
+id = 2
+db_properties.iohandler.filename = matprops.spatialdb
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 3
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+[pylithapp.timedependent.bc.x_pos]
+bc_dof = [0]
+label = face_xpos
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC +x eface
+db_initial.iohandler.filename = axial_disp.spatialdb
+
+[pylithapp.timedependent.bc.x_neg]
+bc_dof = [0]
+label = face_xneg
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC -x face
+db_initial.iohandler.filename = axial_disp.spatialdb
+
+[pylithapp.timedependent.bc.y_neg]
+bc_dof = [1]
+label = face_yneg
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC -y face
+db_initial.iohandler.filename = axial_disp.spatialdb
+
+[pylithapp.timedependent.bc.z_neg]
+bc_dof = [2]
+label = face_zneg
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC -z edge
+db_initial.iohandler.filename = axial_disp.spatialdb
+
+# ----------------------------------------------------------------------
+# PETSc
+# ----------------------------------------------------------------------
+[pylithapp.petsc]
+pc_type = asm
+
+# Change the preconditioner settings.
+sub_pc_factor_shift_type = none
+
+ksp_rtol = 1.0e-8
+ksp_max_it = 100
+ksp_gmres_restart = 50
+ksp_monitor = true
+ksp_view = true
+ksp_converged_reason = true
+#log_summary = true
+# start_in_debugger = true
+
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[pylithapp.problem.formulation.output.domain.writer]
+filename = axialdisp.vtk
+
+[pylithapp.problem.formulation.output.subdomain]
+label = face_zpos
+writer.filename = axialdisp-groundsurf.vtk
+
+[pylithapp.timedependent.materials.elastic.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = axialdisp-elastic.vtk
+
+[pylithapp.timedependent.materials.viscoelastic.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = axialdisp-viscoelastic.vtk
Copied: short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp_gendb.py (from rev 17177, short/3D/PyLith/branches/v1.5-stable/tests_auto/3dnew/hex8/axialdisp_gendb.py)
===================================================================
--- short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp_gendb.py (rev 0)
+++ short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp_gendb.py 2010-09-07 22:59:05 UTC (rev 17178)
@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/3d/hex8/axialdisp_gendb.py
+##
+## @brief Python script to generate spatial database with displacement
+## boundary conditions for the axial displacement test.
+
+import numpy
+
+class GenerateDB(object):
+ """
+ Python object to generate spatial database with displacement
+ boundary conditions for the axial displacement test.
+ """
+
+ def __init__(self):
+ """
+ Constructor.
+ """
+ return
+
+
+ def run(self):
+ """
+ Generate the database.
+ """
+ # Domain
+ x = numpy.arange(-4000.0, 4000.1, 1000.0)
+ y = numpy.arange(-4000.0, 4000.1, 1000.0)
+ z = numpy.arange(-6000.0, 0.1, 1000.0)
+ nxpts = x.shape[0]
+ nypts = y.shape[0]
+ nzpts = z.shape[0]
+
+ xx = numpy.tile(x, (1,nypts*nzpts))
+ yy = numpy.tile(y, (nxpts,nzpts)).transpose()
+ zz = numpy.tile(z, (nxpts*nypts,1)).transpose()
+
+ xyz = numpy.zeros( (nxpts*nypts*nzpts, 3), dtype=numpy.float64)
+ xyz[:,0] = numpy.ravel(xx)
+ xyz[:,1] = numpy.ravel(yy)
+ xyz[:,2] = numpy.ravel(zz)
+
+ from axialdisp_soln import AnalyticalSoln
+ soln = AnalyticalSoln()
+ disp = soln.displacement(xyz)
+
+ from spatialdata.geocoords.CSCart import CSCart
+ cs = CSCart()
+ cs.inventory.spaceDim = 3
+ cs._configure()
+ data = {'points': xyz,
+ 'coordsys': cs,
+ 'data_dim': 3,
+ 'values': [{'name': "displacement-x",
+ 'units': "m",
+ 'data': numpy.ravel(disp[:,0])},
+ {'name': "displacement-y",
+ 'units': "m",
+ 'data': numpy.ravel(disp[:,1])},
+ {'name': "displacement-z",
+ 'units': "m",
+ 'data': numpy.ravel(disp[:,2])}]}
+
+ from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
+ io = SimpleIOAscii()
+ io.inventory.filename = "axial_disp.spatialdb"
+ io._configure()
+ io.write(data)
+ return
+
+
+# ======================================================================
+if __name__ == "__main__":
+ app = GenerateDB()
+ app.run()
+
+
+# End of file
Copied: short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp_soln.py (from rev 17177, short/3D/PyLith/branches/v1.5-stable/tests_auto/3dnew/hex8/axialdisp_soln.py)
===================================================================
--- short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp_soln.py (rev 0)
+++ short/3D/PyLith/trunk/tests_auto/3dnew/hex8/axialdisp_soln.py 2010-09-07 22:59:05 UTC (rev 17178)
@@ -0,0 +1,107 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/axialdisp_soln.py
+##
+## @brief Analytical solution to axial displacement problem.
+
+## 3-D axial extension test with linear hexahedral cells.
+
+import numpy
+
+# Physical properties
+p_density = 2500.0
+p_vs = 3000.0
+p_vp = 5291.502622129181
+
+p_mu = p_density*p_vs**2
+p_lambda = p_density*p_vp**2 - 2*p_mu
+
+# Uniform stress field (plane strain)
+sxx = 1.0e+7
+syy = 0.0
+szz = 0.0
+sxy = 0.0
+syz = 0.0
+sxz = 0.0
+
+# Uniform strain field
+exx = 1.0/(2*p_mu) * (sxx - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+eyy = 1.0/(2*p_mu) * (syy - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+ezz = 1.0/(2*p_mu) * (szz - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+
+exy = 1.0/(2*p_mu) * (sxy)
+eyz = 1.0/(2*p_mu) * (syz)
+exz = 1.0/(2*p_mu) * (sxz)
+
+#print exx,eyy,exy,ezz
+#print -exx*p_lambda/(p_lambda+2*p_mu)
+
+# ----------------------------------------------------------------------
+class AnalyticalSoln(object):
+ """
+ Analytical solution to axial/shear displacement problem.
+ """
+
+ def __init__(self):
+ return
+
+
+ def displacement(self, locs):
+ """
+ Compute displacement field at locations.
+ """
+ (npts, dim) = locs.shape
+ disp = numpy.zeros( (npts, 3), dtype=numpy.float64)
+ disp[:,0] = exx*locs[:,0] + exy*locs[:,1] + exz*locs[:,2]
+ disp[:,1] = eyy*locs[:,1] + exy*locs[:,0] + eyz*locs[:,2]
+ disp[:,2] = ezz*locs[:,2] + exz*locs[:,0] + eyz*locs[:,1]
+ return disp
+
+
+ def strain(self, locs):
+ """
+ Compute strain field at locations.
+ """
+ (npts, dim) = locs.shape
+ strain = numpy.zeros( (npts, 6), dtype=numpy.float64)
+ strain[:,0] = exx
+ strain[:,1] = eyy
+ strain[:,2] = ezz
+ strain[:,3] = exy
+ strain[:,4] = eyz
+ strain[:,5] = exz
+ return strain
+
+
+ def stress(self, locs):
+ """
+ Compute stress field at locations.
+ """
+ (npts, dim) = locs.shape
+ stress = numpy.zeros( (npts, 6), dtype=numpy.float64)
+ stress[:,0] = sxx
+ stress[:,1] = syy
+ stress[:,2] = szz
+ stress[:,3] = sxy
+ stress[:,4] = syz
+ stress[:,5] = sxz
+ return stress
+
+
+# End of file
Modified: short/3D/PyLith/trunk/tests_auto/3dnew/hex8/mesh.exo
===================================================================
(Binary files differ)
Modified: short/3D/PyLith/trunk/tests_auto/3dnew/hex8/mesh.jou
===================================================================
--- short/3D/PyLith/trunk/tests_auto/3dnew/hex8/mesh.jou 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/tests_auto/3dnew/hex8/mesh.jou 2010-09-07 22:59:05 UTC (rev 17178)
@@ -41,7 +41,7 @@
group "face_xpos" add node in surface 31
group "face_xpos" add node in surface 39
nodeset 20 group face_xpos
-nodeset 20 name "face xpos"
+nodeset 20 name "face_xpos"
# ----------------------------------------------------------------------
# Create nodeset for -x face
@@ -49,7 +49,7 @@
group "face_xneg" add node in surface 41
group "face_xneg" add node in surface 49
nodeset 21 group face_xneg
-nodeset 21 name "face xneg"
+nodeset 21 name "face_xneg"
# ----------------------------------------------------------------------
# Create nodeset for +y face
@@ -61,7 +61,7 @@
group "face_ypos" add node in surface 46
group "face_ypos" add node in surface 56
nodeset 22 group face_ypos
-nodeset 22 name "face ypos"
+nodeset 22 name "face_ypos"
# ----------------------------------------------------------------------
# Create nodeset for -y face
@@ -73,7 +73,7 @@
group "face_yneg" add node in surface 48
group "face_yneg" add node in surface 58
nodeset 23 group face_yneg
-nodeset 23 name "face yneg"
+nodeset 23 name "face_yneg"
# ----------------------------------------------------------------------
# Create nodeset for +z face
@@ -82,7 +82,7 @@
group "face_zpos" add node in surface 27
group "face_zpos" add node in surface 23
nodeset 24 group face_zpos
-nodeset 24 name "face zpos"
+nodeset 24 name "face_zpos"
# ----------------------------------------------------------------------
# Create nodeset for -z face
@@ -91,7 +91,7 @@
group "face_zneg" add node in surface 28
group "face_zneg" add node in surface 21
nodeset 25 group face_zneg
-nodeset 25 name "face zneg"
+nodeset 25 name "face_zneg"
# ----------------------------------------------------------------------
# Create nodeset for -z face
@@ -99,7 +99,7 @@
group "face_zneg_nofault" add node in face_zneg
group "face_zneg_nofault" remove node in fault
nodeset 26 group face_zneg_nofault
-nodeset 26 name "face zneg nofault"
+nodeset 26 name "face_zneg_nofault"
# ----------------------------------------------------------------------
# Export exodus file
Copied: short/3D/PyLith/trunk/tests_auto/3dnew/tet10 (from rev 17177, short/3D/PyLith/branches/v1.5-stable/tests_auto/3dnew/tet10)
Copied: short/3D/PyLith/trunk/tests_auto/3dnew/tet4 (from rev 17177, short/3D/PyLith/branches/v1.5-stable/tests_auto/3dnew/tet4)
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -167,52 +167,29 @@
const int numCorners = _data->numCorners;
const int spaceDim = _data->spaceDim;
const ALE::Obj<SieveSubMesh::label_sequence>& cells = submesh->heightStratum(1);
- const int numBoundaryVertices = submesh->depthStratum(0)->size();
- const int numBoundaryCells = cells->size();
+ const int numVertices = submesh->depthStratum(0)->size();
+ const int numCells = cells->size();
+ const int boundaryDepth = submesh->depth()-1; // depth of boundary cells
CPPUNIT_ASSERT_EQUAL(_data->cellDim, cellDim);
- CPPUNIT_ASSERT_EQUAL(_data->numBoundaryVertices, numBoundaryVertices);
- CPPUNIT_ASSERT_EQUAL(_data->numBoundaryCells, numBoundaryCells);
+ CPPUNIT_ASSERT_EQUAL(_data->numVertices, numVertices);
+ CPPUNIT_ASSERT_EQUAL(_data->numCells, numCells);
- const int boundaryDepth = submesh->depth()-1; // depth of boundary cells
- const ALE::Obj<SubRealSection>& coordinates =
- submesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates, numCorners*spaceDim);
- // coordinates->view("Mesh coordinates from TestNeumann::testInitialize");
-
- const int numBasis = bc._quadrature->numBasis();
- const int cellVertSize = _data->numCorners * spaceDim;
- double_array cellVertices(cellVertSize);
-
- const double tolerance = 1.0e-06;
-
- // check cell vertices
- int iCell = 0;
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = submesh->getSieve();
+ ALE::ISieveVisitor::PointRetriever<SieveSubMesh::sieve_type> pV(sieve->getMaxConeSize());
+ int dp = 0;
for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
++c_iter) {
const int numCorners = submesh->getNumCellCorners(*c_iter, boundaryDepth);
CPPUNIT_ASSERT_EQUAL(_data->numCorners, numCorners);
- coordsVisitor.clear();
- submesh->restrictClosure(*c_iter, coordsVisitor);
- double vert =0.0;
- double vertE =0.0;
- const double* cellVertices = coordsVisitor.getValues();
- // std::cout << "c_iter " << *c_iter << " vertex coords:" << std::endl;
- for(int iVert = 0; iVert < numCorners; ++iVert) {
- for(int iDim = 0; iDim < spaceDim; ++iDim) {
- vertE = _data->cellVertices[iDim+spaceDim*iVert+iCell*cellVertSize];
- vert = cellVertices[iDim+spaceDim*iVert];
- // std::cout << " " << cellVertices[iDim+spaceDim*iVert];
- if (fabs(vertE) > 1.0)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vert/vertE, tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(vert, vertE, tolerance);
- } // for
- // std::cout << std::endl;
- } // for
- iCell++;
+ sieve->cone(*c_iter, pV);
+ const SieveSubMesh::point_type *cone = pV.getPoints();
+ for(int p = 0; p < pV.getSize(); ++p, ++dp) {
+ CPPUNIT_ASSERT_EQUAL(_data->cells[dp], cone[p]);
+ }
+ pV.clear();
} // for
// Check traction values
@@ -224,6 +201,7 @@
const ALE::Obj<SubRealSection>& tractionSection =
bc._parameters->get("initial").section();
+ const double tolerance = 1.0e-06;
for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
++c_iter) {
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependent.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependent.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependent.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -125,13 +125,8 @@
PointForce bc;
bc.dbTimeHistory(&th);
- bool caught = false;
- try {
- bc.TimeDependent::verifyConfiguration(mesh);
- } catch ( const std::exception& err) {
- caught = true;
- } // catch
- CPPUNIT_ASSERT(caught);
+ CPPUNIT_ASSERT_THROW(bc.TimeDependent::verifyConfiguration(mesh),
+ std::runtime_error);
} // change (missing change)
} // testVerifyConfiguration
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataHex8.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -112,8 +112,8 @@
const int pylith::bc::AbsorbingDampersDataHex8::_numCells = 2;
const int pylith::bc::AbsorbingDampersDataHex8::_numCorners = 4;
const int pylith::bc::AbsorbingDampersDataHex8::_cells[] = {
- 4, 10, 8, 2,
- 6, 12, 10, 4,
+ 2, 4, 10, 8,
+ 4, 6, 12, 10,
};
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersDataTet4.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -35,9 +35,9 @@
0.3333333333333333,
};
const double pylith::bc::AbsorbingDampersDataTet4::_basisDerivRef[] = {
- -1.0, -1.0,
- 1.0, 0.0,
- 0.0, 1.0,
+ -0.5, -0.5,
+ 0.5, 0.0,
+ 0.0, 0.5,
};
const char* pylith::bc::AbsorbingDampersDataTet4::_spatialDBFilename =
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/BoundaryMeshDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/BoundaryMeshDataHex8.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/BoundaryMeshDataHex8.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -33,8 +33,8 @@
2, 4, 6, 8, 10, 12
};
const int pylith::bc::BoundaryMeshDataHex8::_cellsNoFault[] = {
- 4, 10, 8, 2,
- 6, 12, 10, 4
+ 2, 4, 10, 8,
+ 4, 6, 12, 10,
};
const int pylith::bc::BoundaryMeshDataHex8::_numVerticesFault = 8;
@@ -42,8 +42,8 @@
2, 4, 6, 8, 10, 12, 14, 16
};
const int pylith::bc::BoundaryMeshDataHex8::_cellsFault[] = {
- 14, 16, 8, 2,
- 6, 12, 10, 4,
+ 2, 14, 16, 8,
+ 4, 6, 12, 10,
};
pylith::bc::BoundaryMeshDataHex8::BoundaryMeshDataHex8(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -33,10 +33,10 @@
label(""),
spaceDim(0),
cellDim(0),
- numBoundaryVertices(0),
- numBoundaryCells(0),
+ numVertices(0),
+ numCells(0),
numCorners(0),
- cellVertices(0),
+ cells(0),
tractionsCell(0),
valsResidual(0)
{ // constructor
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.hh 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannData.hh 2010-09-07 22:59:05 UTC (rev 17178)
@@ -61,12 +61,12 @@
/// @name Boundary mesh information
//@{
- int spaceDim; ///< Number of dimensions of vertex coordinates
- int cellDim; ///< Dimension of surface cells.
- int numBoundaryVertices; ///< Number of boundary vertices in the mesh.
- int numBoundaryCells; ///< Number of cells on Neumann boundary.
- int numCorners; ///< Number of vertices for each boundary cell.
- double* cellVertices; ///< Vertex coordinates for boundary cells.
+ int spaceDim; ///< Number of dimensions in vertex coordinates
+ int cellDim; ///< Number of dimensions associated with cell
+ int numVertices; ///< Number of vertices
+ int numCells; ///< Number of cells
+ int numCorners; ///< Number of vertices in cell
+ int* cells; ///< Indices of vertices in cells
//@}
/// @name Calculated values.
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -78,18 +78,14 @@
const int pylith::bc::NeumannDataHex8::_spaceDim = 3;
const int pylith::bc::NeumannDataHex8::_cellDim = 2;
-
-const int pylith::bc::NeumannDataHex8::_numBoundaryVertices = 6;
-const int pylith::bc::NeumannDataHex8::_numBoundaryCells = 2;
+const int pylith::bc::NeumannDataHex8::_numVertices = 6;
+const int pylith::bc::NeumannDataHex8::_numCells = 2;
const int pylith::bc::NeumannDataHex8::_numCorners = 4;
-const double pylith::bc::NeumannDataHex8::_cellVertices[] = { 0.0,-1.0,-1.0,
- 0.0,-1.0, 1.0,
- -2.0,-1.0, 1.0,
- -2.0,-1.0,-1.0,
- 2.0,-1.0,-1.0,
- 2.0,-1.0, 1.0,
- 0.0,-1.0, 1.0,
- 0.0,-1.0,-1.0};
+const int pylith::bc::NeumannDataHex8::_cells[] = {
+ 4, 2, 6, 8,
+ 8, 6, 10, 12,
+};
+
const double pylith::bc::NeumannDataHex8::_tractionsCell[] = { 4.0, 0.0, 0.0,
4.0, 0.0, 0.0,
4.0, 0.0, 0.0,
@@ -129,11 +125,11 @@
spaceDim = _spaceDim;
cellDim = _cellDim;
- numBoundaryVertices = _numBoundaryVertices;
- numBoundaryCells = _numBoundaryCells;
+ numVertices = _numVertices;
+ numCells = _numCells;
numCorners = _numCorners;
+ cells = const_cast<int*>(_cells);
- cellVertices = const_cast<double*>(_cellVertices);
tractionsCell = const_cast<double*>(_tractionsCell);
valsResidual = const_cast<double*>(_valsResidual);
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.hh 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataHex8.hh 2010-09-07 22:59:05 UTC (rev 17178)
@@ -60,10 +60,10 @@
// Mesh information
static const int _spaceDim;
static const int _cellDim;
- static const int _numBoundaryVertices;
- static const int _numBoundaryCells;
+ static const int _numVertices;
+ static const int _numCells;
static const int _numCorners;
- static const double _cellVertices[];
+ static const int _cells[];
// Calculated values.
static const double _tractionsCell[];
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataLine2.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataLine2.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -52,11 +52,13 @@
const int pylith::bc::NeumannDataLine2::_spaceDim = 1;
const int pylith::bc::NeumannDataLine2::_cellDim = 0;
+const int pylith::bc::NeumannDataLine2::_numVertices = 1;
+const int pylith::bc::NeumannDataLine2::_numCells = 1;
+const int pylith::bc::NeumannDataLine2::_numCorners = 1;
+const int pylith::bc::NeumannDataLine2::_cells[] = {
+ 2,
+};
-const int pylith::bc::NeumannDataLine2::_numBoundaryVertices = 1;
-const int pylith::bc::NeumannDataLine2::_numBoundaryCells = 1;
-const int pylith::bc::NeumannDataLine2::_numCorners = 1;
-const double pylith::bc::NeumannDataLine2::_cellVertices[] = { -1.0};
const double pylith::bc::NeumannDataLine2::_tractionsCell[] = {
1.0,
};
@@ -84,11 +86,11 @@
spaceDim = _spaceDim;
cellDim = _cellDim;
- numBoundaryVertices = _numBoundaryVertices;
- numBoundaryCells = _numBoundaryCells;
+ numVertices = _numVertices;
+ numCells = _numCells;
numCorners = _numCorners;
+ cells = const_cast<int*>(_cells);
- cellVertices = const_cast<double*>(_cellVertices);
tractionsCell = const_cast<double*>(_tractionsCell);
valsResidual = const_cast<double*>(_valsResidual);
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataLine2.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataLine2.hh 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataLine2.hh 2010-09-07 22:59:05 UTC (rev 17178)
@@ -60,10 +60,10 @@
// Mesh information
static const int _spaceDim;
static const int _cellDim;
- static const int _numBoundaryVertices;
- static const int _numBoundaryCells;
+ static const int _numVertices;
+ static const int _numCells;
static const int _numCorners;
- static const double _cellVertices[];
+ static const int _cells[];
// Calculated values.
static const double _tractionsCell[];
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataQuad4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataQuad4.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataQuad4.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -70,14 +70,14 @@
const int pylith::bc::NeumannDataQuad4::_spaceDim = 2;
const int pylith::bc::NeumannDataQuad4::_cellDim = 1;
+const int pylith::bc::NeumannDataQuad4::_numVertices = 3;
+const int pylith::bc::NeumannDataQuad4::_numCells = 2;
+const int pylith::bc::NeumannDataQuad4::_numCorners = 2;
+const int pylith::bc::NeumannDataQuad4::_cells[] = {
+ 2, 4,
+ 4, 6,
+};
-const int pylith::bc::NeumannDataQuad4::_numBoundaryVertices = 3;
-const int pylith::bc::NeumannDataQuad4::_numBoundaryCells = 2;
-const int pylith::bc::NeumannDataQuad4::_numCorners = 2;
-const double pylith::bc::NeumannDataQuad4::_cellVertices[] = {-1.0,-1.0,
- 0.0,-1.0,
- 0.0,-1.0,
- 1.0,-1.0};
const double pylith::bc::NeumannDataQuad4::_tractionsCell[] = {
0.0, -0.1056624327,
0.0, -0.3943375673,
@@ -111,11 +111,11 @@
spaceDim = _spaceDim;
cellDim = _cellDim;
- numBoundaryVertices = _numBoundaryVertices;
- numBoundaryCells = _numBoundaryCells;
+ numVertices = _numVertices;
+ numCells = _numCells;
numCorners = _numCorners;
+ cells = const_cast<int*>(_cells);
- cellVertices = const_cast<double*>(_cellVertices);
tractionsCell = const_cast<double*>(_tractionsCell);
valsResidual = const_cast<double*>(_valsResidual);
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataQuad4.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataQuad4.hh 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataQuad4.hh 2010-09-07 22:59:05 UTC (rev 17178)
@@ -60,10 +60,10 @@
// Mesh information
static const int _spaceDim;
static const int _cellDim;
- static const int _numBoundaryVertices;
- static const int _numBoundaryCells;
+ static const int _numVertices;
+ static const int _numCells;
static const int _numCorners;
- static const double _cellVertices[];
+ static const int _cells[];
// Calculated values.
static const double _tractionsCell[];
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTet4.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTet4.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -53,9 +53,9 @@
0.3333333333333333,
};
const double pylith::bc::NeumannDataTet4::_basisDerivRef[] = {
- -1.0, -1.0,
- 1.0, 0.0,
- 0.0, 1.0,
+ -0.5, -0.5,
+ 0.5, 0.0,
+ 0.0, 0.5,
};
const char* pylith::bc::NeumannDataTet4::_spatialDBFilename =
@@ -65,13 +65,13 @@
const int pylith::bc::NeumannDataTet4::_spaceDim = 3;
const int pylith::bc::NeumannDataTet4::_cellDim = 2;
+const int pylith::bc::NeumannDataTet4::_numVertices = 3;
+const int pylith::bc::NeumannDataTet4::_numCells = 1;
+const int pylith::bc::NeumannDataTet4::_numCorners = 3;
+const int pylith::bc::NeumannDataTet4::_cells[] = {
+ 3, 4, 5,
+};
-const int pylith::bc::NeumannDataTet4::_numBoundaryVertices = 3;
-const int pylith::bc::NeumannDataTet4::_numBoundaryCells = 1;
-const int pylith::bc::NeumannDataTet4::_numCorners = 3;
-const double pylith::bc::NeumannDataTet4::_cellVertices[] = { 1.0, 0.0, 0.0,
- 0.0, 1.0, 0.0,
- 0.0, 0.0, 1.0};
const double pylith::bc::NeumannDataTet4::_tractionsCell[] = {
-0.5380048025, 0.87620875991, 1.3938468501
};
@@ -101,11 +101,11 @@
spaceDim = _spaceDim;
cellDim = _cellDim;
- numBoundaryVertices = _numBoundaryVertices;
- numBoundaryCells = _numBoundaryCells;
+ numVertices = _numVertices;
+ numCells = _numCells;
numCorners = _numCorners;
+ cells = const_cast<int*>(_cells);
- cellVertices = const_cast<double*>(_cellVertices);
tractionsCell = const_cast<double*>(_tractionsCell);
valsResidual = const_cast<double*>(_valsResidual);
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTet4.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTet4.hh 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTet4.hh 2010-09-07 22:59:05 UTC (rev 17178)
@@ -60,10 +60,10 @@
// Mesh information
static const int _spaceDim;
static const int _cellDim;
- static const int _numBoundaryVertices;
- static const int _numBoundaryCells;
+ static const int _numVertices;
+ static const int _numCells;
static const int _numCorners;
- static const double _cellVertices[];
+ static const int _cells[];
// Calculated values.
static const double _tractionsCell[];
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTri3.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTri3.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -60,12 +60,13 @@
const int pylith::bc::NeumannDataTri3::_spaceDim = 2;
const int pylith::bc::NeumannDataTri3::_cellDim = 1;
+const int pylith::bc::NeumannDataTri3::_numVertices = 2;
+const int pylith::bc::NeumannDataTri3::_numCells = 1;
+const int pylith::bc::NeumannDataTri3::_numCorners = 2;
+const int pylith::bc::NeumannDataTri3::_cells[] = {
+ 3, 5,
+};
-const int pylith::bc::NeumannDataTri3::_numBoundaryVertices = 2;
-const int pylith::bc::NeumannDataTri3::_numBoundaryCells = 1;
-const int pylith::bc::NeumannDataTri3::_numCorners = 2;
-const double pylith::bc::NeumannDataTri3::_cellVertices[] = { 0.0,-1.0,
- 1.0, 0.0};
const double pylith::bc::NeumannDataTri3::_tractionsCell[] = {
1.4142135624, 0.0,
};
@@ -94,11 +95,11 @@
spaceDim = _spaceDim;
cellDim = _cellDim;
- numBoundaryVertices = _numBoundaryVertices;
- numBoundaryCells = _numBoundaryCells;
+ numVertices = _numVertices;
+ numCells = _numCells;
numCorners = _numCorners;
+ cells = const_cast<int*>(_cells);
- cellVertices = const_cast<double*>(_cellVertices);
tractionsCell = const_cast<double*>(_tractionsCell);
valsResidual = const_cast<double*>(_valsResidual);
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTri3.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTri3.hh 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/NeumannDataTri3.hh 2010-09-07 22:59:05 UTC (rev 17178)
@@ -60,10 +60,10 @@
// Mesh information
static const int _spaceDim;
static const int _cellDim;
- static const int _numBoundaryVertices;
- static const int _numBoundaryCells;
+ static const int _numVertices;
+ static const int _numCells;
static const int _numCorners;
- static const double _cellVertices[];
+ static const int _cells[];
// Calculated values.
static const double _tractionsCell[];
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -418,7 +418,7 @@
{ // testAdjustTopologyHex8h
CohesiveDataHex8h data;
FaultCohesiveTract fault;
- _testAdjustTopology(&fault, data, false);
+ _testAdjustTopology(&fault, data, true);
} // testAdjustTopologyHex8h
// ----------------------------------------------------------------------
@@ -551,6 +551,8 @@
tolerance);
} // for
+ mesh.view("MESH");
+
// check cells
const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
CPPUNIT_ASSERT(!sieve.isNull());
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8b.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8b.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8b.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -72,7 +72,7 @@
const int pylith::faults::CohesiveDataHex8b::_cells[] = {
2, 14, 15, 3, 4, 16, 17, 5,
6, 10, 11, 7, 8, 12, 13, 9,
- 8, 9, 7, 6, 16, 17, 15, 14,
+ 7, 6, 8, 9, 15, 14, 16, 17,
};
const int pylith::faults::CohesiveDataHex8b::_materialIds[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8d.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8d.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8d.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -72,7 +72,7 @@
const int pylith::faults::CohesiveDataHex8d::_cells[] = {
4, 8, 6, 2, 5, 9, 7, 3,
16, 12, 10, 14, 17, 13, 11, 15,
- 6, 7, 9, 8, 14, 15, 17, 16,
+ 8, 6, 7, 9, 16, 14, 15, 17,
};
const int pylith::faults::CohesiveDataHex8d::_materialIds[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8e.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8e.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8e.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -72,7 +72,7 @@
const int pylith::faults::CohesiveDataHex8e::_cells[] = {
5, 9, 8, 4, 3, 7, 6, 2,
17, 13, 12, 16, 15, 11, 10, 14,
- 8, 6, 7, 9, 16, 14, 15, 17,
+ 9, 8, 6, 7, 17, 16, 14, 15,
};
const int pylith::faults::CohesiveDataHex8e::_materialIds[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8f.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8f.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8f.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -72,7 +72,7 @@
const int pylith::faults::CohesiveDataHex8f::_cells[] = {
3, 15, 17, 5, 2, 14, 16, 4,
7, 11, 13, 9, 6, 10, 12, 8,
- 6, 8, 9, 7, 14, 16, 17, 15,
+ 9, 7, 6, 8, 17, 15, 14, 16,
};
const int pylith::faults::CohesiveDataHex8f::_materialIds[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8g.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8g.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8g.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -98,8 +98,8 @@
6, 24, 25, 7, 8, 26, 27, 9,
10, 16, 17, 11, 12, 18, 19, 13,
12, 18, 19, 13, 14, 20, 21, 15,
- 12, 13, 11, 10, 24, 25, 23, 22,
- 14, 15, 13, 12, 26, 27, 25, 24,
+ 11, 10, 12, 13, 23, 22, 24, 25,
+ 13, 12, 14, 15, 25, 24, 26, 27,
};
const int pylith::faults::CohesiveDataHex8g::_materialIds[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8h.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8h.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8h.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -97,8 +97,8 @@
19, 18, 20, 21, 13, 12, 14, 15,
11, 17, 19, 13, 10, 16, 18, 12,
27, 26, 8, 9, 25, 24, 6, 7,
- 11, 10, 12, 13, 23, 22, 24, 25,
- 12, 14, 15, 13, 24, 26, 27, 25,
+ 13, 11, 10, 12, 25, 23, 22, 24,
+ 15, 13, 12, 14, 27, 25, 24, 26,
};
const int pylith::faults::CohesiveDataHex8h::_materialIds[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8i.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8i.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDataHex8i.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -90,8 +90,8 @@
41, 35, 30, 38, 42, 36, 31, 39,
27, 28, 38, 37, 32, 33, 41, 40,
15, 14, 13, 24, 19, 18, 17, 26,
- 24, 23, 20, 21, 41, 40, 37, 38,
- 21, 22, 25, 24, 38, 39, 42, 41,
+ 20, 21, 24, 23, 37, 38, 41, 40,
+ 24, 21, 22, 25, 41, 38, 39, 42,
};
const int pylith::faults::CohesiveDataHex8i::_materialIds[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -60,9 +60,9 @@
3.33333333e-01,};
const double pylith::faults::CohesiveDynDataTet4::_basisDeriv[] = {
- -1.00000000e+00, -1.00000000e+00,
- 1.00000000e+00, 0.00000000e+00,
- 0.00000000e+00, 1.00000000e+00,
+ -0.50000000e+00, -0.50000000e+00,
+ 0.50000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 0.50000000e+00,
};
const double pylith::faults::CohesiveDynDataTet4::_verticesRef[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -59,9 +59,9 @@
3.33333333e-01,};
const double pylith::faults::CohesiveKinDataTet4::_basisDeriv[] = {
- -1.00000000e+00, -1.00000000e+00,
- 1.00000000e+00, 0.00000000e+00,
- 0.00000000e+00, 1.00000000e+00,
+ -0.50000000e+00, -0.50000000e+00,
+ 0.50000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 0.50000000e+00,
};
const double pylith::faults::CohesiveKinDataTet4::_verticesRef[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4e.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4e.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4e.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -59,9 +59,9 @@
3.33333333e-01,};
const double pylith::faults::CohesiveKinDataTet4e::_basisDeriv[] = {
- -1.00000000e+00, -1.00000000e+00,
- 1.00000000e+00, 0.00000000e+00,
- 0.00000000e+00, 1.00000000e+00,
+ -0.50000000e+00, -0.50000000e+00,
+ 0.50000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 0.50000000e+00,
};
const double pylith::faults::CohesiveKinDataTet4e::_verticesRef[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -59,9 +59,9 @@
3.33333333e-01,};
const double pylith::faults::CohesiveKinDataTet4f::_basisDeriv[] = {
- -1.00000000e+00, -1.00000000e+00,
- 1.00000000e+00, 0.00000000e+00,
- 0.00000000e+00, 1.00000000e+00,
+ -0.50000000e+00, -0.50000000e+00,
+ 0.50000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 0.50000000e+00,
};
const double pylith::faults::CohesiveKinDataTet4f::_verticesRef[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinSrcsDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinSrcsDataTet4.cc 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinSrcsDataTet4.cc 2010-09-07 22:59:05 UTC (rev 17178)
@@ -59,9 +59,9 @@
3.33333333e-01,};
const double pylith::faults::CohesiveKinSrcsDataTet4::_basisDeriv[] = {
- -1.00000000e+00, -1.00000000e+00,
- 1.00000000e+00, 0.00000000e+00,
- 0.00000000e+00, 1.00000000e+00,
+ -0.50000000e+00, -0.50000000e+00,
+ 0.50000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 0.50000000e+00,
};
const double pylith::faults::CohesiveKinSrcsDataTet4::_verticesRef[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_cell_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_cell_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_cell_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -12,8 +12,8 @@
0.000000e+00 -1.000000e+00 1.000000e+00
0.000000e+00 1.000000e+00 1.000000e+00
CELLS 2 10
-4 6 7 1 0
-4 4 5 3 2
+4 0 6 7 1
+4 2 4 5 3
CELL_TYPES 2
9
9
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -12,8 +12,8 @@
0.000000e+00 -1.000000e+00 1.000000e+00
0.000000e+00 1.000000e+00 1.000000e+00
CELLS 2 10
-4 6 7 1 0
-4 4 5 3 2
+4 0 6 7 1
+4 2 4 5 3
CELL_TYPES 2
9
9
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_vertex_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_vertex_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_bc_vertex_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -12,8 +12,8 @@
0.000000e+00 -1.000000e+00 1.000000e+00
0.000000e+00 1.000000e+00 1.000000e+00
CELLS 2 10
-4 6 7 1 0
-4 4 5 3 2
+4 0 6 7 1
+4 2 4 5 3
CELL_TYPES 2
9
9
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_cell_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_cell_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_cell_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -8,7 +8,7 @@
0.000000e+00 -1.000000e+00 1.000000e+00
0.000000e+00 1.000000e+00 1.000000e+00
CELLS 1 5
-4 2 3 1 0
+4 1 0 2 3
CELL_TYPES 1
9
CELL_DATA 1
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -8,6 +8,6 @@
0.000000e+00 -1.000000e+00 1.000000e+00
0.000000e+00 1.000000e+00 1.000000e+00
CELLS 1 5
-4 2 3 1 0
+4 1 0 2 3
CELL_TYPES 1
9
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_vertex_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_vertex_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_fault_vertex_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -8,7 +8,7 @@
0.000000e+00 -1.000000e+00 1.000000e+00
0.000000e+00 1.000000e+00 1.000000e+00
CELLS 1 5
-4 2 3 1 0
+4 1 0 2 3
CELL_TYPES 1
9
POINT_DATA 4
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_cell_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_cell_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_cell_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -10,8 +10,8 @@
1.000000e+00 -1.000000e+00 1.000000e+00
1.000000e+00 1.000000e+00 1.000000e+00
CELLS 2 10
-4 2 3 1 0
-4 4 5 3 2
+4 0 2 3 1
+4 2 4 5 3
CELL_TYPES 2
9
9
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -10,8 +10,8 @@
1.000000e+00 -1.000000e+00 1.000000e+00
1.000000e+00 1.000000e+00 1.000000e+00
CELLS 2 10
-4 2 3 1 0
-4 4 5 3 2
+4 0 2 3 1
+4 2 4 5 3
CELL_TYPES 2
9
9
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_vertex_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_vertex_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/hex8_surf_vertex_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -10,8 +10,8 @@
1.000000e+00 -1.000000e+00 1.000000e+00
1.000000e+00 1.000000e+00 1.000000e+00
CELLS 2 10
-4 2 3 1 0
-4 4 5 3 2
+4 0 2 3 1
+4 2 4 5 3
CELL_TYPES 2
9
9
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_cell_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_cell_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_cell_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -7,7 +7,7 @@
CELLS 1 2
1 0
CELL_TYPES 1
--1
+1
CELL_DATA 1
VECTORS traction double
1.100000e+00 2.200000e+00 0.0
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -7,4 +7,4 @@
CELLS 1 2
1 0
CELL_TYPES 1
--1
+1
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_vertex_t10.vtk
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_vertex_t10.vtk 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/line2_surf_vertex_t10.vtk 2010-09-07 22:59:05 UTC (rev 17178)
@@ -7,7 +7,7 @@
CELLS 1 2
1 0
CELL_TYPES 1
--1
+1
POINT_DATA 1
VECTORS displacements double
1.100000e+00 2.200000e+00 0.0
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldsNewMesh.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldsNewMesh.hh 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldsNewMesh.hh 2010-09-07 22:59:05 UTC (rev 17178)
@@ -53,9 +53,9 @@
CPPUNIT_TEST( testAllocateSequence );
CPPUNIT_TEST( testAllocateArray );
CPPUNIT_TEST( testAllocateDomain );
-#endif
CPPUNIT_TEST( testGet );
CPPUNIT_TEST( testGetConst );
+#endif
CPPUNIT_TEST( testMesh );
CPPUNIT_TEST( testFieldNames );
Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py 2010-09-07 22:59:05 UTC (rev 17178)
@@ -72,6 +72,7 @@
def N1p(self, p):
return 0.5
+
# ----------------------------------------------------------------------
class Line3(object):
@@ -124,6 +125,7 @@
def N2p(self, p):
return -2.0*p
+
# ----------------------------------------------------------------------
class Quad4(object):
@@ -137,8 +139,8 @@
[-1.0, +1.0]])
quadPts = numpy.array([ [-1.0/3**0.5, -1.0/3**0.5],
[+1.0/3**0.5, -1.0/3**0.5],
- [+1.0/3**0.5, +1.0/3**0.5],
- [-1.0/3**0.5, +1.0/3**0.5] ])
+ [-1.0/3**0.5, +1.0/3**0.5],
+ [+1.0/3**0.5, +1.0/3**0.5] ])
quadWts = numpy.array( [1.0, 1.0, 1.0, 1.0])
# Compute basis fns and derivatives at quadrature points
@@ -203,7 +205,152 @@
def N3q(self, p):
return (1-p[0])/4.0
+
# ----------------------------------------------------------------------
+class Quad9(object):
+
+ def __init__(self):
+ """
+ Setup quad9 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0], # corners
+ [+1.0, -1.0],
+ [+1.0, +1.0],
+ [-1.0, +1.0],
+ [ 0.0, -1.0], # edges
+ [+1.0, 0.0],
+ [ 0.0, +1.0],
+ [-1.0, 0.0],
+ [ 0.0, 0.0], # face
+ ])
+ quadPts = numpy.array([ [-(3.0/5)**0.5, -(3.0/5)**0.5],
+ [ 0.0, -(3.0/5)**0.5],
+ [+(3.0/5)**0.5, -(3.0/5)**0.5],
+ [-(3.0/5)**0.5, 0.0],
+ [ 0.0, 0.0],
+ [+(3.0/5)**0.5, 0.0],
+ [-(3.0/5)**0.5, +(3.0/5)**0.5],
+ [ 0.0, +(3.0/5)**0.5],
+ [+(3.0/5)**0.5, +(3.0/5)**0.5] ])
+ quadWts = numpy.array( [25.0/81, 40.0/81, 25.0/81,
+ 40.0/81, 64.0/81, 40.0/81,
+ 25.0/81, 40.0/81, 25.0/81])
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (9, 9), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (9, 9, 2), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q), self.N2(q),
+ self.N3(q), self.N4(q), self.N5(q),
+ self.N6(q), self.N7(q), self.N8(q)],
+ dtype=numpy.float64).reshape( (9,) )
+ deriv = numpy.array([[self.N0p(q), self.N0q(q)],
+ [self.N1p(q), self.N1q(q)],
+ [self.N2p(q), self.N2q(q)],
+ [self.N3p(q), self.N3q(q)],
+ [self.N4p(q), self.N4q(q)],
+ [self.N5p(q), self.N5q(q)],
+ [self.N6p(q), self.N6q(q)],
+ [self.N7p(q), self.N7q(q)],
+ [self.N8p(q), self.N8q(q)]])
+ basisDeriv[iQuad] = deriv.reshape((9, 2))
+ iQuad += 1
+
+ self.cellDim = 2
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+ def N0(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])
+
+ def N0p(self, p):
+ return (p[0]-0.5)*(-0.5*p[1])*(1.0-p[1])
+
+ def N0q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]-0.5)
+
+ def N1(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])
+
+ def N1p(self, p):
+ return (p[0]+0.5)*(-0.5*p[1])*(1.0-p[1])
+
+ def N1q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]-0.5)
+
+ def N2(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])
+
+ def N2p(self, p):
+ return (p[0]+0.5)*(+0.5*p[1])*(1.0+p[1])
+
+ def N2q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]+0.5)
+
+ def N3(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])
+
+ def N3p(self, p):
+ return (p[0]-0.5)*(+0.5*p[1])*(1.0+p[1])
+
+ def N3q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]+0.5)
+
+ def N4(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])
+
+ def N4p(self, p):
+ return (-2.0*p[0])*(-0.5*p[1])*(1.0-p[1])
+
+ def N4q(self, p):
+ return (1.0-p[0]**2)*(p[1]-0.5)
+
+ def N5(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)
+
+ def N5p(self, p):
+ return (p[0]+0.5)*(1.0-p[1]**2)
+
+ def N5q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-2.0*p[1])
+
+ def N6(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])
+
+ def N6p(self, p):
+ return (-2.0*p[0])*(+0.5*p[1])*(1.0+p[1])
+
+ def N6q(self, p):
+ return (1.0-p[0]**2)*(p[1]+0.5)
+
+ def N7(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)
+
+ def N7p(self, p):
+ return (p[0]-0.5)*(1.0-p[1]**2)
+
+ def N7q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-2.0*p[1])
+
+ def N8(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)
+
+ def N8p(self, p):
+ return (-2.0*p[0])*(1.0-p[1]**2)
+
+ def N8q(self, p):
+ return (1.0-p[0]**2)*(-2.0*p[1])
+
+
+# ----------------------------------------------------------------------
class Hex8(object):
def __init__(self):
@@ -220,12 +367,12 @@
[-1.0, +1.0, +1.0]])
quadPts = numpy.array([ [-1.0/3**0.5, -1.0/3**0.5, -1.0/3**0.5],
[+1.0/3**0.5, -1.0/3**0.5, -1.0/3**0.5],
+ [-1.0/3**0.5, +1.0/3**0.5, -1.0/3**0.5],
[+1.0/3**0.5, +1.0/3**0.5, -1.0/3**0.5],
- [-1.0/3**0.5, +1.0/3**0.5, -1.0/3**0.5],
[-1.0/3**0.5, -1.0/3**0.5, +1.0/3**0.5],
[+1.0/3**0.5, -1.0/3**0.5, +1.0/3**0.5],
- [+1.0/3**0.5, +1.0/3**0.5, +1.0/3**0.5],
- [-1.0/3**0.5, +1.0/3**0.5, +1.0/3**0.5]])
+ [-1.0/3**0.5, +1.0/3**0.5, +1.0/3**0.5],
+ [+1.0/3**0.5, +1.0/3**0.5, +1.0/3**0.5]])
quadWts = numpy.array( [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
# Compute basis fns and derivatives at quadrature points
@@ -356,7 +503,494 @@
def N7r(self, p):
return (1-p[0])*(1+p[1])/8.0
+
# ----------------------------------------------------------------------
+class Hex27(object):
+
+ def __init__(self):
+ """
+ Setup hex8 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0, -1.0], # Corners
+ [+1.0, -1.0, -1.0],
+ [+1.0, +1.0, -1.0],
+ [-1.0, +1.0, -1.0],
+ [-1.0, -1.0, +1.0],
+ [+1.0, -1.0, +1.0],
+ [+1.0, +1.0, +1.0],
+ [-1.0, +1.0, +1.0],
+ [ 0.0, -1.0, -1.0], # Bottom edges
+ [+1.0, 0.0, -1.0],
+ [ 0.0, +1.0, -1.0],
+ [-1.0, 0.0, -1.0],
+ [ 0.0, -1.0, +1.0], # Top edges
+ [+1.0, 0.0, +1.0],
+ [ 0.0, +1.0, +1.0],
+ [-1.0, 0.0, +1.0],
+ [-1.0, -1.0, 0.0], # Middle edges
+ [+1.0, -1.0, 0.0],
+ [+1.0, +1.0, 0.0],
+ [-1.0, +1.0, 0.0],
+ [-1.0, 0.0, 0.0], # Faces
+ [+1.0, 0.0, 0.0],
+ [ 0.0, -1.0, 0.0],
+ [ 0.0, +1.0, 0.0],
+ [ 0.0, 0.0, -1.0],
+ [ 0.0, 0.0, +1.0],
+ [ 0.0, 0.0, 0.0]]) # Interior
+ quadPts = numpy.array([ [-(3.0/5)**0.5, -(3.0/5)**0.5, -(3.0/5)**0.5],
+ [ 0.0, -(3.0/5)**0.5, -(3.0/5)**0.5],
+ [+(3.0/5)**0.5, -(3.0/5)**0.5, -(3.0/5)**0.5],
+ [-(3.0/5)**0.5, 0.0, -(3.0/5)**0.5],
+ [ 0.0, 0.0, -(3.0/5)**0.5],
+ [+(3.0/5)**0.5, 0.0, -(3.0/5)**0.5],
+ [-(3.0/5)**0.5, +(3.0/5)**0.5, -(3.0/5)**0.5],
+ [ 0.0, +(3.0/5)**0.5, -(3.0/5)**0.5],
+ [+(3.0/5)**0.5, +(3.0/5)**0.5, -(3.0/5)**0.5],
+ [-(3.0/5)**0.5, -(3.0/5)**0.5, 0.0],
+ [ 0.0, -(3.0/5)**0.5, 0.0],
+ [+(3.0/5)**0.5, -(3.0/5)**0.5, 0.0],
+ [-(3.0/5)**0.5, 0.0, 0.0],
+ [ 0.0, 0.0, 0.0],
+ [+(3.0/5)**0.5, 0.0, 0.0],
+ [-(3.0/5)**0.5, +(3.0/5)**0.5, 0.0],
+ [ 0.0, +(3.0/5)**0.5, 0.0],
+ [+(3.0/5)**0.5, +(3.0/5)**0.5, 0.0],
+ [-(3.0/5)**0.5, -(3.0/5)**0.5, +(3.0/5)**0.5],
+ [ 0.0, -(3.0/5)**0.5, +(3.0/5)**0.5],
+ [+(3.0/5)**0.5, -(3.0/5)**0.5, +(3.0/5)**0.5],
+ [-(3.0/5)**0.5, 0.0, +(3.0/5)**0.5],
+ [ 0.0, 0.0, +(3.0/5)**0.5],
+ [+(3.0/5)**0.5, 0.0, +(3.0/5)**0.5],
+ [-(3.0/5)**0.5, +(3.0/5)**0.5, +(3.0/5)**0.5],
+ [ 0.0, +(3.0/5)**0.5, +(3.0/5)**0.5],
+ [+(3.0/5)**0.5, +(3.0/5)**0.5, +(3.0/5)**0.5] ])
+ quadWts = numpy.array( [125.0/729, 200.0/729, 125.0/729,
+ 200.0/729, 320.0/729, 200.0/729,
+ 125.0/729, 200.0/729, 125.0/729,
+ 200.0/729, 320.0/729, 200.0/729,
+ 320.0/729, 512.0/729, 320.0/729,
+ 200.0/729, 320.0/729, 200.0/729,
+ 125.0/729, 200.0/729, 125.0/729,
+ 200.0/729, 320.0/729, 200.0/729,
+ 125.0/729, 200.0/729, 125.0/729])
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (27, 27), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (27, 27, 3), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q), self.N2(q), self.N3(q), # Corners
+ self.N4(q), self.N5(q), self.N6(q), self.N7(q),
+ self.N8(q), self.N9(q), self.N10(q), self.N11(q), # Edges
+ self.N12(q), self.N13(q), self.N14(q), self.N15(q),
+ self.N16(q), self.N17(q), self.N18(q), self.N19(q),
+ self.N20(q), # Interior
+ self.N21(q), self.N22(q), # Faces
+ self.N23(q), self.N24(q),
+ self.N25(q), self.N26(q)],
+ dtype=numpy.float64).reshape( (27,) )
+ deriv = numpy.array([[self.N0p(q), self.N0q(q), self.N0r(q)],
+ [self.N1p(q), self.N1q(q), self.N1r(q)],
+ [self.N2p(q), self.N2q(q), self.N2r(q)],
+ [self.N3p(q), self.N3q(q), self.N3r(q)],
+ [self.N4p(q), self.N4q(q), self.N4r(q)],
+ [self.N5p(q), self.N5q(q), self.N5r(q)],
+ [self.N6p(q), self.N6q(q), self.N6r(q)],
+ [self.N7p(q), self.N7q(q), self.N7r(q)],
+ [self.N8p(q), self.N8q(q), self.N8r(q)],
+ [self.N9p(q), self.N9q(q), self.N9r(q)],
+ [self.N10p(q), self.N10q(q), self.N10r(q)],
+ [self.N11p(q), self.N11q(q), self.N11r(q)],
+ [self.N12p(q), self.N12q(q), self.N12r(q)],
+ [self.N13p(q), self.N13q(q), self.N13r(q)],
+ [self.N14p(q), self.N14q(q), self.N14r(q)],
+ [self.N15p(q), self.N15q(q), self.N15r(q)],
+ [self.N16p(q), self.N16q(q), self.N16r(q)],
+ [self.N17p(q), self.N17q(q), self.N17r(q)],
+ [self.N18p(q), self.N18q(q), self.N18r(q)],
+ [self.N19p(q), self.N19q(q), self.N19r(q)],
+ [self.N20p(q), self.N20q(q), self.N20r(q)],
+ [self.N21p(q), self.N21q(q), self.N21r(q)],
+ [self.N22p(q), self.N22q(q), self.N22r(q)],
+ [self.N23p(q), self.N23q(q), self.N23r(q)],
+ [self.N24p(q), self.N24q(q), self.N24r(q)],
+ [self.N25p(q), self.N25q(q), self.N25r(q)],
+ [self.N26p(q), self.N26q(q), self.N26r(q)]])
+ basisDeriv[iQuad] = deriv.reshape((27, 3))
+ iQuad += 1
+
+ self.cellDim = 3
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+ # Corners
+ def N0(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N0p(self, p):
+ return (p[0]-0.5)*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N0q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]-0.5)*(-0.5*p[2])*(1.0-p[2])
+
+ def N0r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(p[2]-0.5)
+
+
+ def N1(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N1p(self, p):
+ return (p[0]+0.5)*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N1q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]-0.5)*(-0.5*p[2])*(1.0-p[2])
+
+ def N1r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(p[2]-0.5)
+
+
+ def N2(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N2p(self, p):
+ return (p[0]+0.5)*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N2q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]+0.5)*(-0.5*p[2])*(1.0-p[2])
+
+ def N2r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(p[2]-0.5)
+
+
+ def N3(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N3p(self, p):
+ return (p[0]-0.5)*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N3q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]+0.5)*(-0.5*p[2])*(1.0-p[2])
+
+ def N3r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(p[2]-0.5)
+
+
+ def N4(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N4p(self, p):
+ return (p[0]-0.5)*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N4q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]-0.5)*(+0.5*p[2])*(1.0+p[2])
+
+ def N4r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(p[2]+0.5)
+
+
+ def N5(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N5p(self, p):
+ return (p[0]+0.5)*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N5q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]-0.5)*(+0.5*p[2])*(1.0+p[2])
+
+ def N5r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(p[2]+0.5)
+
+
+ def N6(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N6p(self, p):
+ return (p[0]+0.5)*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N6q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]+0.5)*(+0.5*p[2])*(1.0+p[2])
+
+ def N6r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(p[2]+0.5)
+
+
+ def N7(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N7p(self, p):
+ return (p[0]-0.5)*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N7q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]+0.5)*(+0.5*p[2])*(1.0+p[2])
+
+ def N7r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(p[2]+0.5)
+
+
+ # Bottom edges
+ def N8(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N8p(self, p):
+ return (-2.0*p[0])*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N8q(self, p):
+ return (1.0-p[0]**2)*(p[1]-0.5)*(-0.5*p[2])*(1.0-p[2])
+
+ def N8r(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(p[2]-0.5)
+
+
+ def N9(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N9p(self, p):
+ return (p[0]+0.5)*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N9q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-2.0*p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N9r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(p[2]-0.5)
+
+
+ def N10(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N10p(self, p):
+ return (-2.0*p[0])*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N10q(self, p):
+ return (1.0-p[0]**2)*(p[1]+0.5)*(-0.5*p[2])*(1.0-p[2])
+
+ def N10r(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(p[2]-0.5)
+
+
+ def N11(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N11p(self, p):
+ return (p[0]-0.5)*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N11q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-2.0*p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N11r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(p[2]-0.5)
+
+
+ # Top edges
+ def N12(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N12p(self, p):
+ return (-2.0*p[0])*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N12q(self, p):
+ return (1.0-p[0]**2)*(p[1]-0.5)*(+0.5*p[2])*(1.0+p[2])
+
+ def N12r(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(p[2]+0.5)
+
+
+ def N13(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N13p(self, p):
+ return (p[0]+0.5)*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N13q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-2.0*p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N13r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(p[2]+0.5)
+
+
+ def N14(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N14p(self, p):
+ return (-2.0*p[0])*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N14q(self, p):
+ return (1.0-p[0]**2)*(p[1]+0.5)*(+0.5*p[2])*(1.0+p[2])
+
+ def N14r(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(p[2]+0.5)
+
+
+ def N15(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N15p(self, p):
+ return (p[0]-0.5)*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N15q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-2.0*p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N15r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(p[2]+0.5)
+
+
+ # Middle edges
+ def N16(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N16p(self, p):
+ return (p[0]-0.5)*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N16q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]-0.5)*(1.0-p[2]**2)
+
+ def N16r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(-2.0*p[2])
+
+
+ def N17(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N17p(self, p):
+ return (p[0]+0.5)*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N17q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]-0.5)*(1.0-p[2]**2)
+
+ def N17r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(-2.0*p[2])
+
+
+ def N18(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N18p(self, p):
+ return (p[0]+0.5)*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N18q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]+0.5)*(1.0-p[2]**2)
+
+ def N18r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(-2.0*p[2])
+
+
+ def N19(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N19p(self, p):
+ return (p[0]-0.5)*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N19q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]+0.5)*(1.0-p[2]**2)
+
+ def N19r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(-2.0*p[2])
+
+
+ # Left/right
+ def N20(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N20p(self, p):
+ return (p[0]-0.5)*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N20q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-2.0*p[1])*(1.0-p[2]**2)
+
+ def N20r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(-2.0*p[2])
+
+
+ def N21(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N21p(self, p):
+ return (p[0]+0.5)*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N21q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-2.0*p[1])*(1.0-p[2]**2)
+
+ def N21r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(-2.0*p[2])
+
+
+ # Front/back
+ def N22(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N22p(self, p):
+ return (-2.0*p[0])*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N22q(self, p):
+ return (1.0-p[0]**2)*(p[1]-0.5)*(1.0-p[2]**2)
+
+ def N22r(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(-2.0*p[2])
+
+
+ def N23(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N23p(self, p):
+ return (-2.0*p[0])*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N23q(self, p):
+ return (1.0-p[0]**2)*(p[1]+0.5)*(1.0-p[2]**2)
+
+ def N23r(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(-2.0*p[2])
+
+
+ # Bottom/top
+ def N24(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N24p(self, p):
+ return (-2.0*p[0])*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N24q(self, p):
+ return (1.0-p[0]**2)*(-2.0*p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N24r(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(p[2]-0.5)
+
+
+ def N25(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N25p(self, p):
+ return (-2.0*p[0])*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N25q(self, p):
+ return (1.0-p[0]**2)*(-2.0*p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N25r(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(p[2]+0.5)
+
+
+ # Interior
+ def N26(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N26p(self, p):
+ return (-2.0*p[0])*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N26q(self, p):
+ return (1.0-p[0]**2)*(-2.0*p[1])*(1.0-p[2]**2)
+
+ def N26r(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(-2.0*p[2])
+
+
+# ----------------------------------------------------------------------
class TestFIATLagrange(unittest.TestCase):
"""
Unit testing of FIATLagrange object.
@@ -366,11 +1000,11 @@
"""
Test initialize() with line2 cell.
"""
- from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATLagrange()
- cell.cellDim = 1
- cell.degree = 1
- cell.order = 1
+ cell.inventory.dimension = 1
+ cell.inventory.degree = 1
+ cell.inventory.order = 1
+ cell._configure()
cell.initialize(spaceDim=1)
cellE = Line2()
@@ -384,11 +1018,11 @@
"""
Test initialize() with line3 cell.
"""
- from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATLagrange()
- cell.cellDim = 1
- cell.degree = 2
- cell.order = 2
+ cell.inventory.dimension = 1
+ cell.inventory.degree = 2
+ cell.inventory.order = 2
+ cell._configure()
cell.initialize(spaceDim=2)
cellE = Line3()
@@ -402,11 +1036,11 @@
"""
Test initialize() with quad4 cell.
"""
- from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATLagrange()
- cell.cellDim = 2
- cell.degree = 1
- cell.order = 2
+ cell.inventory.dimension = 2
+ cell.inventory.degree = 1
+ cell.inventory.order = 2
+ cell._configure()
cell.initialize(spaceDim=2)
cellE = Quad4()
@@ -416,15 +1050,33 @@
return
+ def test_initialize_quad9(self):
+ """
+ Test initialize() with quad9 cell.
+ """
+ cell = FIATLagrange()
+ cell.inventory.dimension = 2
+ cell.inventory.degree = 2
+ cell.inventory.order = 5
+ cell._configure()
+ cell.initialize(spaceDim=2)
+
+ cellE = Quad9()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryQuad2D
+ self.failUnless(isinstance(cell.geometry, GeometryQuad2D))
+ return
+
+
def test_initialize_hex8(self):
"""
Test initialize() with hex8 cell.
"""
- from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATLagrange()
- cell.cellDim = 3
- cell.degree = 1
- cell.order = 2
+ cell.inventory.dimension = 3
+ cell.inventory.degree = 1
+ cell.inventory.order = 2
+ cell._configure()
cell.initialize(spaceDim=3)
cellE = Hex8()
@@ -434,6 +1086,24 @@
return
+ def test_initialize_hex27(self):
+ """
+ Test initialize() with hex27 cell.
+ """
+ cell = FIATLagrange()
+ cell.inventory.dimension = 3
+ cell.inventory.degree = 2
+ cell.inventory.order = 5
+ cell._configure()
+ cell.initialize(spaceDim=3)
+
+ cellE = Hex27()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryHex3D
+ self.failUnless(isinstance(cell.geometry, GeometryHex3D))
+ return
+
+
def test_factory(self):
"""
Test factory method.
@@ -445,7 +1115,7 @@
def _checkVals(self, cellE, cell):
"""
- Check known values against those generated by FIATSimplex.
+ Check known values against those generated by FIATLagrange.
"""
# Check basic attributes
Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATSimplex.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATSimplex.py 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATSimplex.py 2010-09-07 22:59:05 UTC (rev 17178)
@@ -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,7 +195,116 @@
def N2q(self, p):
return 0.5
+
# ----------------------------------------------------------------------
+class Tri6(object):
+
+ def __init__(self):
+ """
+ Setup tri33 cell.
+ """
+ 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]])
+ quadPts = numpy.array([ [-0.64288254, -0.68989795],
+ [-0.84993778, 0.28989795],
+ [ 0.33278049, -0.68989795],
+ [-0.43996017, 0.28989795]])
+ quadWts = numpy.array( [0.63608276, 0.36391724, 0.63608276, 0.36391724])
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (4, 6), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (4, 6, 2), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q), self.N2(q),
+ self.N3(q), self.N4(q), self.N5(q)],
+ dtype=numpy.float64).reshape( (6,) )
+ deriv = numpy.array([[self.N0p(q), self.N0q(q)],
+ [self.N1p(q), self.N1q(q)],
+ [self.N2p(q), self.N2q(q)],
+ [self.N3p(q), self.N3q(q)],
+ [self.N4p(q), self.N4q(q)],
+ [self.N5p(q), self.N5q(q)]])
+ print deriv
+ basisDeriv[iQuad] = deriv.reshape((6, 2))
+ iQuad += 1
+
+ self.cellDim = 2
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+ def N0(self, p):
+ return (-p[0]-p[1])*(1+p[0])
+
+ def N0p(self, p):
+ return -1.0-2*p[0]-p[1]
+
+ def N0q(self, p):
+ return -(1+p[0])
+
+
+ def N1(self, p):
+ return (1.0+p[0])*(1+p[1])
+
+ def N1p(self, p):
+ return (1+p[1])
+
+ def N1q(self, p):
+ return (1.0+p[0])
+
+
+ def N2(self, p):
+ return (-p[0]-p[1])*(1+p[1])
+
+ def N2p(self, p):
+ return -(1+p[1])
+
+ def N2q(self, p):
+ return -1.0-p[0]-2*p[1]
+
+
+ def N3(self, p):
+ return 0.5*(-p[0]-p[1])*(-1.0-p[0]-p[1])
+
+ def N3p(self, p):
+ return 0.5+p[0]+p[1]
+
+ def N3q(self, p):
+ return 0.5+p[0]+p[1]
+
+
+ def N4(self, p):
+ return 0.5*(1.0+p[0])*(p[0])
+
+ def N4p(self, p):
+ return 0.5+p[0]
+
+ def N4q(self, p):
+ return 0
+
+
+ def N5(self, p):
+ return 0.5*(1.0+p[1])*(p[1])
+
+ def N5p(self, p):
+ return 0
+
+ def N5q(self, p):
+ return 0.5+p[1]
+
+
+# ----------------------------------------------------------------------
class Tet4(object):
def __init__(self):
@@ -276,6 +392,7 @@
def N3r(self, p):
return 0.5
+
# ----------------------------------------------------------------------
class TestFIATSimplex(unittest.TestCase):
"""
@@ -344,9 +461,9 @@
Test initialize() with tri3 cell.
"""
cell = FIATSimplex()
- cell.shape = "triangle"
- cell.degree = 1
- cell.order = 1
+ cell.inventory.shape = "triangle"
+ cell.inventory.degree = 1
+ cell._configure()
cell.initialize(spaceDim=2)
cellE = Tri3()
@@ -356,14 +473,31 @@
return
+ def test_initialize_tri6(self):
+ """
+ Test initialize() with tri6 cell.
+ """
+ cell = FIATSimplex()
+ cell.inventory.shape = "triangle"
+ cell.inventory.degree = 2
+ cell._configure()
+ cell.initialize(spaceDim=2)
+
+ cellE = Tri6()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryTri2D
+ self.failUnless(isinstance(cell.geometry, GeometryTri2D))
+ return
+
+
def test_initialize_tet4(self):
"""
Test initialize() with tet4 cell.
"""
cell = FIATSimplex()
- cell.shape = "tetrahedron"
- cell.degree = 1
- cell.order = 1
+ cell.inventory.shape = "tetrahedron"
+ cell.inventory.degree = 1
+ cell._configure()
cell.initialize(spaceDim=3)
cellE = Tet4()
Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestMeshQuadrature.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestMeshQuadrature.py 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestMeshQuadrature.py 2010-09-07 22:59:05 UTC (rev 17178)
@@ -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/trunk/unittests/pytests/feassemble/TestSubMeshQuadrature.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestSubMeshQuadrature.py 2010-09-07 22:35:37 UTC (rev 17177)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestSubMeshQuadrature.py 2010-09-07 22:59:05 UTC (rev 17178)
@@ -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