[cig-commits] r15011 - in short/3D/PyLith/trunk: . libsrc/bc libsrc/faults libsrc/feassemble libsrc/materials libsrc/meshio libsrc/problems libsrc/topology modulesrc/meshio modulesrc/topology pylith pylith/materials pylith/meshio pylith/problems pylith/utils tests/1d/line2 unittests/libtests/bc unittests/libtests/bc/data unittests/libtests/faults unittests/libtests/faults/data unittests/libtests/feassemble unittests/libtests/feassemble/data unittests/libtests/materials unittests/libtests/materials/data unittests/libtests/meshio unittests/libtests/topology unittests/pytests/bc unittests/pytests/faults unittests/pytests/feassemble/data unittests/pytests/meshio/data unittests/pytests/topology
brad at geodynamics.org
brad at geodynamics.org
Mon May 18 16:02:10 PDT 2009
Author: brad
Date: 2009-05-18 16:02:05 -0700 (Mon, 18 May 2009)
New Revision: 15011
Added:
short/3D/PyLith/trunk/pylith/utils/NullComponent.py
Modified:
short/3D/PyLith/trunk/README
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/materials/Material.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.hh
short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh
short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshIOCubit.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.cc
short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc
short/3D/PyLith/trunk/libsrc/topology/Mesh.cc
short/3D/PyLith/trunk/libsrc/topology/Mesh.hh
short/3D/PyLith/trunk/libsrc/topology/SolutionFields.cc
short/3D/PyLith/trunk/libsrc/topology/SolutionFields.hh
short/3D/PyLith/trunk/modulesrc/meshio/MeshIOObj.i
short/3D/PyLith/trunk/modulesrc/topology/Mesh.i
short/3D/PyLith/trunk/modulesrc/topology/SolutionFields.i
short/3D/PyLith/trunk/pylith/Makefile.am
short/3D/PyLith/trunk/pylith/materials/ElasticMaterial.py
short/3D/PyLith/trunk/pylith/materials/Material.py
short/3D/PyLith/trunk/pylith/meshio/CellFilter.py
short/3D/PyLith/trunk/pylith/meshio/CellFilterAvgMesh.py
short/3D/PyLith/trunk/pylith/meshio/CellFilterAvgSubMesh.py
short/3D/PyLith/trunk/pylith/meshio/MeshIOObj.py
short/3D/PyLith/trunk/pylith/meshio/OutputManager.py
short/3D/PyLith/trunk/pylith/meshio/OutputManagerMesh.py
short/3D/PyLith/trunk/pylith/meshio/OutputManagerSubMesh.py
short/3D/PyLith/trunk/pylith/meshio/OutputSolnSubset.py
short/3D/PyLith/trunk/pylith/meshio/VertexFilter.py
short/3D/PyLith/trunk/pylith/problems/Explicit.py
short/3D/PyLith/trunk/pylith/problems/Formulation.py
short/3D/PyLith/trunk/pylith/problems/Implicit.py
short/3D/PyLith/trunk/pylith/problems/Problem.py
short/3D/PyLith/trunk/pylith/problems/TimeStepUniform.py
short/3D/PyLith/trunk/pylith/problems/TimeStepUser.py
short/3D/PyLith/trunk/pylith/utils/__init__.py
short/3D/PyLith/trunk/tests/1d/line2/matprops.spatialdb
short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryMesh.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBoundary.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticisotropic3d.spatialdb
short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticplanestrain.spatialdb
short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticstrain1d.spatialdb
short/3D/PyLith/trunk/unittests/libtests/bc/test_bc.cc
short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_1d.spatialdb
short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_2d.spatialdb
short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_3d.spatialdb
short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticisotropic3d.spatialdb
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticplanestrain.spatialdb
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticstrain1d.spatialdb
short/3D/PyLith/trunk/unittests/libtests/feassemble/test_feassemble.cc
short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/matinitialize.spatialdb
short/3D/PyLith/trunk/unittests/libtests/materials/test_materials.cc
short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.hh
short/3D/PyLith/trunk/unittests/libtests/topology/TestSolutionFields.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestSolutionFields.hh
short/3D/PyLith/trunk/unittests/pytests/bc/TestAbsorbingDampers.py
short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py
short/3D/PyLith/trunk/unittests/pytests/feassemble/data/elasticplanestrain.spatialdb
short/3D/PyLith/trunk/unittests/pytests/meshio/data/mesh2Din3D.txt
short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py
short/3D/PyLith/trunk/unittests/pytests/topology/TestSolutionFields.py
Log:
Cleaned up nondimensionalization of coordinates (create coordinates_dimensioned section). Cleaned up SolutionFields (remove time history stuff). Updated tests for use of units in spatial databases (Python syntax for units strings). Setup integrators to not precompute geometry (use define to allow switching back).
Modified: short/3D/PyLith/trunk/README
===================================================================
--- short/3D/PyLith/trunk/README 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/README 2009-05-18 23:02:05 UTC (rev 15011)
@@ -59,9 +59,14 @@
(5) The DirichletPoints boundary condition has been renamed to
DirichletBC.
+(6) Procedure for enabling certain features no longer involves setting
+a "use" property to True. Instead, the features is enable when the
+user sets the component to a facility. This applies to gravity,
+initial stresses, initial strains, and initial state variables.
+
STILL ON THE TODO LIST
-(5) Nondimensionalization of the problem eliminates the need to
+(7) Nondimensionalization of the problem eliminates the need to
condition the fault constraints. The "mat_db" facility was removed.
======================================================================
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/TODO 2009-05-18 23:02:05 UTC (rev 15011)
@@ -40,21 +40,15 @@
1. SNES
- Cleanup SolutionFields.
- Remove time history stuff. [TEST]
- Replace solveSolnVec with solution. [TEST]
+ Cleanup SolutionFields. [full scale test]
Power-law rheology
2. Nondimensionalization
- Add units to spatialdata. [UPDATE TESTING]
-
Add Quadrature::computeGeometry(coordinatesCell)
Use #define PRECOMPUTE_GEOMETRY [TEST]
- MeshIO: coordinates and coordinates_dimensioned [UPDATE TESTING]
-
DataWriteVTK needs dimensionalizer
Done coordinates in output (nondimensional instead of dimensioned)
@@ -96,6 +90,8 @@
4. Add missing unit tests
+ libtests/topology/TestMesh::testNondimensionalize()
+
libtests/topology/Field add constraints to field in unit tests
copy
+=
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -29,11 +29,14 @@
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
+//#define PRECOMPUTE_GEOMETRY
+
// ----------------------------------------------------------------------
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
typedef pylith::topology::SubMesh::RealSection SubRealSection;
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
// ----------------------------------------------------------------------
// Default constructor.
@@ -82,21 +85,21 @@
const int numCorners = _quadrature->numBasis();
// Get 'surface' cells (1 dimension lower than top-level cells)
- const ALE::Obj<SieveSubMesh>& submesh = _boundaryMesh->sieveMesh();
- assert(!submesh.isNull());
+ const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
+ assert(!sieveSubMesh.isNull());
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- submesh->heightStratum(1);
+ sieveSubMesh->heightStratum(1);
assert(!cells.isNull());
const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
// Make sure surface cells are compatible with quadrature.
- const int boundaryDepth = submesh->depth()-1; // depth of bndry cells
+ const int boundaryDepth = sieveSubMesh->depth()-1; // depth of bndry cells
for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
++c_iter) {
const int cellNumCorners =
- submesh->getNumCellCorners(*c_iter, boundaryDepth);
+ sieveSubMesh->getNumCellCorners(*c_iter, boundaryDepth);
if (numCorners != cellNumCorners) {
std::ostringstream msg;
msg << "Quadrature is incompatible with cell for absorbing boundary "
@@ -128,7 +131,6 @@
double_array jacobian(jacobianSize);
double jacobianDet = 0;
double_array orientation(orientationSize);
- double_array cellVertices(numCorners*spaceDim);
// open database with material property information
_db->open();
@@ -152,15 +154,18 @@
double_array dampingConstsLocal(fiberDim);
double_array dampingConstsGlobal(fiberDim);
+ double_array coordinatesCell(numCorners*spaceDim);
const ALE::Obj<RealSection>& coordinates =
- submesh->getRealSection("coordinates");
+ sieveSubMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- numCorners*spaceDim);
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
const double densityScale = _normalizer->densityScale();
+ assert(_normalizer->timeScale() > 0);
const double velocityScale =
_normalizer->lengthScale() / _normalizer->timeScale();
@@ -171,19 +176,27 @@
// Compute quadrature information
_quadrature->initializeGeometry();
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->computeGeometry(*_boundaryMesh, cells);
+#endif
for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
c_iter != cellsEnd;
++c_iter) {
+ // Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
+
const double_array& quadPtsNondim = _quadrature->quadPts();
const double_array& quadPtsRef = _quadrature->quadPtsRef();
quadPtsGlobal = quadPtsNondim;
_normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
lengthScale);
- coordsVisitor.clear();
- submesh->restrictClosure(*c_iter, coordsVisitor);
dampingConstsGlobal = 0.0;
for(int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
@@ -217,14 +230,15 @@
dampingConstsLocal[spaceDim-1] = constNormal;
// Compute normal/tangential orientation
- submesh->restrictClosure(coordinates, *c_iter,
- &cellVertices[0], cellVertices.size());
memcpy(&quadPtRef[0], &quadPtsRef[iQuad*cellDim],
cellDim*sizeof(double));
- memcpy(&cellVertices[0], coordsVisitor.getValues(),
- cellVertices.size()*sizeof(double));
- cellGeometry.jacobian(&jacobian, &jacobianDet, cellVertices, quadPtRef);
+#if defined(PRECOMPUTE_GEOMETRY)
+ coordsVisitor.clear();
+ sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
+#endif
+ cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell, quadPtRef);
cellGeometry.orientation(&orientation, jacobian, jacobianDet, up);
+ assert(jacobianDet > 0.0);
orientation /= jacobianDet;
for (int iDim=0; iDim < spaceDim; ++iDim) {
@@ -267,10 +281,10 @@
double_array dampersCell(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveSubMesh>& submesh = _boundaryMesh->sieveMesh();
- assert(!submesh.isNull());
+ const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
+ assert(!sieveSubMesh.isNull());
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- submesh->heightStratum(1);
+ sieveSubMesh->heightStratum(1);
assert(!cells.isNull());
const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
@@ -297,6 +311,14 @@
topology::Mesh::RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
numBasis*spaceDim);
+#if !defined(PRECOMPUTE_GEOMETRY)
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveSubMesh->getRealSection("coordinates");
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
+#endif
+
// Get parameters used in integration.
const double dt = _dt;
assert(dt > 0);
@@ -305,18 +327,24 @@
c_iter != cellsEnd;
++c_iter) {
// Get geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
// Reset element vector to zero
_resetCellVector();
// Restrict input fields to cell
dispTVisitor.clear();
- submesh->restrictClosure(*c_iter, dispTVisitor);
+ sieveSubMesh->restrictClosure(*c_iter, dispTVisitor);
const double* dispTCell = dispTVisitor.getValues();
dispTmdtVisitor.clear();
- submesh->restrictClosure(*c_iter, dispTmdtVisitor);
+ sieveSubMesh->restrictClosure(*c_iter, dispTmdtVisitor);
const double* dispTmdtCell = dispTmdtVisitor.getValues();
dampersSection->restrictPoint(*c_iter, &dampersCell[0], dampersCell.size());
@@ -347,7 +375,7 @@
// Assemble cell contribution into field
residualVisitor.clear();
- submesh->updateAdd(*c_iter, residualVisitor);
+ sieveSubMesh->updateAdd(*c_iter, residualVisitor);
} // for
} // integrateResidual
@@ -372,10 +400,10 @@
const int spaceDim = _quadrature->spaceDim();
// Get cell information
- const ALE::Obj<SieveSubMesh>& submesh = _boundaryMesh->sieveMesh();
- assert(!submesh.isNull());
+ const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
+ assert(!sieveSubMesh.isNull());
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- submesh->heightStratum(1);
+ sieveSubMesh->heightStratum(1);
assert(!cells.isNull());
const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
@@ -397,6 +425,14 @@
(int) pow(sieveMesh->getSieve()->getMaxConeSize(),
sieveMesh->depth())*spaceDim);
+#if !defined(PRECOMPUTE_GEOMETRY)
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveSubMesh->getRealSection("coordinates");
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
+#endif
+
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
assert(0 != jacobianMat);
@@ -411,8 +447,14 @@
for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
++c_iter) {
- // Get geometry information for current cell
+ // Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
// Reset element vector to zero
_resetCellMatrix();
@@ -445,7 +487,7 @@
// Assemble cell contribution into PETSc Matrix
jacobianVisitor.clear();
- PetscErrorCode err = updateOperator(jacobianMat, *submesh->getSieve(),
+ PetscErrorCode err = updateOperator(jacobianMat, *sieveSubMesh->getSieve(),
jacobianVisitor, *c_iter,
&_cellMatrix[0], ADD_VALUES);
CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -27,10 +27,13 @@
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
+//#define PRECOMPUTE_GEOMETRY
+
// ----------------------------------------------------------------------
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
typedef pylith::topology::SubMesh::RealSection SubRealSection;
typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
// ----------------------------------------------------------------------
// Default constructor.
@@ -79,20 +82,21 @@
const int numCorners = _quadrature->numBasis();
// Get 'surface' cells (1 dimension lower than top-level cells)
- const ALE::Obj<SieveSubMesh>& submesh = _boundaryMesh->sieveMesh();
- assert(!submesh.isNull());
+ const ALE::Obj<SieveSubMesh>& subSieveMesh = _boundaryMesh->sieveMesh();
+ assert(!subSieveMesh.isNull());
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- submesh->heightStratum(1);
+ subSieveMesh->heightStratum(1);
assert(!cells.isNull());
+ const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
// Make sure surface cells are compatible with quadrature.
- const int boundaryDepth = submesh->depth()-1; // depth of bndry cells
- for (SieveSubMesh::label_sequence::iterator c_iter=cells->begin();
+ const int boundaryDepth = subSieveMesh->depth()-1; // depth of bndry cells
+ for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
++c_iter) {
const int cellNumCorners =
- submesh->getNumCellCorners(*c_iter, boundaryDepth);
+ subSieveMesh->getNumCellCorners(*c_iter, boundaryDepth);
if (numCorners != cellNumCorners) {
std::ostringstream msg;
msg << "Quadrature is incompatible with cell for Neumann traction "
@@ -123,7 +127,6 @@
double_array jacobian(jacobianSize);
double jacobianDet = 0;
double_array orientation(orientationSize);
- double_array cellVertices(numCorners*spaceDim);
// Set names based on dimension of problem.
// 1-D problem = {'normal-traction'}
@@ -165,11 +168,13 @@
double_array cellTractionsGlobal(fiberDim);
// Get sections.
+ double_array coordinatesCell(numCorners*spaceDim);
const ALE::Obj<RealSection>& coordinates =
- submesh->getRealSection("coordinates");
+ subSieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- numCorners*spaceDim);
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
const ALE::Obj<SubRealSection>& tractSection = _tractions->section();
assert(!tractSection.isNull());
@@ -182,21 +187,28 @@
// Compute quadrature information
_quadrature->initializeGeometry();
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->computeGeometry(*_boundaryMesh, cells);
+#endif
// Loop over cells in boundary mesh, compute orientations, and then
// compute corresponding traction vector in global coordinates
// (store values in _tractionGlobal).
- for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
+ for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
c_iter != cellsEnd;
++c_iter) {
+ // Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ subSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
const double_array& quadPtsNondim = _quadrature->quadPts();
quadPtsGlobal = quadPtsNondim;
_normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
lengthScale);
- coordsVisitor.clear();
- submesh->restrictClosure(*c_iter, coordsVisitor);
cellTractionsGlobal = 0.0;
for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts;
@@ -219,10 +231,14 @@
// Compute Jacobian and determinant at quadrature point, then get
// orientation.
memcpy(&quadPtRef[0], &quadPtsRef[iRef], cellDim*sizeof(double));
- memcpy(&cellVertices[0], coordsVisitor.getValues(),
- cellVertices.size()*sizeof(double));
- cellGeometry.jacobian(&jacobian, &jacobianDet, cellVertices, quadPtRef);
+#if defined(PRECOMPUTE_GEOMETRY)
+ coordsVisitor.clear();
+ subSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+#endif
+ cellGeometry.jacobian(&jacobian, &jacobianDet,
+ coordinatesCell, quadPtRef);
cellGeometry.orientation(&orientation, jacobian, jacobianDet, up);
+ assert(jacobianDet > 0.0);
orientation /= jacobianDet;
// Rotate traction vector from local coordinate system to global
@@ -265,10 +281,10 @@
double_array tractionsCell(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveSubMesh>& submesh = _boundaryMesh->sieveMesh();
- assert(!submesh.isNull());
+ const ALE::Obj<SieveSubMesh>& subSieveMesh = _boundaryMesh->sieveMesh();
+ assert(!subSieveMesh.isNull());
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- submesh->heightStratum(1);
+ subSieveMesh->heightStratum(1);
assert(!cells.isNull());
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
@@ -279,12 +295,25 @@
topology::SubMesh::UpdateAddVisitor residualVisitor(*residualSection,
&_cellVector[0]);
+#if !defined(PRECOMPUTE_GEOMETRY)
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ subSieveMesh->getRealSection("coordinates");
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
+#endif
+
// Loop over faces and integrate contribution from each face
for (SieveSubMesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
++c_iter) {
- // Get geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ subSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
// Reset element vector to zero
_resetCellVector();
@@ -312,7 +341,7 @@
} // for
// Assemble cell contribution into field
residualVisitor.clear();
- submesh->updateAdd(*c_iter, residualVisitor);
+ subSieveMesh->updateAdd(*c_iter, residualVisitor);
PetscLogFlops(numQuadPts*(1+numBasis*(1+numBasis*(1+2*spaceDim))));
} // for
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -39,6 +39,8 @@
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
+//#define PRECOMPUTE_GEOMETRY
+
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
@@ -145,7 +147,9 @@
const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
_quadrature->initializeGeometry();
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->computeGeometry(*_faultMesh, cells);
+#endif
// Compute orientation at vertices in fault mesh.
_calcOrientation(upDir, normalDir);
@@ -196,7 +200,7 @@
// Allocate vectors for cell values
double_array orientationCell(numConstraintVert*orientationSize);
double_array stiffnessCell(numConstraintVert);
- double_array solutionCell(numCorners*spaceDim);
+ double_array dispTCell(numCorners*spaceDim);
double_array residualCell(numCorners*spaceDim);
// Tributary area for the current for each vertex.
@@ -242,25 +246,40 @@
topology::Mesh::RestrictVisitor areaVisitor(*areaSection,
areaCell.size(), &areaCell[0]);
- topology::Field<topology::Mesh>& solution = fields->solution();
- const ALE::Obj<RealSection>& solutionSection = solution.section();
- assert(!solutionSection.isNull());
- topology::Mesh::RestrictVisitor solutionVisitor(*solutionSection,
- solutionCell.size(),
- &solutionCell[0]);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ const ALE::Obj<RealSection>& dispTSection = dispT.section();
+ assert(!dispTSection.isNull());
+ topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
+ dispTCell.size(),
+ &dispTCell[0]);
const ALE::Obj<RealSection>& residualSection = residual.section();
topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
&residualCell[0]);
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ faultSieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
for (SieveMesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
c_iter != cellsCohesiveEnd;
++c_iter) {
const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
- _quadrature->retrieveGeometry(c_fault);
areaVertexCell = 0.0;
residualCell = 0.0;
+ // Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(c_fault);
+#else
+ coordsVisitor.clear();
+ faultSieveMesh->restrictClosure(c_fault, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, c_fault);
+#endif
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
const double_array& jacobianDet = _quadrature->jacobianDet();
@@ -286,9 +305,9 @@
areaVisitor.clear();
faultSieveMesh->restrictClosure(c_fault, areaVisitor);
- // Get solution at cohesive cell's vertices.
- solutionVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, solutionVisitor);
+ // Get disp(t) at cohesive cell's vertices.
+ dispTVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTVisitor);
for (int iConstraint=0; iConstraint < numConstraintVert; ++iConstraint) {
// Blocks in cell matrix associated with normal cohesive
@@ -311,7 +330,7 @@
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int kDim=0; kDim < spaceDim; ++kDim)
residualCell[indexI*spaceDim+iDim] -=
- solutionCell[indexK*spaceDim+kDim] *
+ dispTCell[indexK*spaceDim+kDim] *
-orientationVertex[kDim*spaceDim+iDim] * wt;
} // for
@@ -319,7 +338,7 @@
for (int jDim=0; jDim < spaceDim; ++jDim) {
for (int kDim=0; kDim < spaceDim; ++kDim)
residualCell[indexJ*spaceDim+jDim] -=
- solutionCell[indexK*spaceDim+kDim] *
+ dispTCell[indexK*spaceDim+kDim] *
orientationVertex[kDim*spaceDim+jDim] * wt;
} // for
} // for
@@ -333,7 +352,7 @@
std::cout << " slip["<<i<<"]: " << cellSlip[i] << std::endl;
}
for(int i = 0; i < numCorners*spaceDim; ++i) {
- std::cout << " soln["<<i<<"]: " << solutionCell[i] << std::endl;
+ std::cout << " soln["<<i<<"]: " << dispTCell[i] << std::endl;
}
for(int i = 0; i < numCorners*spaceDim; ++i) {
std::cout << " v["<<i<<"]: " << residualCell[i] << std::endl;
@@ -743,9 +762,9 @@
} else if (0 == strcasecmp("traction_change", name)) {
assert(0 != fields);
- const topology::Field<topology::Mesh>& solution = fields->solution();
+ const topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
_allocateBufferVectorField();
- _calcTractionsChange(_bufferVectorField, solution);
+ _calcTractionsChange(_bufferVectorField, dispT);
_bufferVectorField->label("traction_change");
return *_bufferVectorField;
@@ -828,7 +847,7 @@
orientation.allocate();
orientation.zero();
- // Get fault cells (1 dimension lower than top-level cells)
+ // Get fault cells.
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
faultSieveMesh->heightStratum(0);
assert(!cells.isNull());
@@ -1047,7 +1066,6 @@
assert(quadWts.size() == numQuadPts);
double jacobianDet = 0;
double_array areaCell(numBasis);
- double_array verticesCell(numBasis*spaceDim);
// Get vertices in fault mesh.
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
@@ -1068,6 +1086,14 @@
assert(!areaSection.isNull());
topology::Mesh::UpdateAddVisitor areaVisitor(*areaSection, &areaCell[0]);
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ faultSieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
faultSieveMesh->heightStratum(0);
assert(!cells.isNull());
@@ -1075,12 +1101,20 @@
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
// Loop over cells in fault mesh, compute area
- for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
+ for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
c_iter != cellsEnd;
++c_iter) {
- _quadrature->retrieveGeometry(*c_iter);
areaCell = 0.0;
+ // Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
+
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
const double_array& jacobianDet = _quadrature->jacobianDet();
@@ -1115,14 +1149,14 @@
void
pylith::faults::FaultCohesiveKin::_calcTractionsChange(
topology::Field<topology::SubMesh>* tractions,
- const topology::Field<topology::Mesh>& solution)
+ const topology::Field<topology::Mesh>& dispT)
{ // _calcTractionsChange
assert(0 != tractions);
assert(0 != _faultMesh);
assert(0 != _fields);
// Get vertices from mesh of domain.
- const ALE::Obj<SieveMesh>& sieveMesh = solution.mesh().sieveMesh();
+ const ALE::Obj<SieveMesh>& sieveMesh = dispT.mesh().sieveMesh();
assert(!sieveMesh.isNull());
const ALE::Obj<SieveMesh::label_sequence>& vertices =
sieveMesh->depthStratum(0);
@@ -1145,8 +1179,8 @@
const ALE::Obj<RealSection>& areaSection =
_fields->get("area").section();
assert(!areaSection.isNull());
- const ALE::Obj<RealSection>& solutionSection = solution.section();
- assert(!solutionSection.isNull());
+ const ALE::Obj<RealSection>& dispTSection = dispT.section();
+ assert(!dispTSection.isNull());
const int numFaultVertices = fvertices->size();
Mesh::renumbering_type& renumbering = faultSieveMesh->getRenumbering();
@@ -1196,13 +1230,13 @@
if (renumbering.find(*v_iter) != renumberingEnd) {
const int vertexMesh = *v_iter;
const int vertexFault = renumbering[*v_iter];
- assert(fiberDim == solutionSection->getFiberDimension(vertexMesh));
+ assert(fiberDim == dispTSection->getFiberDimension(vertexMesh));
assert(fiberDim == tractionsSection->getFiberDimension(vertexFault));
assert(1 == stiffnessSection->getFiberDimension(vertexFault));
assert(1 == areaSection->getFiberDimension(vertexFault));
- const double* solutionVertex = solutionSection->restrictPoint(vertexMesh);
- assert(0 != solutionVertex);
+ const double* dispTVertex = dispTSection->restrictPoint(vertexMesh);
+ assert(0 != dispTVertex);
const double* stiffnessVertex = stiffnessSection->restrictPoint(vertexFault);
assert(0 != stiffnessVertex);
const double* areaVertex = areaSection->restrictPoint(vertexFault);
@@ -1210,7 +1244,7 @@
const double scale = stiffnessVertex[0] / areaVertex[0];
for (int i=0; i < fiberDim; ++i)
- tractionsVertex[i] = solutionVertex[i] * scale;
+ tractionsVertex[i] = dispTVertex[i] * scale;
tractionsSection->updatePoint(vertexFault, &tractionsVertex[0]);
} // if
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -34,6 +34,8 @@
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
+//#define PRECOMPUTE_GEOMETRY
+
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
@@ -180,6 +182,14 @@
topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
&_cellVector[0]);
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
const double gravityScale =
@@ -199,7 +209,13 @@
++c_iter) {
// Compute geometry information for current cell
_logger->eventBegin(geometryEvent);
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
_logger->eventEnd(geometryEvent);
// Get state variables for cell.
@@ -332,6 +348,14 @@
(int) pow(sieveMesh->getSieve()->getMaxConeSize(),
sieveMesh->depth())*spaceDim);
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
_logger->eventEnd(setupEvent);
// Loop over cells
@@ -340,7 +364,13 @@
++c_iter) {
// Compute geometry information for current cell
_logger->eventBegin(geometryEvent);
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
_logger->eventEnd(geometryEvent);
// Get state variables for cell.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -36,6 +36,8 @@
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
+//#define PRECOMPUTE_GEOMETRY
+
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
@@ -185,6 +187,14 @@
topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
&_cellVector[0]);
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
const double gravityScale =
@@ -199,7 +209,13 @@
++c_iter) {
// Compute geometry information for current cell
_logger->eventBegin(geometryEvent);
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
_logger->eventEnd(geometryEvent);
// Get state variables for cell.
@@ -397,6 +413,14 @@
(int) pow(sieveMesh->getSieve()->getMaxConeSize(),
sieveMesh->depth())*spaceDim);
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
_logger->eventEnd(setupEvent);
// Loop over cells
@@ -405,7 +429,13 @@
++c_iter) {
// Compute geometry information for current cell
_logger->eventBegin(geometryEvent);
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
_logger->eventEnd(geometryEvent);
// Get state variables for cell.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -31,6 +31,8 @@
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
+//#define PRECOMPUTE_GEOMETRY
+
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
@@ -102,7 +104,9 @@
// Compute geometry for quadrature operations.
_quadrature->initializeGeometry();
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->computeGeometry(mesh, cells);
+#endif
// Initialize material.
_material->initialize(mesh, _quadrature);
@@ -200,18 +204,34 @@
const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
// Get fields
- const topology::Field<topology::Mesh>& solution = fields->solution();
- const ALE::Obj<RealSection>& disp = solution.section();
- assert(!disp.isNull());
- topology::Mesh::RestrictVisitor dispVisitor(*disp,
+ const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
+ const ALE::Obj<RealSection>& dispSection = disp.section();
+ assert(!dispSection.isNull());
+ topology::Mesh::RestrictVisitor dispVisitor(*dispSection,
dispCell.size(), &dispCell[0]);
+#if !defined(PRECOMPUTE_GEOMETRY)
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+#endif
+
// Loop over cells
for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
++c_iter) {
// Retrieve geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
// Restrict input fields to cell
dispVisitor.clear();
@@ -410,12 +430,22 @@
const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
// Get field
- const topology::Field<topology::Mesh>& solution = fields->solution();
- const ALE::Obj<RealSection>& disp = solution.section();
- assert(!disp.isNull());
- topology::Mesh::RestrictVisitor dispVisitor(*disp,
+ const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
+ const ALE::Obj<RealSection>& dispSection = disp.section();
+ assert(!dispSection.isNull());
+ topology::Mesh::RestrictVisitor dispVisitor(*dispSection,
dispCell.size(), &dispCell[0]);
+#if !defined(PRECOMPUTE_GEOMETRY)
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+#endif
+
const ALE::Obj<RealSection>& fieldSection = field->section();
assert(!fieldSection.isNull());
@@ -424,7 +454,13 @@
c_iter != cellsEnd;
++c_iter) {
// Retrieve geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
// Restrict input fields to cell
dispVisitor.clear();
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -28,6 +28,8 @@
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
+//#define PRECOMPUTE_GEOMETRY
+
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
@@ -335,15 +337,26 @@
// Get quadrature information
const int numQuadPts = quadrature->numQuadPts();
const int spaceDim = quadrature->spaceDim();
+ const int numBasis = quadrature->numBasis();
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+#if !defined(PRECOMPUTE_GEOMETRY)
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+#endif
+
// Create arrays for querying
const int tensorSize = _tensorSize;
double_array quadPtsGlobal(numQuadPts*spaceDim);
double_array stressCell(numQuadPts*tensorSize);
// Get cells associated with material
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->getLabelStratum("material-id", id());
assert(!cells.isNull());
@@ -408,7 +421,13 @@
c_iter != cellsEnd;
++c_iter) {
// Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
// Dimensionalize coordinates for querying
const double_array& quadPtsNonDim = quadrature->quadPts();
@@ -461,15 +480,26 @@
// Get quadrature information
const int numQuadPts = quadrature->numQuadPts();
const int spaceDim = quadrature->spaceDim();
+ const int numBasis = quadrature->numBasis();
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+#if !defined(PRECOMPUTE_GEOMETRY)
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+#endif
+
// Create arrays for querying
const int tensorSize = _tensorSize;
double_array quadPtsGlobal(numQuadPts*spaceDim);
double_array strainCell(numQuadPts*tensorSize);
// Get cells associated with material
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->getLabelStratum("material-id", id());
assert(!cells.isNull());
@@ -534,7 +564,13 @@
c_iter != cellsEnd;
++c_iter) {
// Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
// Dimensionalize coordinates for querying
const double_array& quadPtsNonDim = quadrature->quadPts();
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -28,9 +28,12 @@
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
+//#define PRECOMPUTE_GEOMETRY
+
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
// ----------------------------------------------------------------------
// Default constructor.
@@ -104,6 +107,7 @@
// Get quadrature information
const int numQuadPts = quadrature->numQuadPts();
+ const int numBasis = quadrature->numBasis();
const int spaceDim = quadrature->spaceDim();
// Get cells associated with material
@@ -127,6 +131,14 @@
const ALE::Obj<RealSection>& propertiesSection = _properties->section();
assert(!propertiesSection.isNull());
+#if !defined(PRECOMPUTE_GEOMETRY)
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
+#endif
+
// Create arrays for querying.
const int numDBProperties = _metadata.numDBProperties();
double_array quadPtsGlobal(numQuadPts*spaceDim);
@@ -180,7 +192,13 @@
c_iter != cellsEnd;
++c_iter) {
// Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
const double_array& quadPtsNonDim = quadrature->quadPts();
quadPtsGlobal = quadPtsNonDim;
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -39,8 +39,7 @@
const int numCells,
const int numCorners,
const int meshDim,
- const bool interpolate,
- const spatialdata::units::Nondimensional& normalizer)
+ const bool interpolate)
{ // buildMesh
assert(0 != mesh);
@@ -175,10 +174,6 @@
logger.stagePop();
} // if/else
- const double lengthScale = normalizer.lengthScale();
- normalizer.nondimensionalize(&(*coordinates)[0], coordinates->size(),
- lengthScale);
-
logger.stagePush("MeshCoordinates");
ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
&(*coordinates)[0]);
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.hh 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshBuilder.hh 2009-05-18 23:02:05 UTC (rev 15011)
@@ -53,8 +53,7 @@
const int numCells,
const int numCorners,
const int meshDim,
- const bool interpolate,
- const spatialdata::units::Nondimensional& normalizer);
+ const bool interpolate);
/** Build fault mesh topology.
*
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -17,6 +17,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/utils/array.hh" // USES double_array, int_array
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include "Selection.hh" // USES boundary()
@@ -34,7 +35,6 @@
// Constructor
pylith::meshio::MeshIO::MeshIO(void) :
_mesh(0),
- _normalizer(new spatialdata::units::Nondimensional),
_debug(false),
_interpolate(false)
{ // constructor
@@ -44,7 +44,6 @@
// Destructor
pylith::meshio::MeshIO::~MeshIO(void)
{ // destructor
- delete _normalizer; _normalizer = 0;
} // destructor
// ----------------------------------------------------------------------
@@ -56,18 +55,6 @@
} // getMeshDim
// ----------------------------------------------------------------------
-// Set scales used to nondimensionalize mesh.
-void
-pylith::meshio::MeshIO::normalizer(
- const spatialdata::units::Nondimensional& dim)
-{ // normalizer
- if (0 == _normalizer)
- _normalizer = new spatialdata::units::Nondimensional(dim);
- else
- *_normalizer = dim;
-} // normalizer
-
-// ----------------------------------------------------------------------
// Read mesh from file.
void
pylith::meshio::MeshIO::read(topology::Mesh* mesh)
@@ -77,6 +64,7 @@
_mesh = mesh;
_mesh->debug(_debug);
_read();
+
_mesh = 0;
} // read
@@ -115,6 +103,8 @@
const SieveMesh::label_sequence::iterator verticesEnd =
vertices->end();
const ALE::Obj<RealSection>& coordsField =
+ sieveMesh->hasRealSection("coordinates_dimensioned") ?
+ sieveMesh->getRealSection("coordinates_dimensioned") :
sieveMesh->getRealSection("coordinates");
assert(!coordsField.isNull());
@@ -128,7 +118,7 @@
coordinates->resize(size);
int i = 0;
- for(SieveMesh::label_sequence::iterator v_iter=verticesBegin;
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter) {
const RealSection::value_type *vertexCoords =
@@ -136,11 +126,6 @@
for (int iDim=0; iDim < *spaceDim; ++iDim)
(*coordinates)[i++] = vertexCoords[iDim];
} // for
-
- assert(0 != _normalizer);
- const double lengthScale = _normalizer->lengthScale();
- _normalizer->dimensionalize(&(*coordinates)[0], coordinates->size(),
- lengthScale);
} // _getVertices
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh 2009-05-18 23:02:05 UTC (rev 15011)
@@ -76,12 +76,6 @@
*/
bool interpolate(void) const;
- /** Set scales used to nondimensionalize mesh.
- *
- * @param dim Nondimensionalizer.
- */
- void normalizer(const spatialdata::units::Nondimensional& dim);
-
/** Read mesh from file.
*
* @param mesh PyLith finite-element mesh.
@@ -193,7 +187,6 @@
protected :
topology::Mesh* _mesh; ///< Pointer to finite-element mesh.
- spatialdata::units::Nondimensional* _normalizer; ///< Nondimensionalizer.
bool _debug; ///< True to turn of mesh debugging output.
bool _interpolate; ///< True if building intermediate topology elements.
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -54,8 +54,6 @@
void
pylith::meshio::MeshIOAscii::_read(void)
{ // _read
- assert(0 != _normalizer);
-
MPI_Comm comm = _mesh->comm();
int rank = 0;
int meshDim = 0;
@@ -141,7 +139,7 @@
// Can now build mesh
MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim,
cells, numCells, numCorners, meshDim,
- _interpolate, *_normalizer);
+ _interpolate);
_setMaterials(materialIds);
builtMesh = true;
} // if
@@ -169,7 +167,7 @@
} else {
MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim,
cells, numCells, numCorners, meshDim,
- _interpolate, *_normalizer);
+ _interpolate);
_setMaterials(materialIds);
} // if/else
_distributeGroups();
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIOCubit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIOCubit.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIOCubit.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -79,7 +79,7 @@
_orientCells(&cells, numCells, numCorners, meshDim);
MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim,
cells, numCells, numCorners, meshDim,
- _interpolate, *_normalizer);
+ _interpolate);
_setMaterials(materialIds);
_readGroups(ncfile);
@@ -97,7 +97,7 @@
} else {
MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim,
cells, numCells, numCorners, meshDim,
- _interpolate, *_normalizer);
+ _interpolate);
_setMaterials(materialIds);
}
_distributeGroups();
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -75,7 +75,7 @@
}
MeshBuilder::buildMesh(_mesh, &coordinates, numVertices, spaceDim,
cells, numCells, numCorners, meshDim,
- _interpolate, *_normalizer);
+ _interpolate);
_setMaterials(materialIds);
if (0 == rank) {
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -91,14 +91,14 @@
// Reform system residual.
void
pylith::problems::Formulation::reformResidual(const PetscVec* tmpResidualVec,
- const PetscVec* tmpSolveSolnVec)
+ const PetscVec* tmpSolutionVec)
{ // reformResidual
assert(0 != _fields);
// Update section view of field.
- if (0 != tmpSolveSolnVec) {
- topology::Field<topology::Mesh>& solveSoln = _fields->solveSoln();
- solveSoln.scatterVectorToSection(*tmpSolveSolnVec);
+ if (0 != tmpSolutionVec) {
+ topology::Field<topology::Mesh>& solution = _fields->solution();
+ solution.scatterVectorToSection(*tmpSolutionVec);
} // if
// Set residual to zero.
@@ -144,15 +144,15 @@
// ----------------------------------------------------------------------
// Reform system Jacobian.
void
-pylith::problems::Formulation::reformJacobian(const PetscVec* tmpSolveSolnVec)
+pylith::problems::Formulation::reformJacobian(const PetscVec* tmpSolutionVec)
{ // reformJacobian
assert(0 != _jacobian);
assert(0 != _fields);
// Update section view of field.
- if (0 != tmpSolveSolnVec) {
- topology::Field<topology::Mesh>& solveSoln = _fields->solveSoln();
- solveSoln.scatterVectorToSection(*tmpSolveSolnVec);
+ if (0 != tmpSolutionVec) {
+ topology::Field<topology::Mesh>& solution = _fields->solution();
+ solution.scatterVectorToSection(*tmpSolutionVec);
} // if
// Set residual to zero.
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -83,19 +83,19 @@
// Solve the system.
void
pylith::problems::SolverNonlinear::solve(
- topology::Field<topology::Mesh>* solveSoln,
+ topology::Field<topology::Mesh>* solution,
const topology::Jacobian& jacobian,
const topology::Field<topology::Mesh>& residual)
{ // solve
- assert(0 != solveSoln);
+ assert(0 != solution);
PetscErrorCode err = 0;
- const PetscVec solveSolnVec = solveSoln->vector();
- err = SNESSolve(_snes, PETSC_NULL, solveSolnVec); CHECK_PETSC_ERROR(err);
+ const PetscVec solutionVec = solution->vector();
+ err = SNESSolve(_snes, PETSC_NULL, solutionVec); CHECK_PETSC_ERROR(err);
// Update section view of field.
- solveSoln->scatterVectorToSection();
+ solution->scatterVectorToSection();
} // solve
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/topology/Mesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Mesh.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/topology/Mesh.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -15,6 +15,8 @@
#include "Mesh.hh" // implementation of class methods
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+#include "pylith/utils/array.hh" // USES double_array
// ----------------------------------------------------------------------
// Default constructor
@@ -68,12 +70,57 @@
} // coordsys
// ----------------------------------------------------------------------
-// Initialize the finite-element mesh.
+// Nondimensionalizer the finite-element mesh.
void
-pylith::topology::Mesh::initialize(void)
+pylith::topology::Mesh::nondimensionalize(const spatialdata::units::Nondimensional& normalizer)
{ // initialize
-} // initialize
+ // Get coordinates (currently dimensioned).
+ assert(!_mesh.isNull());
+ const ALE::Obj<RealSection>& coordsSection =
+ _mesh->getRealSection("coordinates");
+ assert(!coordsSection.isNull());
+ // Get field for dimensioned coordinates.
+ const ALE::Obj<RealSection>& coordsDimSection =
+ _mesh->getRealSection("coordinates_dimensioned");
+ assert(!coordsDimSection.isNull());
+ coordsDimSection->setAtlas(coordsSection->getAtlas());
+ coordsDimSection->allocateStorage();
+ coordsDimSection->setBC(coordsSection->getBC());
+
+ const double lengthScale = normalizer.lengthScale();
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ _mesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const SieveMesh::label_sequence::iterator verticesBegin =
+ vertices->begin();
+ const SieveMesh::label_sequence::iterator verticesEnd =
+ vertices->end();
+
+ assert(0 != _coordsys);
+ const int spaceDim = _coordsys->spaceDim();
+ double_array coordsVertex(spaceDim);
+ double_array coordsDimVertex(spaceDim);
+
+ int i = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter) {
+ coordsSection->restrictPoint(*v_iter,
+ &coordsVertex[0], coordsVertex.size());
+ coordsDimSection->restrictPoint(*v_iter, &coordsDimVertex[0],
+ coordsDimVertex.size());
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ coordsDimVertex[iDim] = coordsVertex[iDim];
+
+ // Nondimensionalize original coordinates.
+ normalizer.nondimensionalize(&coordsVertex[0], spaceDim, lengthScale);
+
+ coordsSection->updatePoint(*v_iter, &coordsVertex[0]);
+ coordsDimSection->updatePoint(*v_iter, &coordsDimVertex[0]);
+ } // for
+} // nondimensionalize
+
// ----------------------------------------------------------------------
// Return the names of all vertex groups.
void
Modified: short/3D/PyLith/trunk/libsrc/topology/Mesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Mesh.hh 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/topology/Mesh.hh 2009-05-18 23:02:05 UTC (rev 15011)
@@ -25,6 +25,7 @@
// Include directives ---------------------------------------------------
#include "topologyfwd.hh" // forward declarations
#include "spatialdata/geocoords/geocoordsfwd.hh" // forward declarations
+#include "spatialdata/units/unitsfwd.hh" // forward declarations
#include <petscmesh.hh> // HASA ALE::IMesh
@@ -149,8 +150,11 @@
*/
const MPI_Comm comm(void) const;
- /// Initialize the finite-element mesh.
- void initialize(void);
+ /** Initialize the finite-element mesh.
+ *
+ * @param normalizer Nondimensionalizer.
+ */
+ void nondimensionalize(const spatialdata::units::Nondimensional& normalizer);
/** Print mesh to stdout.
*
Modified: short/3D/PyLith/trunk/libsrc/topology/SolutionFields.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/SolutionFields.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/topology/SolutionFields.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -20,8 +20,7 @@
// Default constructor.
pylith::topology::SolutionFields::SolutionFields(const Mesh& mesh) :
Fields<Field<Mesh> >(mesh),
- _solutionName(""),
- _solveSolnName("")
+ _solutionName("")
{ // constructor
} // constructor
@@ -68,76 +67,5 @@
return get(_solutionName.c_str());
} // solution
-// ----------------------------------------------------------------------
-// Set field used in the solve.
-void
-pylith::topology::SolutionFields::solveSolnName(const char* name)
-{ // solveSolnName
- map_type::const_iterator iter = _fields.find(name);
- if (iter == _fields.end()) {
- std::ostringstream msg;
- msg << "Cannot use unknown field '" << name
- << "' when setting name of field used in solve.";
- throw std::runtime_error(msg.str());
- } // if
- _solveSolnName = name;
-} // solveSolnName
-// ----------------------------------------------------------------------
-// Get solveSoln field.
-const pylith::topology::Field<pylith::topology::Mesh>&
-pylith::topology::SolutionFields::solveSoln(void) const
-{ // solveSoln
- if (_solveSolnName == "")
- throw std::runtime_error("Cannot retrieve solve field. Name of solve "
- "field has not been specified.");
- return get(_solveSolnName.c_str());
-} // solveSoln
-
-// ----------------------------------------------------------------------
-// Get solveSoln field.
-pylith::topology::Field<pylith::topology::Mesh>&
-pylith::topology::SolutionFields::solveSoln(void)
-{ // solveSoln
- if (_solveSolnName == "")
- throw std::runtime_error("Cannot retrieve solve field. Name of solve "
- "field has not been specified.");
- return get(_solveSolnName.c_str());
-} // solveSoln
-
-// ----------------------------------------------------------------------
-// Create history manager for a subset of the managed fields.
-void
-pylith::topology::SolutionFields::createHistory(const char* const* fields,
- const int size)
-{ // createHistory
- if (size > 0 && 0 != fields) {
- _history.resize(size);
- for (int i=0; i < size; ++i) {
- map_type::const_iterator iter = _fields.find(fields[i]);
- if (iter == _fields.end()) {
- std::ostringstream msg;
- msg << "Cannot use unknown field '" << fields[i]
- << "' when creating history.";
- throw std::runtime_error(msg.str());
- } // if
- _history[i] = fields[i];
- } // for
- } // if
-} // createHistory
-
-// ----------------------------------------------------------------------
-// Shift fields in history.
-void
-pylith::topology::SolutionFields::shiftHistory(void)
-{ // shiftHistory
- assert(_history.size() > 0);
- const int size = _history.size();
- Field<Mesh>* tmp = _fields[_history[size-1]];
- for (int i=size-1; i > 0; --i)
- _fields[_history[i]] = _fields[_history[i-1]];
- _fields[_history[0]] = tmp;
-} // shiftHistory
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/SolutionFields.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/SolutionFields.hh 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/libsrc/topology/SolutionFields.hh 2009-05-18 23:02:05 UTC (rev 15011)
@@ -46,56 +46,24 @@
/// Destructor.
~SolutionFields(void);
- /** Set name of field for the current solution.
+ /** Set name of field used in the solve.
*
* @param name Name of field that holds the solution.
*/
void solutionName(const char* name);
- /** Get current solution field.
+ /** Get field used in solve.
*
* @returns Solution field.
*/
const Field<Mesh>& solution(void) const;
- /** Get current solution field.
+ /** Get field used in solve.
*
* @returns Solution field.
*/
Field<Mesh>& solution(void);
- /** Set name of field that will be used in the solve.
- *
- * @param name Name of field used in the solve.
- */
- void solveSolnName(const char* name);
-
- /** Get field used in the solve.
- *
- * @returns Field used in the solve.
- */
- const Field<Mesh>& solveSoln(void) const;
-
- /** Get field used in the solve.
- *
- * @returns Field used in the solve.
- */
- Field<Mesh>& solveSoln(void);
-
- /** Create history manager for a subset of the managed fields.
- *
- * @param fields Fields in history (first is most recent).
- * @param size Number of fields in history.
- */
- void createHistory(const char* const* fields,
- const int size);
-
- /** Shift fields in history. Handles to fields are shifted so that
- * the most recent values become associated with the second most
- * recent item in the history, etc.
- */
- void shiftHistory(void);
-
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
@@ -105,14 +73,6 @@
/// problem.
std::string _solutionName;
- /// Name of field used in the solve (solution of the solve). This
- /// may be an increment that is applied to the "working" solutio to
- /// form the complete solution.
- std::string _solveSolnName;
-
- /// History manager for a subset of the fields
- std::vector<std::string> _history;
-
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/modulesrc/meshio/MeshIOObj.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/meshio/MeshIOObj.i 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/modulesrc/meshio/MeshIOObj.i 2009-05-18 23:02:05 UTC (rev 15011)
@@ -71,12 +71,6 @@
*/
bool interpolate(void) const;
- /** Set scales used to nondimensionalize mesh.
- *
- * @param dim Nondimensionalizer.
- */
- void normalizer(const spatialdata::units::Nondimensional& dim);
-
/** Read mesh from file.
*
* @param mesh PyLith finite-element mesh.
Modified: short/3D/PyLith/trunk/modulesrc/topology/Mesh.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/Mesh.i 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/modulesrc/topology/Mesh.i 2009-05-18 23:02:05 UTC (rev 15011)
@@ -106,9 +106,12 @@
*/
const MPI_Comm comm(void) const;
- /// Initialize the finite-element mesh.
- void initialize(void);
-
+ /** Initialize the finite-element mesh.
+ *
+ * @param normalizer Nondimensionalizer.
+ */
+ void nondimensionalize(const spatialdata::units::Nondimensional& normalizer);
+
/** Print mesh to stdout.
*
* @param label Label for mesh.
Modified: short/3D/PyLith/trunk/modulesrc/topology/SolutionFields.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/SolutionFields.i 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/modulesrc/topology/SolutionFields.i 2009-05-18 23:02:05 UTC (rev 15011)
@@ -54,43 +54,6 @@
*/
Field<Mesh>& solution(void);
- /** Set name of field that will be used in the solve.
- *
- * @param name Name of field used in the solve.
- */
- void solveSolnName(const char* name);
-
- /** Get field used in the solve.
- *
- * @returns Field used in the solve.
- */
- const Field<Mesh>& solveSoln(void) const;
-
- /** Get field used in the solve.
- *
- * @returns Field used in the solve.
- */
- Field<Mesh>& solveSoln(void);
-
- /** Create history manager for a subset of the managed fields.
- *
- * @param fields Fields in history (first is most recent).
- * @param size Number of fields in history.
- */
- %apply(const char* const* string_list, const int list_len){
- (const char* const* fields,
- const int size)
- };
- void createHistory(const char* const* fields,
- const int size);
- %clear(const char* const* fields, const int size);
-
- /** Shift fields in history. Handles to fields are shifted so that
- * the most recent values become associated with the second most
- * recent item in the history, etc.
- */
- void shiftHistory(void);
-
}; // SolutionFields
} // topology
Modified: short/3D/PyLith/trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/pylith/Makefile.am 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/Makefile.am 2009-05-18 23:02:05 UTC (rev 15011)
@@ -122,6 +122,7 @@
utils/CppData.py \
utils/EmptyBin.py \
utils/EventLogger.py \
+ utils/NullComponent.py \
utils/PetscComponent.py \
utils/PetscManager.py \
utils/VTKDataReader.py \
Modified: short/3D/PyLith/trunk/pylith/materials/ElasticMaterial.py
===================================================================
--- short/3D/PyLith/trunk/pylith/materials/ElasticMaterial.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/materials/ElasticMaterial.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -39,8 +39,7 @@
## Python object for managing FaultCohesiveKin facilities and properties.
##
## \b Properties
- ## @li \b use_initial_stress Use initial stress (true) or not (false).
- ## @li \b use_initial_strain Use initial strain (true) or not (false).
+ ## @li None
##
## \b Facilities
## @li \b output Output manager associated with fault data.
@@ -49,26 +48,20 @@
import pyre.inventory
- useInitialStress = pyre.inventory.bool("use_initial_stress", default=False)
- useInitialStress.meta['tip'] = "Use initial stress for material."
-
- useInitialStrain = pyre.inventory.bool("use_initial_strain", default=False)
- useInitialStrain.meta['tip'] = "Use initial strain for material."
-
from pylith.meshio.OutputMatElastic import OutputMatElastic
output = pyre.inventory.facility("output", family="output_manager",
factory=OutputMatElastic)
output.meta['tip'] = "Output manager for elastic material information."
- from spatialdata.spatialdb.SimpleDB import SimpleDB
+ from pylith.utils.NullComponent import NullComponent
dbInitialStress = pyre.inventory.facility("initial_stress_db",
family="spatial_database",
- factory=SimpleDB)
+ factory=NullComponent)
dbInitialStress.meta['tip'] = "Database for initial stress."
dbInitialStrain = pyre.inventory.facility("initial_strain_db",
family="spatial_database",
- factory=SimpleDB)
+ factory=NullComponent)
dbInitialStrain.meta['tip'] = "Database for initial strain."
@@ -90,9 +83,10 @@
"""
Material._configure(self)
self.output = self.inventory.output
- if self.inventory.useInitialStress:
+ from pylith.utils.NullComponent import NullComponent
+ if not isinstance(self.inventory.dbInitialStress, NullComponent):
self.dbInitialStress(self.inventory.dbInitialStress)
- if self.inventory.useInitialStrain:
+ if not isinstance(self.inventory.dbInitialStrain, NullComponent):
self.dbInitialStrain(self.inventory.dbInitialStrain)
return
Modified: short/3D/PyLith/trunk/pylith/materials/Material.py
===================================================================
--- short/3D/PyLith/trunk/pylith/materials/Material.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/materials/Material.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -51,7 +51,6 @@
## \b Properties
## @li \b id Material identifier (from mesh generator)
## @li \b name Name of material
- ## @li \b use_initial_state Use initial state (true) or not (false).
##
## \b Facilities
## @li \b db Database of material property parameters
@@ -60,9 +59,6 @@
import pyre.inventory
- useInitialState = pyre.inventory.bool("use_initial_state", default=False)
- useInitialState.meta['tip'] = "Use initial state for material."
-
id = pyre.inventory.int("id", default=0)
id.meta['tip'] = "Material identifier (from mesh generator)."
@@ -75,9 +71,10 @@
factory=SimpleDB)
dbProperties.meta['tip'] = "Database for physical property parameters."
+ from pylith.utils.NullComponent import NullComponent
dbInitialState = pyre.inventory.facility("initial_state_db",
family="spatial_database",
- factory=SimpleDB)
+ factory=NullComponent)
dbInitialState.meta['tip'] = "Database for initial state variables."
from pylith.feassemble.Quadrature import MeshQuadrature
@@ -149,7 +146,8 @@
self.id(self.inventory.id)
self.label(self.inventory.label)
self.dbProperties(self.inventory.dbProperties)
- if self.inventory.useInitialState:
+ from pylith.utils.NullComponent import NullComponent
+ if not isinstance(self.inventory.dbInitialState, NullComponent):
self.dbInitialState(self.inventory.dbInitialState)
self.quadrature = self.inventory.quadrature
Modified: short/3D/PyLith/trunk/pylith/meshio/CellFilter.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/CellFilter.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/meshio/CellFilter.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -35,7 +35,6 @@
Constructor.
"""
PetscComponent.__init__(self, name, facility="cellfilter")
- self.filter = None
return
@@ -57,9 +56,10 @@
def output_cell_filter():
"""
- Factory associated with MeshCellFilter.
+ Factory associated with CellFilter.
"""
- return CellFilter()
+ # Abstract object (so return None).
+ return None
# End of file
Modified: short/3D/PyLith/trunk/pylith/meshio/CellFilterAvgMesh.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/CellFilterAvgMesh.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/meshio/CellFilterAvgMesh.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -37,7 +37,6 @@
"""
CellFilter.__init__(self, name)
ModuleCellFilterAvg.__init__(self)
- self.filter = True
return
Modified: short/3D/PyLith/trunk/pylith/meshio/CellFilterAvgSubMesh.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/CellFilterAvgSubMesh.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/meshio/CellFilterAvgSubMesh.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -37,7 +37,6 @@
"""
CellFilter.__init__(self, name)
ModuleCellFilterAvg.__init__(self)
- self.filter = True
return
Modified: short/3D/PyLith/trunk/pylith/meshio/MeshIOObj.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/MeshIOObj.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/meshio/MeshIOObj.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -46,7 +46,6 @@
self._info.log("Reading finite-element mesh")
# Set flags
- self.normalizer(normalizer)
self.debug(debug)
self.interpolate(interpolate)
@@ -59,10 +58,11 @@
mesh = Mesh()
mesh.setComm(petsc_comm_world())
mesh.coordsys(self.coordsys)
- mesh.initialize()
# Read mesh
ModuleMeshIO.read(self, mesh)
+
+ mesh.nondimensionalize(normalizer)
return mesh
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputManager.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputManager.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputManager.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -18,6 +18,7 @@
## Factory: output_manager
from pylith.utils.PetscComponent import PetscComponent
+from pylith.utils.NullComponent import NullComponent
# OutputManager class
class OutputManager(PetscComponent):
@@ -59,16 +60,14 @@
factory=CSCart)
coordsys.meta['tip'] = "Coordinate system for output."
- from VertexFilter import VertexFilter
vertexFilter = pyre.inventory.facility("vertex_filter",
- factory=VertexFilter,
- family="output_vertex_filter")
+ family="output_vertex_filter",
+ factory=NullComponent)
vertexFilter.meta['tip'] = "Filter for vertex data."
- from CellFilter import CellFilter
cellFilter = pyre.inventory.facility("cell_filter",
- factory=CellFilter,
- family="output_cell_filter")
+ family="output_cell_filter",
+ factory=NullComponent)
cellFilter.meta['tip'] = "Filter for cell data."
@@ -140,7 +139,8 @@
raise ValueError, "Coordinate system for output is unknown."
self.coordsys.initialize()
- self.cellFilter.initialize(quadrature)
+ if not isinstance(self.inventory.cellFilter, NullComponent):
+ self.cellFilter.initialize(quadrature)
self.writer.initialize(normalizer)
self._logger.eventEnd(logEvent)
@@ -241,6 +241,7 @@
PetscComponent._configure(self)
return
+
def _createModuleObj(self):
"""
Create handle to C++ object.
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputManagerMesh.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputManagerMesh.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputManagerMesh.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -61,12 +61,14 @@
"""
Set members based using inventory.
"""
+ from pylith.utils.NullComponent import NullComponent
+
OutputManager._configure(self)
ModuleOutputManager.coordsys(self, self.inventory.coordsys)
ModuleOutputManager.writer(self, self.inventory.writer)
- if None != self.vertexFilter.filter:
+ if not isinstance(self.inventory.vertexFilter, NullComponent):
ModuleOutputManager.vertexFilter(self, self.inventory.vertexFilter)
- if None != self.cellFilter.filter:
+ if not isinstance(self.inventory.cellFilter, NullComponent):
ModuleOutputManager.cellFilter(self, self.inventory.cellFilter)
return
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputManagerSubMesh.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputManagerSubMesh.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputManagerSubMesh.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -60,12 +60,14 @@
"""
Set members based using inventory.
"""
+ from pylith.utils.NullComponent import NullComponent
+
OutputManager._configure(self)
ModuleOutputManager.coordsys(self, self.inventory.coordsys)
ModuleOutputManager.writer(self, self.inventory.writer)
- if None != self.vertexFilter.filter:
+ if not isinstance(self.inventory.vertexFilter, NullComponent):
ModuleOutputManager.vertexFilter(self, self.inventory.vertexFilter)
- if None != self.cellFilter.filter:
+ if not isinstance(self.inventory.cellFilter, NullComponent):
ModuleOutputManager.cellFilter(self, self.inventory.cellFilter)
return
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputSolnSubset.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputSolnSubset.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputSolnSubset.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -134,9 +134,10 @@
ModuleOutputSolnSubset.label(self, self.label)
ModuleOutputSolnSubset.coordsys(self, self.inventory.coordsys)
ModuleOutputSolnSubset.writer(self, self.inventory.writer)
- if None != self.vertexFilter.filter:
+ from pylith.utils.NullComponent import NullComponent
+ if not isinstance(self.inventory.vertexFilter, NullComponent):
ModuleOutputSolnSubset.vertexFilter(self, self.inventory.vertexFilter)
- if None != self.cellFilter.filter:
+ if not isinstance(self.inventory.cellFilter, NullComponent):
ModuleOutputSolnSubset.cellFilter(self, self.inventory.cellFilter)
return
Modified: short/3D/PyLith/trunk/pylith/meshio/VertexFilter.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/VertexFilter.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/meshio/VertexFilter.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -39,7 +39,6 @@
Constructor.
"""
PetscComponent.__init__(self, name, facility="vertexfilter")
- self.filter = None
return
@@ -63,7 +62,8 @@
"""
Factory associated with VertexFilter.
"""
- return VertexFilter()
+ # Abstract object (so return None).
+ return None
# End of file
Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -74,7 +74,7 @@
self.fields.add("disp(t-dt)", "displacement")
self.fields.add("residual", "residual")
self.fields.copyLayout("disp(t)")
- self.fields.solveSolnName("dispIncr(t->t+dt)")
+ self.fields.solutionName("dispIncr(t->t+dt)")
self._debug.log(resourceUsageString())
# Set fields to zero
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -182,13 +182,14 @@
self.fields = SolutionFields(self.mesh)
self._debug.log(resourceUsageString())
- if self.gravityField != None:
+ if not self.gravityField is None:
self._info.log("Initializing gravity field.")
self.gravityField.initialize()
self._info.log("Initializing integrators.")
for integrator in self.integratorsMesh + self.integratorsSubMesh:
- integrator.gravityField = self.gravityField
+ if not self.gravityField is None:
+ integrator.gravityField(self.gravityField)
integrator.initialize(totalTime, numTimeSteps, normalizer)
ModuleFormulation.meshIntegrators(self, self.integratorsMesh)
ModuleFormulation.submeshIntegrators(self, self.integratorsSubMesh)
Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -94,7 +94,7 @@
self.fields.add("dispIncr(t->t+dt)", "displacement increment")
self.fields.add("residual", "residual")
self.fields.copyLayout("disp(t)")
- self.fields.solveSolnName("dispIncr(t->t+dt)")
+ self.fields.solutionName("dispIncr(t->t+dt)")
# Set fields to zero
disp = self.fields.get("disp(t)")
@@ -178,7 +178,7 @@
"""
Advance to next time step.
"""
- dispIncr = self.fields.solveSoln()
+ dispIncr = self.fields.solution()
dispIncr.zero()
### NONLINEAR: This moves under SNES control as IntegrateResidual()
Modified: short/3D/PyLith/trunk/pylith/problems/Problem.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Problem.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/problems/Problem.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -17,6 +17,7 @@
## Factory: problem.
from pylith.utils.PetscComponent import PetscComponent
+from pylith.utils.NullComponent import NullComponent
# ITEM FACTORIES ///////////////////////////////////////////////////////
@@ -67,7 +68,6 @@
##
## \b Properties
## @li \b dimension Spatial dimension of problem space.
- ## @li \b useGravity Gravity on (true) or off (false).
##
## \b Facilities
## @li \b normalizer Nondimensionalizer for problem.
@@ -110,10 +110,9 @@
interfaces.meta['tip'] = "Interior surfaces with constraints or " \
"constitutive models."
- from spatialdata.spatialdb.GravityField import GravityField
gravityField = pyre.inventory.facility("gravity_field",
- factory=GravityField,
- family="spatial_database")
+ family="spatial_database",
+ factory=NullComponent)
gravityField.meta['tip'] = "Database used for gravity field."
from pylith.perf.MemoryLogger import MemoryLogger
@@ -231,10 +230,10 @@
self.materials = self.inventory.materials
self.bc = self.inventory.bc
self.interfaces = self.inventory.interfaces
- if self.inventory.useGravity:
- self.gravityField = self.inventory.gravityField
- else:
+ if isinstance(self.inventory.gravityField, NullComponent):
self.gravityField = None
+ else:
+ self.gravityField = self.inventory.gravityField
self.perfLogger = self.inventory.perfLogger
return
Modified: short/3D/PyLith/trunk/pylith/problems/TimeStepUniform.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/TimeStepUniform.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/problems/TimeStepUniform.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -71,7 +71,7 @@
"""
Get number of total time steps (or best guess if adaptive).
"""
- nsteps = int(1 + self.totalTimeN / self.dtN)
+ nsteps = int(1.0 + self.totalTimeN / self.dtN)
return nsteps
Modified: short/3D/PyLith/trunk/pylith/problems/TimeStepUser.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/TimeStepUser.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/problems/TimeStepUser.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -144,6 +144,8 @@
else:
index -= 1
nsteps += 1
+ if nsteps <= 0:
+ raise ValueError("No time steps found.")
return nsteps
Added: short/3D/PyLith/trunk/pylith/utils/NullComponent.py
===================================================================
--- short/3D/PyLith/trunk/pylith/utils/NullComponent.py (rev 0)
+++ short/3D/PyLith/trunk/pylith/utils/NullComponent.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/utils/NullComponent.py
+##
+## @brief Python NullComponent object that is an empty component.
+
+from PetscComponent import PetscComponent
+
+# NullComponent class
+class NullComponent(PetscComponent):
+ """
+ Python NullComponent object that is an empty component.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self):
+ """
+ Constructor.
+ """
+ PetscComponent.__init__(self, name="nullcomponent", facility="nullcomponent")
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _cleanup(self):
+ """
+ Deallocate locally managed data structures.
+ """
+ return
+
+
+# End of file
Modified: short/3D/PyLith/trunk/pylith/utils/__init__.py
===================================================================
--- short/3D/PyLith/trunk/pylith/utils/__init__.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/pylith/utils/__init__.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -12,6 +12,7 @@
__all__ = ['CheckpointTimer',
'CppData',
+ 'NullComponent',
'PetscComponent',
'PetscManager',
'testarray',
Modified: short/3D/PyLith/trunk/tests/1d/line2/matprops.spatialdb
===================================================================
--- short/3D/PyLith/trunk/tests/1d/line2/matprops.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/tests/1d/line2/matprops.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 3
value-names = density vs vp
- value-units = kg/m^3 m/s m/s
+ value-units = kg/m**3 m/s m/s
num-locs = 1
data-dim = 0
space-dim = 1
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/Makefile.am 2009-05-18 23:02:05 UTC (rev 15011)
@@ -152,8 +152,13 @@
data/NeumannDataTet4.hh \
data/NeumannDataHex8.hh
-AM_CPPFLAGS = $(PETSC_SIEVE_FLAGS) $(PETSC_INCLUDE)
+AM_CPPFLAGS = \
+ $(PYTHON_EGG_CPPFLAGS) -I$(PYTHON_INCDIR) \
+ $(PETSC_SIEVE_FLAGS) $(PETSC_INCLUDE)
+testbc_LDFLAGS = \
+ $(AM_LDFLAGS) $(PYTHON_LA_LDFLAGS)
+
testbc_LDADD = \
-lcppunit \
$(top_builddir)/libsrc/libpylith.la \
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -280,9 +280,11 @@
// Set coordinate system
spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
cs.setSpaceDim(mesh->dimension());
cs.initialize();
mesh->coordsys(&cs);
+ mesh->nondimensionalize(normalizer);
// Set up quadrature
_quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryMesh.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestBoundaryMesh.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -51,15 +51,18 @@
CPPUNIT_ASSERT(0 != _data);
topology::Mesh mesh;
+
meshio::MeshIOAscii iohandler;
iohandler.filename(_data->filename);
iohandler.read(&mesh);
- // Set coordinate system
+ // Set up coordinates
spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
cs.setSpaceDim(mesh.dimension());
cs.initialize();
mesh.coordsys(&cs);
+ mesh.nondimensionalize(normalizer);
// Create submesh
topology::SubMesh submesh(mesh, _data->bcLabel);
@@ -116,11 +119,13 @@
iohandler.filename(_data->filename);
iohandler.read(&mesh);
- // Set coordinate system
+ // Set up coordinates
spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
cs.setSpaceDim(mesh.dimension());
cs.initialize();
mesh.coordsys(&cs);
+ mesh.nondimensionalize(normalizer);
// Adjust topology
faults::FaultCohesiveKin fault;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -405,9 +405,11 @@
iohandler.read(mesh);
spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
cs.setSpaceDim(mesh->dimension());
cs.initialize();
mesh->coordsys(&cs);
+ mesh->nondimensionalize(normalizer);
spatialdata::spatialdb::SimpleDB db("TestDirichletBC initial");
spatialdata::spatialdb::SimpleIOAscii dbIO;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -283,9 +283,11 @@
iohandler.read(mesh);
spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
cs.setSpaceDim(mesh->dimension());
cs.initialize();
mesh->coordsys(&cs);
+ mesh->nondimensionalize(normalizer);
// Setup boundary condition A
spatialdata::spatialdb::SimpleDB db("TestDirichletBCMulti initial");
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBoundary.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBoundary.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBoundary.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -23,6 +23,7 @@
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
#include "spatialdata/spatialdb/UniformDB.hh" // USES UniformDB
@@ -116,9 +117,11 @@
iohandler.read(mesh);
spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
cs.setSpaceDim(mesh->dimension());
cs.initialize();
mesh->coordsys(&cs);
+ mesh->nondimensionalize(normalizer);
spatialdata::spatialdb::SimpleDB db("TestDirichletBoundary initial");
spatialdata::spatialdb::SimpleIOAscii dbIO;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -224,9 +224,11 @@
// Set up coordinates
spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
cs.setSpaceDim(mesh->dimension());
cs.initialize();
mesh->coordsys(&cs);
+ mesh->nondimensionalize(normalizer);
// Set up quadrature
_quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticisotropic3d.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticisotropic3d.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticisotropic3d.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 3
value-names = density vs vp
- value-units = kg/m^3 m/s m/s
+ value-units = kg/m**3 m/s m/s
num-locs = 1
data-dim = 0
space-dim = 3
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticplanestrain.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticplanestrain.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticplanestrain.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 3
value-names = density vs vp
- value-units = kg/m^3 m/s m/s
+ value-units = kg/m**3 m/s m/s
num-locs = 1
data-dim = 0
space-dim = 2
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticstrain1d.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticstrain1d.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/elasticstrain1d.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 3
value-names = density vs vp
- value-units = kg/m^3 m/s m/s
+ value-units = kg/m**3 m/s m/s
num-locs = 2
data-dim = 1
space-dim = 1
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/test_bc.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/test_bc.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/test_bc.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -11,6 +11,7 @@
//
#include "petsc.h"
+#include <Python.h>
#include <cppunit/extensions/TestFactoryRegistry.h>
@@ -31,9 +32,14 @@
try {
// Initialize PETSc
- PetscErrorCode err = PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
+ PetscErrorCode err = PetscInitialize(&argc, &argv,
+ PETSC_NULL, PETSC_NULL);
CHKERRQ(err);
+ // Initialize Python (to eliminate need to initialize when
+ // parsing units in spatial databases).
+ Py_Initialize();
+
// Create event manager and test controller
CppUnit::TestResult controller;
@@ -53,9 +59,11 @@
CppUnit::TextOutputter outputter(&result, std::cerr);
outputter.write();
+ // Finalize Python
+ Py_Finalize();
+
// Finalize PETSc
- err = PetscFinalize();
- CHKERRQ(err);
+ err = PetscFinalize(); CHKERRQ(err);
} catch (...) {
abort();
} // catch
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am 2009-05-18 23:02:05 UTC (rev 15011)
@@ -185,8 +185,13 @@
data/CohesiveKinSrcsDataTet4.hh \
data/CohesiveKinSrcsDataHex8.hh
-AM_CPPFLAGS = $(PETSC_SIEVE_FLAGS) $(PETSC_INCLUDE)
+AM_CPPFLAGS = \
+ $(PYTHON_EGG_CPPFLAGS) -I$(PYTHON_INCDIR) \
+ $(PETSC_SIEVE_FLAGS) $(PETSC_INCLUDE)
+testfaults_LDFLAGS = \
+ $(AM_LDFLAGS) $(PYTHON_LA_LDFLAGS)
+
testfaults_LDADD = \
-lcppunit -ldl \
$(top_builddir)/libsrc/libpylith.la \
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -238,9 +238,9 @@
const ALE::Obj<RealSection>& residualSection = residual.section();
CPPUNIT_ASSERT(!residualSection.isNull());
- const ALE::Obj<RealSection>& solutionSection =
- fields.get("solution").section();
- CPPUNIT_ASSERT(!solutionSection.isNull());
+ const ALE::Obj<RealSection>& dispSection =
+ fields.get("disp(t)").section();
+ CPPUNIT_ASSERT(!dispSection.isNull());
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
CPPUNIT_ASSERT(!sieveMesh.isNull());
@@ -253,12 +253,12 @@
for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex)
- solutionSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ dispSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
const double t = 2.134;
const double dt = 0.01;
fault.timeStep(dt);
- { // Integrate residual with solution (as opposed to solution increment).
+ { // Integrate residual with disp (as opposed to disp increment).
fault.useSolnIncr(false);
fault.integrateResidual(residual, t, &fields);
@@ -294,14 +294,14 @@
CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
} // for
} // for
- } // Integrate residual with solution (as opposed to solution increment).
+ } // Integrate residual with disp (as opposed to disp increment).
residual.zero();
- { // Integrate residual with solution increment.
+ { // Integrate residual with disp increment.
fault.useSolnIncr(true);
fault.integrateResidual(residual, t, &fields);
- residual.view("RESIDUAL"); // DEBUGGING
+ //residual.view("RESIDUAL"); // DEBUGGING
// Check values
const double* valsE = _data->valsResidualIncr;
@@ -333,7 +333,7 @@
CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
} // for
} // for
- } // Integrate residual with solution increment.
+ } // Integrate residual with disp increment.
} // testIntegrateResidual
// ----------------------------------------------------------------------
@@ -346,8 +346,8 @@
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- const ALE::Obj<RealSection>& solutionSection = fields.get("solution").section();
- CPPUNIT_ASSERT(!solutionSection.isNull());
+ const ALE::Obj<RealSection>& dispSection = fields.get("disp(t)").section();
+ CPPUNIT_ASSERT(!dispSection.isNull());
const int spaceDim = _data->spaceDim;
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
@@ -360,7 +360,7 @@
for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex) {
- solutionSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ dispSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
} // for
topology::Jacobian jacobian(fields);
@@ -374,7 +374,7 @@
//MatView(jacobian, PETSC_VIEWER_STDOUT_WORLD); // DEBUGGING
const double* valsE = _data->valsJacobian;
- const int nrowsE = solutionSection->sizeWithBC();
+ const int nrowsE = dispSection->sizeWithBC();
const int ncolsE = nrowsE;
int nrows = 0;
@@ -431,8 +431,8 @@
const ALE::Obj<RealSection>& residualSection = residual.section();
CPPUNIT_ASSERT(!residualSection.isNull());
- const ALE::Obj<RealSection>& solutionSection = fields.get("solution").section();
- CPPUNIT_ASSERT(!solutionSection.isNull());
+ const ALE::Obj<RealSection>& dispSection = fields.get("disp(t)").section();
+ CPPUNIT_ASSERT(!dispSection.isNull());
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
CPPUNIT_ASSERT(!sieveMesh.isNull());
@@ -444,12 +444,12 @@
for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex)
- solutionSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ dispSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
const double t = 2.134;
const double dt = 0.01;
fault.timeStep(dt);
- { // Integrate residual with solution (as opposed to solution increment).
+ { // Integrate residual with disp (as opposed to disp increment).
fault.useSolnIncr(false);
fault.integrateResidualAssembled(residual, t, &fields);
@@ -485,10 +485,10 @@
} // for
} // if/else
} // for
- } // Integrate residual with solution (as opposed to solution increment).
+ } // Integrate residual with disp (as opposed to disp increment).
residual.zero();
- { // Integrate residual with solution increment.
+ { // Integrate residual with disp increment.
fault.useSolnIncr(true);
fault.integrateResidualAssembled(residual, t, &fields);
@@ -524,7 +524,7 @@
} // for
} // if/else
} // for
- } // Integrate residual with solution increment.
+ } // Integrate residual with disp increment.
} // testIntegrateResidualAssembled
// ----------------------------------------------------------------------
@@ -538,8 +538,8 @@
_initialize(&mesh, &fault, &fields);
const int spaceDim = _data->spaceDim;
- const ALE::Obj<RealSection>& solutionSection = fields.get("solution").section();
- CPPUNIT_ASSERT(!solutionSection.isNull());
+ const ALE::Obj<RealSection>& dispSection = fields.get("disp(t)").section();
+ CPPUNIT_ASSERT(!dispSection.isNull());
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
CPPUNIT_ASSERT(!sieveMesh.isNull());
@@ -551,7 +551,7 @@
for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex)
- solutionSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ dispSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
topology::Jacobian jacobian(fields);
@@ -564,7 +564,7 @@
//MatView(jacobian, PETSC_VIEWER_STDOUT_WORLD); // DEBUGGING
const double* valsE = _data->valsJacobian;
- const int nrowsE = solutionSection->sizeWithBC();
+ const int nrowsE = dispSection->sizeWithBC();
const int ncolsE = nrowsE;
PetscMat jacobianMat = jacobian.matrix();
@@ -623,10 +623,10 @@
_initialize(&mesh, &fault, &fields);
const int spaceDim = _data->spaceDim;
- const ALE::Obj<RealSection>& solutionSection = fields.get("solution").section();
- CPPUNIT_ASSERT(!solutionSection.isNull());
- { // setup solution
- solutionSection->zero();
+ const ALE::Obj<RealSection>& dispSection = fields.get("disp(t)").section();
+ CPPUNIT_ASSERT(!dispSection.isNull());
+ { // setup disp
+ dispSection->zero();
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
CPPUNIT_ASSERT(!sieveMesh.isNull());
@@ -638,8 +638,8 @@
for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex)
- solutionSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
- } // setup solution
+ dispSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ } // setup disp
topology::Field<topology::Mesh>& residual = fields.get("residual");
const double t = 2.134;
@@ -725,9 +725,9 @@
_initialize(&mesh, &fault, &fields);
const int spaceDim = _data->spaceDim;
- const ALE::Obj<RealSection>& solutionSection = fields.get("solution").section();
- CPPUNIT_ASSERT(!solutionSection.isNull());
- { // setup solution
+ const ALE::Obj<RealSection>& dispSection = fields.get("disp(t)").section();
+ CPPUNIT_ASSERT(!dispSection.isNull());
+ { // setup disp
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
CPPUNIT_ASSERT(!sieveMesh.isNull());
const ALE::Obj<SieveMesh::label_sequence>& vertices =
@@ -739,9 +739,9 @@
for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
++v_iter, ++iVertex) {
- solutionSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ dispSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
} // for
- } // setup solution
+ } // setup disp
CPPUNIT_ASSERT(0 != fault._faultMesh);
topology::Field<topology::SubMesh> tractions(*fault._faultMesh);
@@ -753,7 +753,7 @@
const double t = 0;
fault.updateStateVars(t, &fields);
- fault._calcTractionsChange(&tractions, fields.get("solution"));
+ fault._calcTractionsChange(&tractions, fields.get("disp(t)"));
int iVertex = 0;
const double tolerance = 1.0e-06;
@@ -792,14 +792,14 @@
const double* tractionsVertex = tractionsSection->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != tractionsVertex);
- fiberDim = solutionSection->getFiberDimension(meshVertex);
+ fiberDim = dispSection->getFiberDimension(meshVertex);
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const double* solutionVertex = solutionSection->restrictPoint(meshVertex);
- CPPUNIT_ASSERT(0 != solutionVertex);
+ const double* dispVertex = dispSection->restrictPoint(meshVertex);
+ CPPUNIT_ASSERT(0 != dispVertex);
const double scale = _data->pseudoStiffness / _data->area[iVertex];
for (int iDim=0; iDim < spaceDim; ++iDim) {
- const double tractionE = solutionVertex[iDim] * scale;
+ const double tractionE = dispVertex[iDim] * scale;
if (tractionE > 1.0)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, tractionsVertex[iDim]/tractionE,
tolerance);
@@ -831,9 +831,11 @@
//(*mesh)->setDebug(true); // DEBUGGING
spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
cs.setSpaceDim(mesh->dimension());
cs.initialize();
mesh->coordsys(&cs);
+ mesh->nondimensionalize(normalizer);
_quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
_data->basisDeriv,
@@ -898,8 +900,8 @@
// Setup fields
fields->add("residual", "residual");
- fields->add("solution", "displacement");
- fields->solutionName("solution");
+ fields->add("disp(t)", "displacement");
+ fields->solutionName("disp(t)");
const int spaceDim = _data->spaceDim;
topology::Field<topology::Mesh>& residual = fields->get("residual");
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_1d.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_1d.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_1d.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 2
value-names = density vs
- value-units = none none
+ value-units = kg/m**3 m/s
num-locs = 1
data-dim = 0
space-dim = 1
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_2d.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_2d.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_2d.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 2
value-names = density vs
- value-units = none none
+ value-units = kg/m**3 m/s
num-locs = 1
data-dim = 0
space-dim = 2
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_3d.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_3d.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/bulkprops_3d.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 2
value-names = density vs
- value-units = none none
+ value-units = kg/m**3 m/s
num-locs = 1
data-dim = 0
space-dim = 3
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -11,6 +11,7 @@
//
#include "petsc.h"
+#include <Python.h>
#include <cppunit/extensions/TestFactoryRegistry.h>
@@ -31,9 +32,14 @@
try {
// Initialize PETSc
- PetscErrorCode err = PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
+ PetscErrorCode err = PetscInitialize(&argc, &argv,
+ PETSC_NULL, PETSC_NULL);
CHKERRQ(err);
+ // Initialize Python (to eliminate need to initialize when
+ // parsing units in spatial databases).
+ Py_Initialize();
+
// Create event manager and test controller
CppUnit::TestResult controller;
@@ -53,9 +59,11 @@
CppUnit::TextOutputter outputter(&result, std::cerr);
outputter.write();
+ // Finalize Python
+ Py_Finalize();
+
// Finalize PETSc
- err = PetscFinalize();
- CHKERRQ(err);
+ err = PetscFinalize(); CHKERRQ(err);
} catch (...) {
abort();
} // catch
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2009-05-18 23:02:05 UTC (rev 15011)
@@ -216,8 +216,13 @@
data/QuadratureData3DLinear.hh \
data/QuadratureData3DQuadratic.hh
-AM_CPPFLAGS = $(PETSC_SIEVE_FLAGS) $(PETSC_INCLUDE)
+AM_CPPFLAGS = \
+ $(PYTHON_EGG_CPPFLAGS) -I$(PYTHON_INCDIR) \
+ $(PETSC_SIEVE_FLAGS) $(PETSC_INCLUDE)
+testfeassemble_LDFLAGS = \
+ $(AM_LDFLAGS) $(PYTHON_LA_LDFLAGS)
+
testfeassemble_LDADD = \
-lcppunit -ldl \
$(top_builddir)/libsrc/libpylith.la \
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -292,10 +292,6 @@
CPPUNIT_ASSERT(0 != _material);
// Setup mesh
- spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(_data->spaceDim);
- cs.initialize();
- mesh->coordsys(&cs);
mesh->createSieveMesh(_data->cellDim);
const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
CPPUNIT_ASSERT(!sieveMesh.isNull());
@@ -341,6 +337,13 @@
_data->quadWts, _data->numQuadPts,
_data->spaceDim);
+ spatialdata::units::Nondimensional normalizer;
+ spatialdata::geocoords::CSCart cs;
+ cs.setSpaceDim(_data->spaceDim);
+ cs.initialize();
+ mesh->coordsys(&cs);
+ mesh->nondimensionalize(normalizer);
+
// Setup gravityField
_gravityField = 0;
@@ -350,8 +353,6 @@
spatialdata::spatialdb::SimpleDB dbProperties;
dbProperties.ioHandler(&iohandler);
- spatialdata::units::Nondimensional normalizer;
-
_material->id(_data->matId);
_material->label(_data->matLabel);
_material->dbProperties(&dbProperties);
@@ -370,9 +371,6 @@
fields->add("disp(t)", "displacement");
fields->add("disp(t-dt)", "displacement");
fields->solutionName("dispIncr(t->t+dt)");
- const char* history[] = { "dispIncr(t->t+dt)", "disp(t)", "disp(t-dt)" };
- const int historySize = 3;
- fields->createHistory(history, historySize);
topology::Field<topology::Mesh>& residual = fields->get("residual");
residual.newSection(topology::FieldBase::VERTICES_FIELD, _data->spaceDim);
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticisotropic3d.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticisotropic3d.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticisotropic3d.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 3
value-names = density vs vp
- value-units = kg/m^3 m/s m/s
+ value-units = kg/m**3 m/s m/s
num-locs = 1
data-dim = 0
space-dim = 3
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticplanestrain.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticplanestrain.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticplanestrain.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 3
value-names = density vs vp
- value-units = kg/m^3 m/s m/s
+ value-units = kg/m**3 m/s m/s
num-locs = 1
data-dim = 0
space-dim = 2
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticstrain1d.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticstrain1d.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/elasticstrain1d.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 3
value-names = density vs vp
- value-units = kg/m^3 m/s m/s
+ value-units = kg/m**3 m/s m/s
num-locs = 1
data-dim = 0
space-dim = 1
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/test_feassemble.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/test_feassemble.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/test_feassemble.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -11,6 +11,7 @@
//
#include "petsc.h"
+#include <Python.h>
#include <cppunit/extensions/TestFactoryRegistry.h>
@@ -31,9 +32,14 @@
try {
// Initialize PETSc
- PetscErrorCode err = PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
+ PetscErrorCode err = PetscInitialize(&argc, &argv,
+ PETSC_NULL, PETSC_NULL);
CHKERRQ(err);
+ // Initialize Python (to eliminate need to initialize when
+ // parsing units in spatial databases).
+ Py_Initialize();
+
// Create event manager and test controller
CppUnit::TestResult controller;
@@ -53,9 +59,11 @@
CppUnit::TextOutputter outputter(&result, std::cerr);
outputter.write();
+ // Finalize Python
+ Py_Finalize();
+
// Finalize PETSc
- err = PetscFinalize();
- CHKERRQ(err);
+ err = PetscFinalize(); CHKERRQ(err);
} catch (...) {
abort();
} // catch
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am 2009-05-18 23:02:05 UTC (rev 15011)
@@ -76,8 +76,13 @@
data/MaxwellIsotropic3DTimeDepData.hh \
data/header.hh
-AM_CPPFLAGS = $(PETSC_SIEVE_FLAGS) $(PETSC_INCLUDE)
+AM_CPPFLAGS = \
+ $(PYTHON_EGG_CPPFLAGS) -I$(PYTHON_INCDIR) \
+ $(PETSC_SIEVE_FLAGS) $(PETSC_INCLUDE)
+testmaterials_LDFLAGS = \
+ $(AM_LDFLAGS) $(PYTHON_LA_LDFLAGS)
+
testmaterials_LDADD = \
-lcppunit -ldl \
$(top_builddir)/libsrc/libpylith.la \
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -34,6 +34,8 @@
#include <cstring> // USES memcpy()
+//#define PRECOMPUTE_GEOMETRY
+
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestElasticMaterial );
@@ -672,9 +674,11 @@
// Set up coordinates
spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
cs.setSpaceDim(mesh->dimension());
cs.initialize();
mesh->coordsys(&cs);
+ mesh->nondimensionalize(normalizer);
// Setup quadrature
feassemble::Quadrature<topology::Mesh> quadrature;
@@ -711,7 +715,9 @@
// Compute geometry for cells
quadrature.initializeGeometry();
+#if defined(PRECOMPUTE_GEOMETRY)
quadrature.computeGeometry(*mesh, cells);
+#endif
spatialdata::spatialdb::SimpleDB db;
spatialdata::spatialdb::SimpleIOAscii dbIO;
@@ -731,8 +737,6 @@
dbStrain.ioHandler(&dbIOStrain);
dbStrain.queryType(spatialdata::spatialdb::SimpleDB::NEAREST);
- spatialdata::units::Nondimensional normalizer;
-
material->dbProperties(&db);
material->id(materialId);
material->label("my_material");
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -35,6 +35,8 @@
#include <cassert> // USES assert()
#include <cmath> // USES sqrt()
+//#define PRECOMPUTE_GEOMETRY
+
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestMaterial );
@@ -194,7 +196,9 @@
// Compute geometry for cells
quadrature.initializeGeometry();
+#if defined(PRECOMPUTE_GEOMETRY)
quadrature.computeGeometry(mesh, cells);
+#endif
spatialdata::spatialdb::SimpleDB db;
spatialdata::spatialdb::SimpleIOAscii dbIO;
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/matinitialize.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/matinitialize.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/matinitialize.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 3
value-names = density vs vp
- value-units = kg/m^3 m/s m/s
+ value-units = kg/m**3 m/s m/s
num-locs = 2
data-dim = 1
space-dim = 1
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/test_materials.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/test_materials.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/test_materials.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -11,6 +11,7 @@
//
#include "petsc.h"
+#include <Python.h>
#include <cppunit/extensions/TestFactoryRegistry.h>
@@ -31,9 +32,14 @@
try {
// Initialize PETSc
- PetscErrorCode err = PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
+ PetscErrorCode err = PetscInitialize(&argc, &argv,
+ PETSC_NULL, PETSC_NULL);
CHKERRQ(err);
+ // Initialize Python (to eliminate need to initialize when
+ // parsing units in spatial databases).
+ Py_Initialize();
+
// Create event manager and test controller
CppUnit::TestResult controller;
@@ -53,9 +59,11 @@
CppUnit::TextOutputter outputter(&result, std::cerr);
outputter.write();
+ // Finalize Python
+ Py_Finalize();
+
// Finalize PETSc
- err = PetscFinalize();
- CHKERRQ(err);
+ err = PetscFinalize(); CHKERRQ(err);
} catch (...) {
abort();
} // catch
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -21,11 +21,15 @@
#include "data/MeshData.hh"
+#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
#include <strings.h> // USES strcasecmp()
#include <stdexcept> // USES std::logic_error
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
// ----------------------------------------------------------------------
// Get simple mesh for testing I/O.
@@ -105,6 +109,14 @@
sieveMesh->allocate(groupField);
} // for
+ // Set coordinate system
+ spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
+ cs.setSpaceDim(data.spaceDim);
+ cs.initialize();
+ mesh->coordsys(&cs);
+ mesh->nondimensionalize(normalizer);
+
return mesh;
} // _createMesh
@@ -125,7 +137,7 @@
// Check vertices
const ALE::Obj<SieveMesh::label_sequence>& vertices =
sieveMesh->depthStratum(0);
- const ALE::Obj<SieveMesh::real_section_type>& coordsField =
+ const ALE::Obj<RealSection>& coordsField =
sieveMesh->getRealSection("coordinates");
const int numVertices = vertices->size();
CPPUNIT_ASSERT(!vertices.isNull());
@@ -139,8 +151,7 @@
vertices->begin();
v_iter != vertices->end();
++v_iter) {
- const SieveMesh::real_section_type::value_type *vertexCoords =
- coordsField->restrictPoint(*v_iter);
+ const double* vertexCoords = coordsField->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != vertexCoords);
const double tolerance = 1.0e-06;
for (int iDim=0; iDim < spaceDim; ++iDim)
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -15,7 +15,10 @@
#include "TestMesh.hh" // Implementation of class methods
#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
+
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::topology::TestMesh );
@@ -127,19 +130,91 @@
} // testComm
// ----------------------------------------------------------------------
-// Test initialize().
+// Test nondimensionalize().
void
-pylith::topology::TestMesh::testInitialize(void)
-{ // testInitialize
+pylith::topology::TestMesh::testNondimensionalize(void)
+{ // testNondimensionalizer
+ const double lengthScale = 2.0;
+ const int spaceDim = 2;
+ const int numVertices = 4;
+ const int coordinates[] = {
+ -1.0, 0.0,
+ 0.0, -1.0,
+ 0.0, 1.0,
+ 1.0, 0.0,
+ };
+
Mesh mesh;
+ meshio::MeshIOAscii iohandler;
+ iohandler.filename("data/tri3.mesh");
+ iohandler.read(&mesh);
spatialdata::geocoords::CSCart cs;
cs.setSpaceDim(2);
-
mesh.coordsys(&cs);
- mesh.initialize();
+ spatialdata::units::Nondimensional normalizer;
+ normalizer.lengthScale(lengthScale);
+ mesh.nondimensionalize(normalizer);
-} // testInitialize
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const Mesh::SieveMesh::label_sequence::iterator verticesBegin =
+ vertices->begin();
+ const Mesh::SieveMesh::label_sequence::iterator verticesEnd =
+ vertices->end();
+ CPPUNIT_ASSERT_EQUAL(numVertices, int(vertices->size()));
+ // Check nondimensional coordinates
+ const ALE::Obj<Mesh::RealSection>& coordsField =
+ sieveMesh->getRealSection("coordinates");
+ CPPUNIT_ASSERT(!coordsField.isNull());
+ CPPUNIT_ASSERT_EQUAL(spaceDim,
+ coordsField->getFiberDimension(*verticesBegin));
+ int i = 0;
+ for(Mesh::SieveMesh::label_sequence::iterator v_iter = verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter) {
+ const double* coordsVertex = coordsField->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != coordsVertex);
+ const double tolerance = 1.0e-06;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const double coordE = coordinates[i++] / lengthScale;
+ if (coordE < 1.0)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(coordE, coordsVertex[iDim],
+ tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, coordsVertex[iDim]/coordE,
+ tolerance);
+ } // for
+ } // for
+
+ // Check dimensioned coordinates
+ const ALE::Obj<Mesh::RealSection>& coordsDimField =
+ sieveMesh->getRealSection("coordinates_dimensioned");
+ CPPUNIT_ASSERT(!coordsDimField.isNull());
+ CPPUNIT_ASSERT_EQUAL(spaceDim,
+ coordsDimField->getFiberDimension(*verticesBegin));
+ i = 0;
+ for(Mesh::SieveMesh::label_sequence::iterator v_iter = verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter) {
+ const double* coordsVertex = coordsDimField->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != coordsVertex);
+ const double tolerance = 1.0e-06;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const double coordE = coordinates[i++];
+ if (coordE < 1.0)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(coordE, coordsVertex[iDim],
+ tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, coordsVertex[iDim]/coordE,
+ tolerance);
+ } // for
+ } // for
+} // testNondimensionalize
+
// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.hh 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestMesh.hh 2009-05-18 23:02:05 UTC (rev 15011)
@@ -47,7 +47,7 @@
CPPUNIT_TEST( testDebug );
CPPUNIT_TEST( testDimension );
CPPUNIT_TEST( testComm );
- CPPUNIT_TEST( testInitialize );
+ CPPUNIT_TEST( testNondimensionalize );
CPPUNIT_TEST_SUITE_END();
@@ -75,8 +75,8 @@
/// Test comm().
void testComm(void);
- /// Test initialize().
- void testInitialize(void);
+ /// Test nondimensionalize().
+ void testNondimensionalize(void);
}; // class TestMesh
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestSolutionFields.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestSolutionFields.cc 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestSolutionFields.cc 2009-05-18 23:02:05 UTC (rev 15011)
@@ -85,67 +85,7 @@
} // testSolution
// ----------------------------------------------------------------------
-// Test createHistory().
void
-pylith::topology::TestSolutionFields::testCreateHistory(void)
-{ // testCreateHistory
- Mesh mesh;
- _initialize(&mesh);
- SolutionFields manager(mesh);
-
- const char* labels[] = { "field A", "field B", "field C" };
- const int totalSize = 3;
- const int historySize = 2;
-
- // Add fields
- for (int i=0; i < totalSize; ++i)
- manager.add(labels[i], "displacement");
-
- manager.createHistory(labels, historySize);
- for (int i=0; i < historySize; ++i)
- CPPUNIT_ASSERT_EQUAL(std::string(labels[i]), manager._history[i]);
-} // testCreateHistory
-
-// ----------------------------------------------------------------------
-// Test shiftHistory().
-void
-pylith::topology::TestSolutionFields::testShiftHistory(void)
-{ // testShiftHistory
- Mesh mesh;
- _initialize(&mesh);
- SolutionFields manager(mesh);
-
- const char* fieldNames[] = { "field A", "field B" };
- const int numFields = 2;
- const int fiberDimA = 2;
- const int fiberDimB = 3;
-
- for (int i=0; i < numFields; ++i)
- manager.add(fieldNames[i], "displacement");
- manager.createHistory(fieldNames, numFields);
-
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- Field<Mesh>& fieldA = manager.get(fieldNames[0]);
- Field<Mesh>& fieldB = manager.get(fieldNames[1]);
- fieldA.newSection(vertices, fiberDimA);
- fieldB.newSection(vertices, fiberDimB);
-
- manager.shiftHistory();
- const Field<Mesh>& testA = manager.get(fieldNames[0]);
- const ALE::Obj<Mesh::RealSection>& sectionA = testA.section();
- const Field<Mesh>& testB = manager.get(fieldNames[1]);
- const ALE::Obj<Mesh::RealSection>& sectionB = testB.section();
- CPPUNIT_ASSERT_EQUAL(fiberDimB,
- sectionA->getFiberDimension(*(vertices->begin())));
- CPPUNIT_ASSERT_EQUAL(fiberDimA,
- sectionB->getFiberDimension(*(vertices->begin())));
-
-} // testShiftHistory
-
-// ----------------------------------------------------------------------
-void
pylith::topology::TestSolutionFields::_initialize(Mesh* mesh) const
{ // _initialize
CPPUNIT_ASSERT(0 != mesh);
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestSolutionFields.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestSolutionFields.hh 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestSolutionFields.hh 2009-05-18 23:02:05 UTC (rev 15011)
@@ -42,8 +42,6 @@
CPPUNIT_TEST( testConstructor );
CPPUNIT_TEST( testSolutionName );
CPPUNIT_TEST( testSolution );
- CPPUNIT_TEST( testCreateHistory );
- CPPUNIT_TEST( testShiftHistory );
CPPUNIT_TEST_SUITE_END();
@@ -59,15 +57,6 @@
/// Test solution().
void testSolution(void);
- /// Test createHistory().
- void testCreateHistory(void);
-
- /// Test shiftHistory().
- void testShiftHistory(void);
-
- /// Test getFieldByHistory().
- void testGetFieldByHistory(void);
-
// PRIVATE METHODS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/unittests/pytests/bc/TestAbsorbingDampers.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/bc/TestAbsorbingDampers.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/pytests/bc/TestAbsorbingDampers.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -227,7 +227,6 @@
fields.add("disp(t)", "displacement")
fields.add("disp(t-dt)", "displacement")
fields.solutionName("disp(t+dt)")
- fields.createHistory(["disp(t+dt)", "disp(t)", "disp(t-dt)"])
residual = fields.get("residual")
residual.newSection(residual.VERTICES_FIELD, cs.spaceDim())
Modified: short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -345,9 +345,9 @@
from pylith.topology.SolutionFields import SolutionFields
fields = SolutionFields(mesh)
fields.add("residual", "residual")
- fields.add("solution", "displacement")
- fields.add("disp", "displacement")
- fields.solutionName("solution")
+ fields.add("dispIncr(t->t+dt)", "displacement_increment")
+ fields.add("disp(t)", "displacement")
+ fields.solutionName("dispIncr(t->t+dt)")
residual = fields.get("residual")
residual.newSection(residual.VERTICES_FIELD, cs.spaceDim())
residual.allocate()
Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/data/elasticplanestrain.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/data/elasticplanestrain.spatialdb 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/data/elasticplanestrain.spatialdb 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 3
value-names = density vs vp
- value-units = kg/m^3 m/s m/s
+ value-units = kg/m**3 m/s m/s
num-locs = 1
data-dim = 0
space-dim = 2
Modified: short/3D/PyLith/trunk/unittests/pytests/meshio/data/mesh2Din3D.txt
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/meshio/data/mesh2Din3D.txt 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/pytests/meshio/data/mesh2Din3D.txt 2009-05-18 23:02:05 UTC (rev 15011)
@@ -2,18 +2,18 @@
dimension = 2
use-index-zero = true
vertices = {
- dimension = 2
+ dimension = 3
count = 9
coordinates = {
- 0 -1.000000e+00 3.000000e+00
- 1 2.000000e-01 1.000000e+00
- 2 3.300000e+00 5.000000e-01
- 3 -1.200000e+00 9.000000e-01
- 4 3.000000e-01 9.000000e-01
- 5 1.000000e+00 4.000000e-01
- 6 3.000000e+00 2.900000e+00
- 7 -1.000000e-01 6.000000e+00
- 8 1.200000e+00 -2.000000e-01
+ 0 -1.000000e+00 3.000000e+00 0.000000e+00
+ 1 2.000000e-01 1.000000e+00 0.000000e+00
+ 2 3.300000e+00 5.000000e-01 0.000000e+00
+ 3 -1.200000e+00 9.000000e-01 0.000000e+00
+ 4 3.000000e-01 9.000000e-01 0.000000e+00
+ 5 1.000000e+00 4.000000e-01 0.000000e+00
+ 6 3.000000e+00 2.900000e+00 0.000000e+00
+ 7 -1.000000e-01 6.000000e+00 0.000000e+00
+ 8 1.200000e+00 -2.000000e-01 0.000000e+00
}
}
cells = {
Modified: short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -122,19 +122,6 @@
return
- def test_initialize(self):
- """
- Test initialize().
- """
- from spatialdata.geocoords.CSCart import CSCart
- cs = CSCart()
-
- mesh = Mesh()
- mesh.coordsys(cs)
- mesh.initialize()
- return
-
-
def test_view(self):
"""
Test view().
Modified: short/3D/PyLith/trunk/unittests/pytests/topology/TestSolutionFields.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/topology/TestSolutionFields.py 2009-05-18 15:27:52 UTC (rev 15010)
+++ short/3D/PyLith/trunk/unittests/pytests/topology/TestSolutionFields.py 2009-05-18 23:02:05 UTC (rev 15011)
@@ -82,36 +82,9 @@
return
- def test_createHistory(self):
- """
- Test createHistory().
- """
- fields = self.fields
- fields.add("field A", "A");
- fields.add("field B", "B");
- fields.add("field C", "C");
-
- fields.createHistory(["field B", "field A", "field C"])
- return
-
-
- def test_shiftHistory(self):
- """
- Test shiftHistory().
- """
- fields = self.fields
- fields.add("field A", "A");
- fields.add("field B", "B");
- fields.add("field C", "C");
-
- fields.createHistory(["field B", "field A", "field C"])
- fields.shiftHistory()
- return
-
-
def test_fieldAdd(self):
"""
- Test shiftHistory().
+ Test fieldAdd().
"""
fields = self.fields
fields.add("field A", "A");
More information about the CIG-COMMITS
mailing list