[cig-commits] r16247 - in short/3D/PyLith/trunk: examples/bar_shearwave/tet4 libsrc libsrc/feassemble libsrc/materials modulesrc/feassemble pylith pylith/feassemble pylith/problems unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Tue Feb 9 16:09:42 PST 2010
Author: brad
Date: 2010-02-09 16:09:40 -0800 (Tue, 09 Feb 2010)
New Revision: 16247
Added:
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.hh
short/3D/PyLith/trunk/modulesrc/feassemble/ElasticityExplicitTet4.i
short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicitTet4.py
short/3D/PyLith/trunk/pylith/problems/ExplicitLumpedTet4.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.hh
Modified:
short/3D/PyLith/trunk/examples/bar_shearwave/tet4/pylithapp.cfg
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
short/3D/PyLith/trunk/libsrc/feassemble/feassemblefwd.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.icc
short/3D/PyLith/trunk/modulesrc/feassemble/Makefile.am
short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.i
short/3D/PyLith/trunk/pylith/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
Log:
Created optimized ElasciticyExplicit for tet4.
Modified: short/3D/PyLith/trunk/examples/bar_shearwave/tet4/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/bar_shearwave/tet4/pylithapp.cfg 2010-02-10 00:04:18 UTC (rev 16246)
+++ short/3D/PyLith/trunk/examples/bar_shearwave/tet4/pylithapp.cfg 2010-02-10 00:09:40 UTC (rev 16247)
@@ -42,7 +42,8 @@
# Change to an explicit time stepping formulation
#formulation = pylith.problems.Explicit
-formulation = pylith.problems.ExplicitLumped
+#formulation = pylith.problems.ExplicitLumped
+formulation = pylith.problems.ExplicitLumpedTet4
# Nondimensionalize problem using wave propagation parameters.
normalizer = spatialdata.units.NondimElasticDynamic
@@ -82,7 +83,7 @@
# 3-D simplex cell with 2nd order quadrature
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.shape = tetrahedron
-quadrature.cell.quad_order = 2
+quadrature.cell.quad_order = 1
[pylithapp.timedependent.materials.neg]
@@ -99,7 +100,7 @@
# 3-D simplex cell with 2nd order quadrature
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.shape = tetrahedron
-quadrature.cell.quad_order = 2
+quadrature.cell.quad_order = 1
# ----------------------------------------------------------------------
# boundary conditions
@@ -118,7 +119,7 @@
# 2-D simplex cell in 3-D space with 2nd order quadrature
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.shape = triangle
-quadrature.cell.quad_order = 2
+quadrature.cell.quad_order = 1
[pylithapp.timedependent.bc.x_neg]
# Absorbing boundary condition on +x face of bar
@@ -134,7 +135,7 @@
# 2-D simplex cell in 3-D space with 2nd order quadrature
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.shape = triangle
-quadrature.cell.quad_order = 2
+quadrature.cell.quad_order = 1
[pylithapp.timedependent.bc.y_pos]
# Dirichlet boundary condition on all vertices except fault vertices
@@ -169,7 +170,7 @@
# 2-D simplex cell in 3-D space with 2nd order quadrature
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.shape = triangle
-quadrature.cell.quad_order = 2
+quadrature.cell.quad_order = 1
# Switch to Brune slip time function
eq_srcs.rupture.slip_function = pylith.faults.BruneSlipFn
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2010-02-10 00:04:18 UTC (rev 16246)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2010-02-10 00:09:40 UTC (rev 16247)
@@ -74,6 +74,7 @@
feassemble/IntegratorElasticity.cc \
feassemble/ElasticityImplicit.cc \
feassemble/ElasticityExplicit.cc \
+ feassemble/ElasticityExplicitTet4.cc \
feassemble/IntegratorElasticityLgDeform.cc \
feassemble/ElasticityImplicitLgDeform.cc \
feassemble/ElasticityExplicitLgDeform.cc \
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2010-02-10 00:04:18 UTC (rev 16246)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2010-02-10 00:09:40 UTC (rev 16247)
@@ -889,25 +889,6 @@
const double_array& density = _material->calcDensity();
// Compute Jacobian for inertial terms
-#if 0
- for (int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad]
- / dt2;
- const int iQ = iQuad * numBasis;
- double valJ = 0.0;
- for (int jBasis = 0; jBasis < numBasis; ++jBasis)
- valJ += basis[iQ + jBasis];
- valJ *= wt;
- for (int iBasis = 0; iBasis < numBasis; ++iBasis) {
- // Compute value for DOF 0
- const double valIJ = basis[iQ + iBasis] * valJ;
- _cellVector[iBasis*spaceDim] += valIJ;
- } // for
- } // for
- for (int iBasis = 0; iBasis < numBasis; ++iBasis)
- for (int iDim = 1; iDim < spaceDim; ++iDim)
- _cellVector[iBasis*spaceDim+iDim] += _cellVector[iBasis*spaceDim];
- #else
valuesIJ = 0.0;
for (int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad]
@@ -924,7 +905,6 @@
for (int iBasis = 0; iBasis < numBasis; ++iBasis)
for (int iDim = 0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] += valuesIJ[iBasis];
-#endif
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(4 + numBasis*3) + numBasis*spaceDim);
Added: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc 2010-02-10 00:09:40 UTC (rev 16247)
@@ -0,0 +1,981 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "ElasticityExplicitTet4.hh" // implementation of class methods
+
+#include "Quadrature.hh" // USES Quadrature
+
+#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
+#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+#include "pylith/topology/Jacobian.hh" // USES Jacobian
+
+#include "pylith/utils/array.hh" // USES double_array
+#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
+#include "pylith/utils/lapack.h" // USES LAPACKdgesvd
+
+#include "petscmat.h" // USES PetscMat
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimendional
+
+#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
+#include <cassert> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+
+#define DETAILED_EVENT_LOGGING
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
+const int pylith::feassemble::ElasticityExplicitTet4::_spaceDim = 3;
+const int pylith::feassemble::ElasticityExplicitTet4::_cellDim = 3;
+const int pylith::feassemble::ElasticityExplicitTet4::_tensorSize = 6;
+const int pylith::feassemble::ElasticityExplicitTet4::_numBasis = 4;
+const int pylith::feassemble::ElasticityExplicitTet4::_numQuadPts = 1;
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::ElasticityExplicitTet4::ElasticityExplicitTet4(void) :
+ _dtm1(-1.0)
+{ // constructor
+ _basisDerivArray.resize(_numQuadPts*_numBasis*_spaceDim);
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::ElasticityExplicitTet4::~ElasticityExplicitTet4(void)
+{ // destructor
+ deallocate();
+} // destructor
+
+// ----------------------------------------------------------------------
+// Deallocate PETSc and local data structures.
+void
+pylith::feassemble::ElasticityExplicitTet4::deallocate(void)
+{ // deallocate
+ IntegratorElasticity::deallocate();
+} // deallocate
+
+// ----------------------------------------------------------------------
+// Set time step for advancing from time t to time t+dt.
+void
+pylith::feassemble::ElasticityExplicitTet4::timeStep(const double dt)
+{ // timeStep
+ if (_dt != -1.0)
+ _dtm1 = _dt;
+ else
+ _dtm1 = dt;
+ _dt = dt;
+ assert(_dt == _dtm1); // For now, don't allow variable time step
+ if (0 != _material)
+ _material->timeStep(_dt);
+} // timeStep
+
+// ----------------------------------------------------------------------
+// Set flag for setting constraints for total field solution or
+// incremental field solution.
+void
+pylith::feassemble::ElasticityExplicitTet4::useSolnIncr(const bool flag)
+{ // useSolnIncr
+ if (!flag)
+ throw std::logic_error("Non-incremental solution not supported for "
+ "explicit time integration of elasticity "
+ "equation.");
+} // useSolnIncr
+
+// ----------------------------------------------------------------------
+// Integrate constributions to residual term (r) for operator.
+void
+pylith::feassemble::ElasticityExplicitTet4::integrateResidual(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateResidual
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(0 != _logger);
+ assert(0 != fields);
+
+ const int setupEvent = _logger->eventId("ElIR setup");
+ const int geometryEvent = _logger->eventId("ElIR geometry");
+ const int computeEvent = _logger->eventId("ElIR compute");
+ const int restrictEvent = _logger->eventId("ElIR restrict");
+ const int stateVarsEvent = _logger->eventId("ElIR stateVars");
+ const int stressEvent = _logger->eventId("ElIR stress");
+ const int updateEvent = _logger->eventId("ElIR update");
+
+ _logger->eventBegin(setupEvent);
+
+ // Get cell geometry information that doesn't depend on cell
+ assert(_quadrature->numQuadPts() == _numQuadPts);
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == _numQuadPts);
+ assert(_quadrature->numBasis() == _numBasis);
+ assert(_quadrature->spaceDim() == _spaceDim);
+ assert(_quadrature->cellDim() == _cellDim);
+ assert(_material->tensorSize() == _tensorSize);
+ const int spaceDim = _spaceDim;
+ const int cellDim = _cellDim;
+ const int tensorSize = _tensorSize;
+ const int numBasis = _numBasis;
+ const int numQuadPts = _numQuadPts;
+ /** :TODO:
+ *
+ * If cellDim and spaceDim are different, we need to transform
+ * displacements into cellDim, compute action, and transform result
+ * back into spaceDim. We get this information from the Jacobian and
+ * inverse of the Jacobian.
+ */
+ if (cellDim != spaceDim)
+ throw std::logic_error("Integration for cells with spatial dimensions "
+ "different than the spatial dimension of the "
+ "domain not implemented yet.");
+
+ // Allocate vectors for cell values.
+ double_array dispTCell(numBasis*spaceDim);
+ double_array dispTmdtCell(numBasis*spaceDim);
+ double_array strainCell(numQuadPts*tensorSize);
+ strainCell = 0.0;
+ double_array gravVec(spaceDim);
+ double_array quadPtsGlobal(numQuadPts*spaceDim);
+
+ // Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const int materialId = _material->id();
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
+ assert(!dispTSection.isNull());
+ topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
+ numBasis*spaceDim,
+ &dispTCell[0]);
+ const ALE::Obj<RealSection>& dispTmdtSection =
+ fields->get("disp(t-dt)").section();
+ assert(!dispTmdtSection.isNull());
+ topology::Mesh::RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
+ numBasis*spaceDim,
+ &dispTmdtCell[0]);
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ 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 =
+ _normalizer->pressureScale() / (_normalizer->lengthScale() *
+ _normalizer->densityScale());
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ const double dt2 = dt*dt;
+ assert(dt > 0);
+
+ _logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
+
+ // Loop over cells
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(geometryEvent);
+#endif
+
+ // Compute geometry information for current cell
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const double volume = _volume(coordinatesCell);
+ assert(volume > 0.0);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(geometryEvent);
+ _logger->eventBegin(stateVarsEvent);
+#endif
+
+ // Get state variables for cell.
+ _material->retrievePropsAndVars(*c_iter);
+
+ #if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(stateVarsEvent);
+ _logger->eventBegin(restrictEvent);
+#endif
+
+ // Reset element vector to zero
+ _resetCellVector();
+
+ // Restrict input fields to cell
+ dispTVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+ dispTmdtVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTmdtVisitor);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(restrictEvent);
+ _logger->eventBegin(computeEvent);
+#endif
+
+ const double_array& density = _material->calcDensity();
+ assert(density.size() == 1);
+
+ // Compute body force vector if gravity is being used.
+ if (0 != _gravityField) {
+ const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();
+ assert(0 != cs);
+
+ quadPtsGlobal = 0.0;
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ quadPtsGlobal[iDim] += 0.25 * coordinatesCell[iBasis*spaceDim+iDim];
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
+ lengthScale);
+
+ // Compute action for element body forces
+ const int err = _gravityField->query(&gravVec[0], gravVec.size(),
+ &quadPtsGlobal[0], spaceDim, cs);
+ if (err)
+ throw std::runtime_error("Unable to get gravity vector for point.");
+ _normalizer->nondimensionalize(&gravVec[0], gravVec.size(),
+ gravityScale);
+ const double wtVertex = density[0] * volume / 4.0;
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellVector[iBasis * spaceDim + iDim] += wtVertex * gravVec[iDim];
+ PetscLogFlops(numBasis*spaceDim*2 + numBasis*spaceDim*2);
+ } // if
+
+ // Compute action for inertial terms
+ const double wtVertex = density[0] * volume / (16.0 * dt2);
+ for (int iBasis = 0; iBasis < numBasis; ++iBasis)
+ for (int jBasis = 0; jBasis < numBasis; ++jBasis)
+ for (int iDim = 0; iDim < spaceDim; ++iDim)
+ _cellVector[iBasis * spaceDim + iDim] += wtVertex * (dispTCell[jBasis
+ * spaceDim + iDim] - dispTmdtCell[jBasis * spaceDim + iDim]);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(3 + numBasis*numBasis*spaceDim*3);
+ _logger->eventEnd(computeEvent);
+ _logger->eventBegin(stressEvent);
+#endif
+
+ // Compute B(transpose) * sigma, first computing strains
+ const double x0 = coordinatesCell[0];
+ const double y0 = coordinatesCell[1];
+ const double z0 = coordinatesCell[2];
+
+ const double x1 = coordinatesCell[3];
+ const double y1 = coordinatesCell[4];
+ const double z1 = coordinatesCell[5];
+
+ const double x2 = coordinatesCell[6];
+ const double y2 = coordinatesCell[7];
+ const double z2 = coordinatesCell[8];
+
+ const double x3 = coordinatesCell[9];
+ const double y3 = coordinatesCell[10];
+ const double z3 = coordinatesCell[11];
+
+ const double scaleB = 6.0 * volume;
+ const double b1 = (y1*(z3-z2)-y2*z3+y3*z2-(y3-y2)*z1) / scaleB;
+ const double c1 = (-x1*(z3-z2)+x2*z3-x3*z2-(x2-x3)*z1) / scaleB;
+ const double d1 = (-x2*y3-x1*(y2-y3)+x3*y2+(x2-x3)*y1) / scaleB;
+
+ const double b2 = (-y0*z3-y2*(z0-z3)+(y0-y3)*z2+y3*z0) / scaleB;
+ const double c2 = (x0*z3+x2*(z0-z3)+(x3-x0)*z2-x3*z0) / scaleB;
+ const double d2 = (x2*(y3-y0)-x0*y3-(x3-x0)*y2+x3*y0) / scaleB;
+
+ const double b3 = (-(y1-y0)*z3+y3*(z1-z0)-y0*z1+y1*z0) / scaleB;
+ const double c3 = (-(x0-x1)*z3-x3*(z1-z0)+x0*z1-x1*z0) / scaleB;
+ const double d3 = ((x0-x1)*y3-x0*y1-x3*(y0-y1)+x1*y0) / scaleB;
+
+ const double b4 = (-y0*(z2-z1)+y1*z2-y2*z1+(y2-y1)*z0) / scaleB;
+ const double c4 = (x0*(z2-z1)-x1*z2+x2*z1+(x1-x2)*z0) / scaleB;
+ const double d4 = (x1*y2+x0*(y1-y2)-x2*y1-(x1-x2)*y0) / scaleB;
+
+ assert(strainCell.size() == 6);
+ strainCell[0] = b1 * dispTCell[0] + b2 * dispTCell[3] + b3 * dispTCell[6]
+ + b4 * dispTCell[9];
+ strainCell[1] = c3 * dispTCell[7] + c2 * dispTCell[4] + c4 * dispTCell[10]
+ + c1 * dispTCell[1];
+ strainCell[2] = d3 * dispTCell[8] + d2 * dispTCell[5] + d1 * dispTCell[2]
+ + d4 * dispTCell[11];
+ strainCell[3] = (c4 * dispTCell[9] + b3 * dispTCell[7] + c3 * dispTCell[6]
+ + b2 * dispTCell[4] + c2 * dispTCell[3] + b4 * dispTCell[10] + b1
+ * dispTCell[1] + c1 * dispTCell[0]) / 2.0;
+ strainCell[4] = (c3 * dispTCell[8] + d3 * dispTCell[7] + c2 * dispTCell[5]
+ + d2 * dispTCell[4] + c1 * dispTCell[2] + c4 * dispTCell[11] + d4
+ * dispTCell[10] + d1 * dispTCell[1]) / 2.0;
+ strainCell[5] = (d4 * dispTCell[9] + b3 * dispTCell[8] + d3 * dispTCell[6]
+ + b2 * dispTCell[5] + d2 * dispTCell[3] + b1 * dispTCell[2] + b4
+ * dispTCell[11] + d1 * dispTCell[0]) / 2.0;
+
+ const double_array& stressCell = _material->calcStress(strainCell, true);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(196);
+ _logger->eventEnd(stressEvent);
+ _logger->eventBegin(computeEvent);
+#endif
+
+ assert(_cellVector.size() == 12);
+ assert(stressCell.size() == 6);
+ _cellVector[0] -= (d1*stressCell[5]+c1*stressCell[3]+b1*stressCell[0]) * volume;
+ _cellVector[1] -= (d1*stressCell[4]+b1*stressCell[3]+c1*stressCell[1]) * volume;
+ _cellVector[2] -= (b1*stressCell[5]+c1*stressCell[4]+d1*stressCell[2]) * volume;
+ _cellVector[3] -= (d2*stressCell[5]+c2*stressCell[3]+b2*stressCell[0]) * volume;
+ _cellVector[4] -= (d2*stressCell[4]+b2*stressCell[3]+c2*stressCell[1]) * volume;
+ _cellVector[5] -= (b2*stressCell[5]+c2*stressCell[4]+d2*stressCell[2]) * volume;
+ _cellVector[6] -= (d3*stressCell[5]+c3*stressCell[3]+b3*stressCell[0]) * volume;
+ _cellVector[7] -= (d3*stressCell[4]+b3*stressCell[3]+c3*stressCell[1]) * volume;
+ _cellVector[8] -= (b3*stressCell[5]+c3*stressCell[4]+d3*stressCell[2]) * volume;
+ _cellVector[9] -= (d4*stressCell[5]+c4*stressCell[3]+b4*stressCell[0]) * volume;
+ _cellVector[10] -= (d4*stressCell[4]+b4*stressCell[3]+c4*stressCell[1]) * volume;
+ _cellVector[11] -= (b4*stressCell[5]+c4*stressCell[4]+d4*stressCell[2]) * volume;
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(84);
+ _logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
+
+ // Assemble cell contribution into field
+ residualVisitor.clear();
+ sieveMesh->updateClosure(*c_iter, residualVisitor);
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(updateEvent);
+#endif
+ } // for
+
+#if !defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(cells->size()*(196+84));
+ _logger->eventEnd(computeEvent);
+#endif
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Integrate constributions to residual term (r) for operator.
+void
+pylith::feassemble::ElasticityExplicitTet4::integrateResidualLumped(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateResidualLumped
+ /// Member prototype for _elasticityResidualXD()
+ typedef void (pylith::feassemble::ElasticityExplicitTet4::*elasticityResidual_fn_type)
+ (const double_array&);
+
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(0 != _logger);
+ assert(0 != fields);
+
+ const int setupEvent = _logger->eventId("ElIR setup");
+ const int geometryEvent = _logger->eventId("ElIR geometry");
+ const int computeEvent = _logger->eventId("ElIR compute");
+ const int restrictEvent = _logger->eventId("ElIR restrict");
+ const int stateVarsEvent = _logger->eventId("ElIR stateVars");
+ const int stressEvent = _logger->eventId("ElIR stress");
+ const int updateEvent = _logger->eventId("ElIR update");
+
+ _logger->eventBegin(setupEvent);
+
+ // Get cell geometry information that doesn't depend on cell
+ assert(_quadrature->numQuadPts() == _numQuadPts);
+ assert(_quadrature->numBasis() == _numBasis);
+ assert(_quadrature->spaceDim() == _spaceDim);
+ assert(_quadrature->cellDim() == _cellDim);
+ assert(_material->tensorSize() == _tensorSize);
+ const int spaceDim = _spaceDim;
+ const int cellDim = _cellDim;
+ const int tensorSize = _tensorSize;
+ const int numBasis = _numBasis;
+ const int numQuadPts = _numQuadPts;
+ /** :TODO:
+ *
+ * If cellDim and spaceDim are different, we need to transform
+ * displacements into cellDim, compute action, and transform result
+ * back into spaceDim. We get this information from the Jacobian and
+ * inverse of the Jacobian.
+ */
+ if (cellDim != spaceDim)
+ throw std::logic_error("Integration for cells with spatial dimensions "
+ "different than the spatial dimension of the "
+ "domain not implemented yet.");
+
+ // Allocate vectors for cell values.
+ double_array dispTCell(numBasis*spaceDim);
+ double_array dispTmdtCell(numBasis*spaceDim);
+ double_array strainCell(numQuadPts*tensorSize);
+ strainCell = 0.0;
+ double_array gravVec(spaceDim);
+ double_array quadPtsGlobal(numQuadPts*spaceDim);
+
+ // Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const int materialId = _material->id();
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
+ assert(!dispTSection.isNull());
+ topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
+ numBasis*spaceDim,
+ &dispTCell[0]);
+ const ALE::Obj<RealSection>& dispTmdtSection =
+ fields->get("disp(t-dt)").section();
+ assert(!dispTmdtSection.isNull());
+ topology::Mesh::RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
+ numBasis*spaceDim,
+ &dispTmdtCell[0]);
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ 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 =
+ _normalizer->pressureScale() / (_normalizer->lengthScale() *
+ _normalizer->densityScale());
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ const double dt2 = dt*dt;
+ assert(dt > 0);
+ double_array valuesIJ(numBasis);
+
+ _logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
+
+ // Loop over cells
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(restrictEvent);
+#endif
+
+ // Restrict input fields to cell
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+
+ dispTVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+
+ dispTmdtVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTmdtVisitor);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(restrictEvent);
+ _logger->eventBegin(geometryEvent);
+#endif
+
+ // Compute geometry information for current cell
+ const double volume = _volume(coordinatesCell);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(geometryEvent);
+ _logger->eventBegin(stateVarsEvent);
+#endif
+
+ // Get density at quadrature points for this cell
+ const double_array& density = _material->calcDensity();
+
+ // Get state variables for cell.
+ _material->retrievePropsAndVars(*c_iter);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(stateVarsEvent);
+ _logger->eventBegin(computeEvent);
+#endif
+
+ // Reset element vector to zero
+ _resetCellVector();
+
+ // Compute body force vector if gravity is being used.
+ if (0 != _gravityField) {
+ const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();
+ assert(0 != cs);
+
+ quadPtsGlobal = 0.0;
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ quadPtsGlobal[iDim] += 0.25 * coordinatesCell[iBasis*spaceDim+iDim];
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
+ lengthScale);
+
+ // Compute action for element body forces
+ const int err = _gravityField->query(&gravVec[0], gravVec.size(),
+ &quadPtsGlobal[0], spaceDim, cs);
+ if (err)
+ throw std::runtime_error("Unable to get gravity vector for point.");
+ _normalizer->nondimensionalize(&gravVec[0], gravVec.size(),
+ gravityScale);
+ const double wtVertex = density[0] * volume / 4.0;
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellVector[iBasis * spaceDim + iDim] += wtVertex * gravVec[iDim];
+ PetscLogFlops(numBasis*spaceDim*2 + numBasis*spaceDim*2);
+ } // if
+
+ // Compute action for inertial terms
+ const double wtVertex = density[0] * volume / (4.0 * dt2);
+ _cellVector += wtVertex * (dispTCell - dispTmdtCell);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(2 + numBasis*spaceDim*3);
+ _logger->eventEnd(computeEvent);
+ _logger->eventBegin(stressEvent);
+#endif
+
+ // Compute B(transpose) * sigma, first computing strains
+ const double x0 = coordinatesCell[0];
+ const double y0 = coordinatesCell[1];
+ const double z0 = coordinatesCell[2];
+
+ const double x1 = coordinatesCell[3];
+ const double y1 = coordinatesCell[4];
+ const double z1 = coordinatesCell[5];
+
+ const double x2 = coordinatesCell[6];
+ const double y2 = coordinatesCell[7];
+ const double z2 = coordinatesCell[8];
+
+ const double x3 = coordinatesCell[9];
+ const double y3 = coordinatesCell[10];
+ const double z3 = coordinatesCell[11];
+
+ const double scaleB = 6.0 * volume;
+ const double b1 = (y1*(z3-z2)-y2*z3+y3*z2-(y3-y2)*z1) / scaleB;
+ const double c1 = (-x1*(z3-z2)+x2*z3-x3*z2-(x2-x3)*z1) / scaleB;
+ const double d1 = (-x2*y3-x1*(y2-y3)+x3*y2+(x2-x3)*y1) / scaleB;
+
+ const double b2 = (-y0*z3-y2*(z0-z3)+(y0-y3)*z2+y3*z0) / scaleB;
+ const double c2 = (x0*z3+x2*(z0-z3)+(x3-x0)*z2-x3*z0) / scaleB;
+ const double d2 = (x2*(y3-y0)-x0*y3-(x3-x0)*y2+x3*y0) / scaleB;
+
+ const double b3 = (-(y1-y0)*z3+y3*(z1-z0)-y0*z1+y1*z0) / scaleB;
+ const double c3 = (-(x0-x1)*z3-x3*(z1-z0)+x0*z1-x1*z0) / scaleB;
+ const double d3 = ((x0-x1)*y3-x0*y1-x3*(y0-y1)+x1*y0) / scaleB;
+
+ const double b4 = (-y0*(z2-z1)+y1*z2-y2*z1+(y2-y1)*z0) / scaleB;
+ const double c4 = (x0*(z2-z1)-x1*z2+x2*z1+(x1-x2)*z0) / scaleB;
+ const double d4 = (x1*y2+x0*(y1-y2)-x2*y1-(x1-x2)*y0) / scaleB;
+
+ assert(strainCell.size() == 6);
+ strainCell[0] = b1 * dispTCell[0] + b2 * dispTCell[3] + b3 * dispTCell[6]
+ + b4 * dispTCell[9];
+ strainCell[1] = c3 * dispTCell[7] + c2 * dispTCell[4] + c4 * dispTCell[10]
+ + c1 * dispTCell[1];
+ strainCell[2] = d3 * dispTCell[8] + d2 * dispTCell[5] + d1 * dispTCell[2]
+ + d4 * dispTCell[11];
+ strainCell[3] = (c4 * dispTCell[9] + b3 * dispTCell[7] + c3 * dispTCell[6]
+ + b2 * dispTCell[4] + c2 * dispTCell[3] + b4 * dispTCell[10] + b1
+ * dispTCell[1] + c1 * dispTCell[0]) / 2.0;
+ strainCell[4] = (c3 * dispTCell[8] + d3 * dispTCell[7] + c2 * dispTCell[5]
+ + d2 * dispTCell[4] + c1 * dispTCell[2] + c4 * dispTCell[11] + d4
+ * dispTCell[10] + d1 * dispTCell[1]) / 2.0;
+ strainCell[5] = (d4 * dispTCell[9] + b3 * dispTCell[8] + d3 * dispTCell[6]
+ + b2 * dispTCell[5] + d2 * dispTCell[3] + b1 * dispTCell[2] + b4
+ * dispTCell[11] + d1 * dispTCell[0]) / 2.0;
+
+ const double_array& stressCell = _material->calcStress(strainCell, true);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(196);
+ _logger->eventEnd(stressEvent);
+ _logger->eventBegin(computeEvent);
+#endif
+
+ assert(_cellVector.size() == 12);
+ assert(stressCell.size() == 6);
+ _cellVector[0] -= (d1*stressCell[5]+c1*stressCell[3]+b1*stressCell[0]) * volume;
+ _cellVector[1] -= (d1*stressCell[4]+b1*stressCell[3]+c1*stressCell[1]) * volume;
+ _cellVector[2] -= (b1*stressCell[5]+c1*stressCell[4]+d1*stressCell[2]) * volume;
+ _cellVector[3] -= (d2*stressCell[5]+c2*stressCell[3]+b2*stressCell[0]) * volume;
+ _cellVector[4] -= (d2*stressCell[4]+b2*stressCell[3]+c2*stressCell[1]) * volume;
+ _cellVector[5] -= (b2*stressCell[5]+c2*stressCell[4]+d2*stressCell[2]) * volume;
+ _cellVector[6] -= (d3*stressCell[5]+c3*stressCell[3]+b3*stressCell[0]) * volume;
+ _cellVector[7] -= (d3*stressCell[4]+b3*stressCell[3]+c3*stressCell[1]) * volume;
+ _cellVector[8] -= (b3*stressCell[5]+c3*stressCell[4]+d3*stressCell[2]) * volume;
+ _cellVector[9] -= (d4*stressCell[5]+c4*stressCell[3]+b4*stressCell[0]) * volume;
+ _cellVector[10] -= (d4*stressCell[4]+b4*stressCell[3]+c4*stressCell[1]) * volume;
+ _cellVector[11] -= (b4*stressCell[5]+c4*stressCell[4]+d4*stressCell[2]) * volume;
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(84);
+ _logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
+
+ // Assemble cell contribution into field
+ residualVisitor.clear();
+ sieveMesh->updateClosure(*c_iter, residualVisitor);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(updateEvent);
+#endif
+ } // for
+
+#if !defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(cells->size()*(2 + numBasis*spaceDim*3));
+ _logger->eventEnd(computeEvent);
+#endif
+} // integrateResidualLumped
+
+// ----------------------------------------------------------------------
+// Compute matrix associated with operator.
+void
+pylith::feassemble::ElasticityExplicitTet4::integrateJacobian(
+ topology::Jacobian* jacobian,
+ const double t,
+ topology::SolutionFields* fields)
+{ // integrateJacobian
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(0 != jacobian);
+ assert(0 != fields);
+
+ const int setupEvent = _logger->eventId("ElIJ setup");
+ const int geometryEvent = _logger->eventId("ElIJ geometry");
+ const int computeEvent = _logger->eventId("ElIJ compute");
+ const int restrictEvent = _logger->eventId("ElIJ restrict");
+ const int stateVarsEvent = _logger->eventId("ElIJ stateVars");
+ const int updateEvent = _logger->eventId("ElIJ update");
+
+ _logger->eventBegin(setupEvent);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const int tensorSize = _material->tensorSize();
+ if (cellDim != spaceDim)
+ throw std::logic_error("Don't know how to integrate elasticity " \
+ "contribution to Jacobian matrix for cells with " \
+ "different dimensions than the spatial dimension.");
+
+ // Allocate vectors for cell data.
+ double_array dispTCell(numBasis*spaceDim);
+
+ // Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const int materialId = _material->id();
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<RealSection>& dispTSection =
+ fields->get("disp(t)").section();
+ assert(!dispTSection.isNull());
+
+ // Get sparse matrix
+ const PetscMat jacobianMat = jacobian->matrix();
+ assert(0 != jacobianMat);
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ const double dt2 = dt*dt;
+ assert(dt > 0);
+
+ const ALE::Obj<SieveMesh::order_type>& globalOrder =
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", dispTSection);
+ assert(!globalOrder.isNull());
+ // We would need to request unique points here if we had an interpolated mesh
+ topology::Mesh::IndicesVisitor jacobianVisitor(*dispTSection, *globalOrder,
+ (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);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
+
+ // Loop over cells
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(geometryEvent);
+#endif
+
+ // Compute geometry information for current cell
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const double volume = _volume(coordinatesCell);
+ assert(volume > 0.0);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(geometryEvent);
+ _logger->eventBegin(stateVarsEvent);
+#endif
+
+ // Get state variables for cell.
+ _material->retrievePropsAndVars(*c_iter);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(stateVarsEvent);
+ _logger->eventBegin(computeEvent);
+#endif
+
+ // Reset element matrix to zero
+ _resetCellMatrix();
+
+ // Get material physical properties at quadrature points for this cell
+ const double_array& density = _material->calcDensity();
+ assert(density.size() == 1);
+
+ // Compute Jacobian for inertial terms
+ const double wtVertex = density[0] * volume / (16.0 * dt2);
+ for (int iBasis = 0; iBasis < numBasis; ++iBasis)
+ for (int jBasis = 0; jBasis < numBasis; ++jBasis)
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const int iBlock = (iBasis*spaceDim + iDim) * (numBasis*spaceDim);
+ const int jBlock = (jBasis*spaceDim + iDim);
+ _cellMatrix[iBlock+jBlock] += wtVertex;
+ } // for
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(numQuadPts*(3+numBasis*numBasis*spaceDim*1));
+ _logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
+
+ // Assemble cell contribution into PETSc matrix.
+ jacobianVisitor.clear();
+ PetscErrorCode err = updateOperator(jacobianMat, *sieveMesh->getSieve(),
+ jacobianVisitor, *c_iter,
+ &_cellMatrix[0], ADD_VALUES);
+ CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(updateEvent);
+#endif
+ } // for
+
+#if !defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(cells->size()*(3+numBasis*numBasis*spaceDim*1));
+ _logger->eventEnd(computeEvent);
+#endif
+
+ _needNewJacobian = false;
+ _material->resetNeedNewJacobian();
+} // integrateJacobian
+
+// ----------------------------------------------------------------------
+// Compute matrix associated with operator.
+void
+pylith::feassemble::ElasticityExplicitTet4::integrateJacobian(
+ topology::Field<topology::Mesh>* jacobian,
+ const double t,
+ topology::SolutionFields* fields)
+{ // integrateJacobian
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(0 != jacobian);
+ assert(0 != fields);
+
+ const int setupEvent = _logger->eventId("ElIJ setup");
+ const int geometryEvent = _logger->eventId("ElIJ geometry");
+ const int computeEvent = _logger->eventId("ElIJ compute");
+ const int restrictEvent = _logger->eventId("ElIJ restrict");
+ const int stateVarsEvent = _logger->eventId("ElIJ stateVars");
+ const int updateEvent = _logger->eventId("ElIJ update");
+
+ _logger->eventBegin(setupEvent);
+
+ // Get cell geometry information that doesn't depend on cell
+ assert(_quadrature->numBasis() == _numBasis);
+ assert(_quadrature->spaceDim() == _spaceDim);
+ assert(_quadrature->cellDim() == _cellDim);
+ assert(_material->tensorSize() == _tensorSize);
+ const int spaceDim = _spaceDim;
+ const int cellDim = _cellDim;
+ const int numBasis = _numBasis;
+ if (cellDim != spaceDim)
+ throw std::logic_error("Don't know how to integrate elasticity " \
+ "contribution to Jacobian matrix for cells with " \
+ "different dimensions than the spatial dimension.");
+
+ // Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const int materialId = _material->id();
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ const double dt2 = dt*dt;
+ assert(dt > 0);
+
+ // Get sections
+ const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
+ assert(!jacobianSection.isNull());
+ topology::Mesh::UpdateAddVisitor jacobianVisitor(*jacobianSection,
+ &_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]);
+
+ _logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
+ // Loop over cells
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(geometryEvent);
+#endif
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const double volume = _volume(coordinatesCell);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(geometryEvent);
+ _logger->eventBegin(stateVarsEvent);
+#endif
+
+ // Get state variables for cell.
+ _material->retrievePropsAndVars(*c_iter);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(stateVarsEvent);
+ _logger->eventBegin(computeEvent);
+#endif
+
+ // Compute Jacobian for inertial terms
+ const double_array& density = _material->calcDensity();
+ assert(volume > 0.0);
+ _cellVector = density[0] * volume / (4.0 * dt2);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(3);
+ _logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
+
+ // Assemble cell contribution into lumped matrix.
+ jacobianVisitor.clear();
+ sieveMesh->updateClosure(*c_iter, jacobianVisitor);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(updateEvent);
+#endif
+ } // for
+
+#if !defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(cells->size()*3);
+ _logger->eventEnd(computeEvent);
+#endif
+
+ _needNewJacobian = false;
+ _material->resetNeedNewJacobian();
+} // integrateJacobian
+
+// ----------------------------------------------------------------------
+// Compute volume of tetrahedral cell.
+double
+pylith::feassemble::ElasticityExplicitTet4::_volume(
+ const double_array& coordinatesCell) const
+{ // __volume
+ assert(12 == coordinatesCell.size());
+
+ const double x0 = coordinatesCell[0];
+ const double y0 = coordinatesCell[1];
+ const double z0 = coordinatesCell[2];
+
+ const double x1 = coordinatesCell[3];
+ const double y1 = coordinatesCell[4];
+ const double z1 = coordinatesCell[5];
+
+ const double x2 = coordinatesCell[6];
+ const double y2 = coordinatesCell[7];
+ const double z2 = coordinatesCell[8];
+
+ const double x3 = coordinatesCell[9];
+ const double y3 = coordinatesCell[10];
+ const double z3 = coordinatesCell[11];
+
+ const double det =
+ x1*(y2*z3-y3*z2)-y1*(x2*z3-x3*z2)+(x2*y3-x3*y2)*z1 -
+ x0*((y2*z3-y3*z2)-y1*(z3-z2)+(y3-y2)*z1) +
+ y0*((x2*z3-x3*z2)-x1*(z3-z2)+(x3-x2)*z1) -
+ z0*((x2*y3-x3*y2)-x1*(y3-y2)+(x3-x2)*y1);
+
+ const double volume = det / 6.0;
+ PetscLogFlops(48);
+
+ return volume;
+} // _volume
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.hh 2010-02-10 00:09:40 UTC (rev 16247)
@@ -0,0 +1,180 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/feassemble/ElasticityExplicitTet4.hh
+ *
+ * @brief Explicit time integration of dynamic elasticity equation
+ * using linear tetrahedral finite-elements.
+ */
+
+#if !defined(pylith_feassemble_elasticityexplicittet4_hh)
+#define pylith_feassemble_elasticityexplicittet4_hh
+
+// Include directives ---------------------------------------------------
+#include "IntegratorElasticity.hh" // ISA IntegratorElasticity
+
+// ElasticityExplicitTet4 ---------------------------------------------------
+/**@brief Explicit time integration of the dynamic elasticity equation
+ * using linear tetrahedral finite-elements.
+ *
+ * Note: This object operates on a single finite-element family, which
+ * is defined by the quadrature and a database of material property
+ * parameters.
+ *
+ * Computes contributions to terms A and r in
+ *
+ * A(t+dt) du(t) = b(t+dt, u(t), u(t-dt)) - A(t+dt) u(t),
+ *
+ * r(t+dt) = b(t+dt) - A(t+dt) (u(t) + du(t))
+ *
+ * where A(t) is a sparse matrix or vector, u(t+dt) is the field we
+ * want to compute at time t+dt, b is a vector that depends on the
+ * field at time t and t-dt, and u0 is zero at unknown DOF and set to
+ * the known values at the constrained DOF.
+ *
+ * Contributions from elasticity include the intertial and stiffness
+ * terms, so this object computes the following portions of A and r:
+ *
+ * A = 1/(dt*dt) [M]
+ *
+ * r = (1/(dt*dt) [M])(- {u(t+dt)} + 2/(dt*dt){u(t)} - {u(t-dt)}) - [K]{u(t)}
+ *
+ * Translational inertia.
+ * - Residual action over cell
+ * \f[
+ * \int_{V^e} \rho N^p \sum_q N^q u_i^q \, dV
+ * \f]
+ * - Jacobian action over cell
+ * \f[
+ * \int_{V^e} (\rho N^q N^q)_i \, dV
+ * \f]
+ *
+ * See governing equations section of user manual for more
+ * information.
+*/
+class pylith::feassemble::ElasticityExplicitTet4 : public IntegratorElasticity
+{ // ElasticityExplicitTet4
+ friend class TestElasticityExplicitTet4; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ ElasticityExplicitTet4(void);
+
+ /// Destructor
+ ~ElasticityExplicitTet4(void);
+
+ /// Deallocate PETSc and local data structures.
+ void deallocate(void);
+
+ /** Set time step for advancing from time t to time t+dt.
+ *
+ * @param dt Time step
+ */
+ void timeStep(const double dt);
+
+ /** Set flag for setting constraints for total field solution or
+ * incremental field solution.
+ *
+ * @param flag True if using incremental solution, false otherwise.
+ */
+ void useSolnIncr(const bool flag);
+
+ /** Integrate contributions to residual term (r) for operator.
+ *
+ * @param residual Field containing values for residual
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateResidual(const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields);
+
+ /** Integrate contributions to residual term (r) for operator.
+ *
+ * @param residual Field containing values for residual
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateResidualLumped(const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields);
+
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param jacobian Sparse matrix for Jacobian of system.
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateJacobian(topology::Jacobian* jacobian,
+ const double t,
+ topology::SolutionFields* const fields);
+
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param jacobian Diagonal matrix (as field) for Jacobian of system.
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateJacobian(topology::Field<topology::Mesh>* jacobian,
+ const double t,
+ topology::SolutionFields* const fields);
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /** Compute volume of tetrahedral cell.
+ *
+ * @param coordinatesCell Coordinates of vertices of cell.
+ * @returns Volume of cell.
+ */
+ double _volume(const double_array& coordinatesCell) const;
+
+ /** Compute volume of tetrahedral cell.
+ *
+ * @param coordinatesCell Coordinates of vertices of cell.
+ * @returns Derivatives of basis functions.
+ */
+ const double_array& _basisDeriv(const double_array& coordinatesCell) const;
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ double_array _basisDerivArray; ///< Array of basis derivatives
+
+ double _dtm1; ///< Time step for t-dt1 -> t
+
+ static const int _spaceDim;
+ static const int _cellDim;
+ static const int _tensorSize;
+ static const int _numBasis;
+ static const int _numQuadPts;
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ /// Not implemented.
+ ElasticityExplicitTet4(const ElasticityExplicitTet4&);
+
+ /// Not implemented
+ const ElasticityExplicitTet4& operator=(const ElasticityExplicitTet4&);
+
+}; // ElasticityExplicitTet4
+
+#endif // pylith_feassemble_elasticityexplicittet4_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2010-02-10 00:04:18 UTC (rev 16246)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2010-02-10 00:09:40 UTC (rev 16247)
@@ -19,6 +19,7 @@
Constraint.hh \
Constraint.icc \
ElasticityExplicit.hh \
+ ElasticityExplicitTet4.hh \
ElasticityExplicitLgDeform.hh \
ElasticityImplicit.hh \
ElasticityImplicitLgDeform.hh \
Modified: short/3D/PyLith/trunk/libsrc/feassemble/feassemblefwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/feassemblefwd.hh 2010-02-10 00:04:18 UTC (rev 16246)
+++ short/3D/PyLith/trunk/libsrc/feassemble/feassemblefwd.hh 2010-02-10 00:09:40 UTC (rev 16247)
@@ -56,6 +56,8 @@
class ElasticityImplicit;
class ElasticityExplicit;
+ class ElasticityExplicitTet4;
+
class IntegratorElasticityLgDeform;
class ElasticityImplicitLgDeform;
class ElasticityExplicitLgDeform;
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2010-02-10 00:04:18 UTC (rev 16246)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2010-02-10 00:09:40 UTC (rev 16247)
@@ -131,26 +131,6 @@
} // retrievePropsAndVars
// ----------------------------------------------------------------------
-// Compute density for cell at quadrature points.
-const pylith::double_array&
-pylith::materials::ElasticMaterial::calcDensity(void)
-{ // calcDensity
- const int numQuadPts = _numQuadPts;
- const int numPropsQuadPt = _numPropsQuadPt;
- const int numVarsQuadPt = _numVarsQuadPt;
- assert(_propertiesCell.size() == numQuadPts*numPropsQuadPt);
- assert(_stateVarsCell.size() == numQuadPts*numVarsQuadPt);
- assert(_densityCell.size() == numQuadPts*1);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- _calcDensity(&_densityCell[iQuad],
- &_propertiesCell[iQuad*numPropsQuadPt], numPropsQuadPt,
- &_stateVarsCell[iQuad*numVarsQuadPt], numVarsQuadPt);
-
- return _densityCell;
-} // calcDensity
-
-// ----------------------------------------------------------------------
// Compute stress tensor for cell at quadrature points.
const pylith::double_array&
pylith::materials::ElasticMaterial::calcStress(const double_array& totalStrain,
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.icc 2010-02-10 00:04:18 UTC (rev 16246)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.icc 2010-02-10 00:09:40 UTC (rev 16247)
@@ -14,6 +14,8 @@
#error "ElasticMaterial.icc can only be included from ElasticMaterial.hh"
#endif
+#include <cassert>
+
// Set database for initial stress state.
inline
void
@@ -49,5 +51,26 @@
return _initialFields;
} // initialFields
+// ----------------------------------------------------------------------
+// Compute density for cell at quadrature points.
+inline
+const pylith::double_array&
+pylith::materials::ElasticMaterial::calcDensity(void)
+{ // calcDensity
+ const int numQuadPts = _numQuadPts;
+ const int numPropsQuadPt = _numPropsQuadPt;
+ const int numVarsQuadPt = _numVarsQuadPt;
+ assert(_propertiesCell.size() == numQuadPts*numPropsQuadPt);
+ assert(_stateVarsCell.size() == numQuadPts*numVarsQuadPt);
+ assert(_densityCell.size() == numQuadPts*1);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ _calcDensity(&_densityCell[iQuad],
+ &_propertiesCell[iQuad*numPropsQuadPt], numPropsQuadPt,
+ &_stateVarsCell[iQuad*numVarsQuadPt], numVarsQuadPt);
+
+ return _densityCell;
+} // calcDensity
+
+
// End of file
Added: short/3D/PyLith/trunk/modulesrc/feassemble/ElasticityExplicitTet4.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/ElasticityExplicitTet4.i (rev 0)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/ElasticityExplicitTet4.i 2010-02-10 00:09:40 UTC (rev 16247)
@@ -0,0 +1,88 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file modulesrc/feassemble/ElasticityExplicitTet4.i
+ *
+ * @brief Python interface to C++ ElasticityExplicitTet4 object.
+ */
+
+namespace pylith {
+ namespace feassemble {
+
+ class ElasticityExplicitTet4 : public IntegratorElasticity
+ { // ElasticityExplicitTet4
+
+ // PUBLIC MEMBERS /////////////////////////////////////////////////
+ public :
+
+ /// Constructor
+ ElasticityExplicitTet4(void);
+
+ /// Destructor
+ ~ElasticityExplicitTet4(void);
+
+ /// Deallocate PETSc and local data structures.
+ void deallocate(void);
+
+ /** Set time step for advancing from time t to time t+dt.
+ *
+ * @param dt Time step
+ */
+ void timeStep(const double dt);
+
+ /** Set flag for setting constraints for total field solution or
+ * incremental field solution.
+ *
+ * @param flag True if using incremental solution, false otherwise.
+ */
+ void useSolnIncr(const bool flag);
+
+ /** Integrate contributions to residual term (r) for operator.
+ *
+ * @param residual Field containing values for residual
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateResidual(const pylith::topology::Field<pylith::topology::Mesh>& residual,
+ const double t,
+ pylith::topology::SolutionFields* const fields);
+
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param jacobian Sparse matrix for Jacobian of system.
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateJacobian(pylith::topology::Jacobian* jacobian,
+ const double t,
+ pylith::topology::SolutionFields* const fields);
+
+ /** Integrate contributions to Jacobian matrix (A) associated
+ * with operator that require assembly across cells, vertices,
+ * or processors.
+ *
+ * @param jacobian Diagonal Jacobian matrix as a field.
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateJacobian(pylith::topology::Field<pylith::topology::Mesh>* jacobian,
+ const double t,
+ pylith::topology::SolutionFields* const fields);
+
+ }; // ElasticityExplicitTet4
+
+ } // feassemble
+} // pylith
+
+
+// End of file
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/Makefile.am 2010-02-10 00:04:18 UTC (rev 16246)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/Makefile.am 2010-02-10 00:09:40 UTC (rev 16247)
@@ -38,6 +38,7 @@
IntegratorElasticity.i \
ElasticityImplicit.i \
ElasticityExplicit.i \
+ ElasticityExplicitTet4.i \
IntegratorElasticityLgDeform.i \
ElasticityImplicitLgDeform.i \
ElasticityExplicitLgDeform.i
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.i 2010-02-10 00:04:18 UTC (rev 16246)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.i 2010-02-10 00:09:40 UTC (rev 16247)
@@ -35,6 +35,7 @@
#include "pylith/feassemble/Quadrature.hh"
#include "pylith/feassemble/ElasticityImplicit.hh"
#include "pylith/feassemble/ElasticityExplicit.hh"
+#include "pylith/feassemble/ElasticityExplicitTet4.hh"
#include "pylith/feassemble/ElasticityImplicitLgDeform.hh"
#include "pylith/feassemble/ElasticityExplicitLgDeform.hh"
@@ -81,6 +82,7 @@
%include "IntegratorElasticity.i"
%include "ElasticityImplicit.i"
%include "ElasticityExplicit.i"
+%include "ElasticityExplicitTet4.i"
%include "IntegratorElasticityLgDeform.i"
%include "ElasticityImplicitLgDeform.i"
%include "ElasticityExplicitLgDeform.i"
Modified: short/3D/PyLith/trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/pylith/Makefile.am 2010-02-10 00:04:18 UTC (rev 16246)
+++ short/3D/PyLith/trunk/pylith/Makefile.am 2010-02-10 00:09:40 UTC (rev 16247)
@@ -43,6 +43,7 @@
feassemble/__init__.py \
feassemble/Constraint.py \
feassemble/ElasticityExplicit.py \
+ feassemble/ElasticityExplicitTet4.py \
feassemble/ElasticityExplicitLgDeform.py \
feassemble/ElasticityImplicit.py \
feassemble/ElasticityImplicitLgDeform.py \
@@ -115,6 +116,7 @@
problems/__init__.py \
problems/Explicit.py \
problems/ExplicitLumped.py \
+ problems/ExplicitLumpedTet4.py \
problems/ExplicitLgDeform.py \
problems/Formulation.py \
problems/Implicit.py \
Added: short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicitTet4.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicitTet4.py (rev 0)
+++ short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicitTet4.py 2010-02-10 00:09:40 UTC (rev 16247)
@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/feassemble/ElasticityExplicitTet4.py
+##
+## @brief Python object for explicit time integration of dynamic
+## elasticity equation using finite-elements.
+##
+## Factory: integrator
+
+from IntegratorElasticity import IntegratorElasticity
+from feassemble import ElasticityExplicitTet4 as ModuleElasticityExplicitTet4
+
+# ElasticityExplicitTet4 class
+class ElasticityExplicitTet4(IntegratorElasticity, ModuleElasticityExplicitTet4):
+ """
+ Python object for explicit time integration of dynamic elasticity
+ equation using finite-elements.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="elasticityexplicit"):
+ """
+ Constructor.
+ """
+ IntegratorElasticity.__init__(self, name)
+ ModuleElasticityExplicitTet4.__init__(self)
+ self._loggingPrefix = "ElEx "
+ return
+
+
+ def initialize(self, totalTime, numTimeSteps, normalizer):
+ """
+ Do initialization.
+ """
+ logEvent = "%sinit" % self._loggingPrefix
+ self._eventLogger.eventBegin(logEvent)
+
+ IntegratorElasticity.initialize(self, totalTime, numTimeSteps, normalizer)
+ ModuleElasticityExplicitTet4.initialize(self, self.mesh)
+ self._initializeOutput(totalTime, numTimeSteps, normalizer)
+
+ self._eventLogger.eventEnd(logEvent)
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def integrator():
+ """
+ Factory associated with ElasticityExplicitTet4.
+ """
+ return ElasticityExplicitTet4()
+
+
+# End of file
Added: short/3D/PyLith/trunk/pylith/problems/ExplicitLumpedTet4.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/ExplicitLumpedTet4.py (rev 0)
+++ short/3D/PyLith/trunk/pylith/problems/ExplicitLumpedTet4.py 2010-02-10 00:09:40 UTC (rev 16247)
@@ -0,0 +1,81 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/problems/ExplicitLumpedTet4.py
+##
+## @brief Python ExplicitLumpedTet4 object for solving equations using an
+## explicit formulation with a lumped Jacobian matrix that is stored
+## as a Field.
+##
+## Factory: pde_formulation
+
+from ExplicitLumped import ExplicitLumped
+from pylith.utils.profiling import resourceUsageString
+
+# ExplicitLumpedTet4 class
+class ExplicitLumpedTet4(ExplicitLumped):
+ """
+ Python ExplicitLumpedTet4 object for solving equations using an explicit
+ formulation.
+
+ The formulation has the general form, [A(t)] {u(t+dt)} = {b(t)},
+ where we want to solve for {u(t+dt)}, A(t) is usually constant
+ (i.e., independent of time), and {b(t)} usually depends on {u(t)}
+ and {u(t-dt)}.
+
+ Jacobian: A(t)
+ solution: u(t+dt)
+ residual: b(t) - A(t) \hat u(t+dt)
+ constant: b(t)
+
+ Factory: pde_formulation.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="explicitlumpedtet4"):
+ """
+ Constructor.
+ """
+ ExplicitLumped.__init__(self, name)
+ self._loggingPrefix = "TSEx "
+ return
+
+
+ def elasticityIntegrator(self):
+ """
+ Get integrator for elastic material.
+ """
+ from pylith.feassemble.ElasticityExplicitTet4 import ElasticityExplicitTet4
+ return ElasticityExplicitTet4()
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based using inventory.
+ """
+ ExplicitLumped._configure(self)
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def pde_formulation():
+ """
+ Factory associated with ExplicitLumped.
+ """
+ return ExplicitLumpedTet4()
+
+
+# End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2010-02-10 00:04:18 UTC (rev 16246)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2010-02-10 00:09:40 UTC (rev 16247)
@@ -59,6 +59,7 @@
TestElasticityExplicitGrav2DQuadratic.cc \
TestElasticityExplicitGrav3DLinear.cc \
TestElasticityExplicitGrav3DQuadratic.cc \
+ TestElasticityExplicitTet4.cc \
TestElasticityImplicit.cc \
TestElasticityImplicit1DLinear.cc \
TestElasticityImplicit1DQuadratic.cc \
@@ -118,6 +119,7 @@
TestElasticityExplicitGrav2DQuadratic.hh \
TestElasticityExplicitGrav3DLinear.hh \
TestElasticityExplicitGrav3DQuadratic.hh \
+ TestElasticityExplicitTet4.hh \
TestElasticityImplicit.hh \
TestElasticityImplicit1DLinear.hh \
TestElasticityImplicit1DQuadratic.hh \
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc 2010-02-10 00:09:40 UTC (rev 16247)
@@ -0,0 +1,476 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestElasticityExplicitTet4.hh" // Implementation of class methods
+
+#include "pylith/feassemble/ElasticityExplicitTet4.hh" // USES ElasticityExplicitTet4
+#include "data/ElasticityExplicitData3DLinear.hh"
+#include "pylith/feassemble/GeometryTet3D.hh" // USES GeometryTet3D
+
+#include "pylith/materials/ElasticIsotropic3D.hh" // USES ElasticIsotropic3D
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+#include "pylith/topology/Jacobian.hh" // USES Jacobian
+
+#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
+#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+#include "spatialdata/spatialdb/GravityField.hh" // USES GravityField
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
+#include <math.h> // USES fabs()
+
+#include <stdexcept> // USES std::exception
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestElasticityExplicitTet4 );
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::feassemble::TestElasticityExplicitTet4::setUp(void)
+{ // setUp
+ _quadrature = new Quadrature<topology::Mesh>();
+ CPPUNIT_ASSERT(0 != _quadrature);
+ GeometryTet3D geometry;
+ _quadrature->refGeometry(&geometry);
+
+ _data = new ElasticityExplicitData3DLinear();
+ CPPUNIT_ASSERT(0 != _data);
+ _material = new materials::ElasticIsotropic3D;
+ CPPUNIT_ASSERT(0 != _material);
+ CPPUNIT_ASSERT_EQUAL(std::string("ElasticIsotropic3D"),
+ std::string(_data->matType));
+ _gravityField = 0;
+} // setUp
+
+// ----------------------------------------------------------------------
+// Tear down testing data.
+void
+pylith::feassemble::TestElasticityExplicitTet4::tearDown(void)
+{ // tearDown
+ delete _data; _data = 0;
+ delete _quadrature; _quadrature = 0;
+ delete _material; _material = 0;
+ delete _gravityField; _gravityField = 0;
+} // tearDown
+
+// ----------------------------------------------------------------------
+// Test constructor.
+void
+pylith::feassemble::TestElasticityExplicitTet4::testConstructor(void)
+{ // testConstructor
+ ElasticityExplicitTet4 integrator;
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test timeStep().
+void
+pylith::feassemble::TestElasticityExplicitTet4::testTimeStep(void)
+{ // testTimeStep
+ ElasticityExplicitTet4 integrator;
+
+ const double dt1 = 2.0;
+ integrator.timeStep(dt1);
+ CPPUNIT_ASSERT_EQUAL(dt1, integrator._dt);
+ integrator.timeStep(dt1);
+ CPPUNIT_ASSERT_EQUAL(dt1, integrator._dtm1);
+ CPPUNIT_ASSERT_EQUAL(dt1, integrator._dt);
+} // testTimeStep
+
+// ----------------------------------------------------------------------
+// Test material().
+void
+pylith::feassemble::TestElasticityExplicitTet4::testMaterial(void)
+{ // testMaterial
+ ElasticityExplicitTet4 integrator;
+
+ materials::ElasticIsotropic3D material;
+ const int id = 3;
+ const std::string label("my material");
+ material.id(id);
+ material.label(label.c_str());
+ integrator.material(&material);
+ CPPUNIT_ASSERT_EQUAL(id, integrator._material->id());
+ CPPUNIT_ASSERT_EQUAL(label, std::string(integrator._material->label()));
+ CPPUNIT_ASSERT_EQUAL(integrator._dt, integrator._material->timeStep());
+ const double dt = 2.0;
+ integrator.timeStep(dt);
+ CPPUNIT_ASSERT_EQUAL(dt, integrator._material->timeStep());
+} // testMaterial
+
+// ----------------------------------------------------------------------
+// Test needNewJacobian().
+void
+pylith::feassemble::TestElasticityExplicitTet4::testNeedNewJacobian(void)
+{ // testNeedNewJacobian
+ ElasticityExplicitTet4 integrator;
+
+ materials::ElasticIsotropic3D material;
+ integrator.material(&material);
+ CPPUNIT_ASSERT_EQUAL(true, integrator.needNewJacobian());
+ integrator._needNewJacobian = false;
+ CPPUNIT_ASSERT_EQUAL(false, integrator.needNewJacobian());
+} // testNeedNewJacobian
+
+// ----------------------------------------------------------------------
+// Test useSolnIncr().
+void
+pylith::feassemble::TestElasticityExplicitTet4::testUseSolnIncr(void)
+{ // testUseSolnIncr
+ ElasticityExplicitTet4 integrator;
+
+ materials::ElasticIsotropic3D material;
+ integrator.material(&material);
+ CPPUNIT_ASSERT_EQUAL(false, integrator._useSolnIncr);
+ try {
+ integrator.useSolnIncr(false);
+
+ // Should have thrown exception, so don't make it here.
+ CPPUNIT_ASSERT(false);
+ } catch (const std::logic_error& err) {
+ // Expect logic error so don't do anything.
+ } catch (...) {
+ CPPUNIT_ASSERT(false);
+ } // try/catch
+} // testUseSolnIncr
+
+// ----------------------------------------------------------------------
+// Test initialize().
+void
+pylith::feassemble::TestElasticityExplicitTet4::testInitialize(void)
+{ // testInitialize
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ ElasticityExplicitTet4 integrator;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+
+} // testInitialize
+
+// ----------------------------------------------------------------------
+// Test integrateResidual().
+void
+pylith::feassemble::TestElasticityExplicitTet4::testIntegrateResidual(void)
+{ // testIntegrateResidual
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ ElasticityExplicitTet4 integrator;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+
+ topology::Field<topology::Mesh>& residual = fields.get("residual");
+ const double t = 1.0;
+ integrator.integrateResidual(residual, t, &fields);
+
+ const double* valsE = _data->valsResidual;
+ const int sizeE = _data->spaceDim * _data->numVertices;
+
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ CPPUNIT_ASSERT(!residualSection.isNull());
+ const double* vals = residualSection->restrictSpace();
+ const int size = residualSection->sizeWithBC();
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+
+#if 0 // DEBUGGING
+ residual.view("RESIDUAL");
+ std::cout << "EXPECTED RESIDUAL" << std::endl;
+ for (int i=0; i < size; ++i)
+ std::cout << " " << valsE[i] << std::endl;
+#endif // DEBUGGING
+
+ const double tolerance = 1.0e-06;
+ for (int i=0; i < size; ++i)
+ if (fabs(valsE[i]) > 1.0)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+} // testIntegrateResidual
+
+// ----------------------------------------------------------------------
+// Test integrateJacobian().
+void
+pylith::feassemble::TestElasticityExplicitTet4::testIntegrateJacobian(void)
+{ // testIntegrateJacobian
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ ElasticityExplicitTet4 integrator;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+ integrator._needNewJacobian = true;
+
+ topology::Jacobian jacobian(fields);
+
+ const double t = 1.0;
+ integrator.integrateJacobian(&jacobian, t, &fields);
+ CPPUNIT_ASSERT_EQUAL(false, integrator.needNewJacobian());
+ jacobian.assemble("final_assembly");
+
+ const double* valsE = _data->valsJacobian;
+ const int nrowsE = _data->numVertices * _data->spaceDim;
+ const int ncolsE = _data->numVertices * _data->spaceDim;
+
+ const PetscMat jacobianMat = jacobian.matrix();
+
+ int nrows = 0;
+ int ncols = 0;
+ MatGetSize(jacobianMat, &nrows, &ncols);
+ CPPUNIT_ASSERT_EQUAL(nrowsE, nrows);
+ CPPUNIT_ASSERT_EQUAL(ncolsE, ncols);
+
+ PetscMat jDense;
+ PetscMat jSparseAIJ;
+ MatConvert(jacobianMat, MATSEQAIJ, MAT_INITIAL_MATRIX, &jSparseAIJ);
+ MatConvert(jSparseAIJ, MATSEQDENSE, MAT_INITIAL_MATRIX, &jDense);
+
+ double_array vals(nrows*ncols);
+ int_array rows(nrows);
+ int_array cols(ncols);
+ for (int iRow=0; iRow < nrows; ++iRow)
+ rows[iRow] = iRow;
+ for (int iCol=0; iCol < ncols; ++iCol)
+ cols[iCol] = iCol;
+ MatGetValues(jDense, nrows, &rows[0], ncols, &cols[0], &vals[0]);
+ const double tolerance = 1.0e-06;
+ for (int iRow=0; iRow < nrows; ++iRow)
+ for (int iCol=0; iCol < ncols; ++iCol) {
+ const int index = ncols*iRow+iCol;
+ if (fabs(valsE[index]) > 1.0)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[index]/valsE[index], tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[index], vals[index], tolerance);
+ } // for
+ MatDestroy(jDense);
+ MatDestroy(jSparseAIJ);
+} // testIntegrateJacobian
+
+// ----------------------------------------------------------------------
+// Test integrateJacobian().
+void
+pylith::feassemble::TestElasticityExplicitTet4::testIntegrateJacobianLumped(void)
+{ // testIntegrateJacobian
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ ElasticityExplicitTet4 integrator;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+ integrator._needNewJacobian = true;
+
+ topology::Field<topology::Mesh> jacobian(mesh);
+ jacobian.label("Jacobian");
+ jacobian.vectorFieldType(topology::FieldBase::VECTOR);
+ jacobian.newSection(topology::FieldBase::VERTICES_FIELD, _data->spaceDim);
+ jacobian.allocate();
+
+ const double t = 1.0;
+ integrator.integrateJacobian(&jacobian, t, &fields);
+ CPPUNIT_ASSERT_EQUAL(false, integrator.needNewJacobian());
+ jacobian.complete();
+
+ const double* valsMatrixE = _data->valsJacobian;
+ const int sizeE = _data->numVertices * _data->spaceDim;
+ double_array valsE(sizeE);
+ const int spaceDim = _data->spaceDim;
+ const int numBasis = _data->numVertices;
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const int indexRow = (iBasis*spaceDim+iDim)*numBasis*spaceDim;
+ double value = 0.0;
+ for (int jBasis=0; jBasis < numBasis; ++jBasis)
+ value += valsMatrixE[indexRow + jBasis*spaceDim+iDim];
+ valsE[iBasis*spaceDim+iDim] = value;
+ } // for
+
+#if 0 // DEBUGGING
+ jacobian.view("JACOBIAN");
+ std::cout << "\n\nJACOBIAN FULL" << std::endl;
+ const int n = numBasis*spaceDim;
+ for (int i=0; i < n; ++i)
+ std::cout << " " << valsE[i] << "\n";
+#endif // DEBUGGING
+
+ const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
+ CPPUNIT_ASSERT(!jacobianSection.isNull());
+ const double* vals = jacobianSection->restrictSpace();
+ const int size = jacobianSection->sizeWithBC();
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+
+ const double tolerance = 1.0e-06;
+ for (int i=0; i < size; ++i)
+ if (fabs(valsE[i]) > 1.0)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+} // testIntegrateJacobianLumped
+
+// ----------------------------------------------------------------------
+// Test updateStateVars().
+void
+pylith::feassemble::TestElasticityExplicitTet4::testUpdateStateVars(void)
+{ // testUpdateStateVars
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ ElasticityExplicitTet4 integrator;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+
+ const double t = 1.0;
+ integrator.updateStateVars(t, &fields);
+} // testUpdateStateVars
+
+// ----------------------------------------------------------------------
+// Test StableTimeStep().
+void
+pylith::feassemble::TestElasticityExplicitTet4::testStableTimeStep(void)
+{ // testStableTimeStep
+ topology::Mesh mesh;
+ ElasticityExplicitTet4 integrator;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+
+ const double stableTimeStep = integrator.stableTimeStep(mesh);
+ CPPUNIT_ASSERT_EQUAL(1.0e+30, stableTimeStep);
+} // testStableTimeStep
+
+// ----------------------------------------------------------------------
+// Initialize elasticity integrator.
+void
+pylith::feassemble::TestElasticityExplicitTet4::_initialize(
+ topology::Mesh* mesh,
+ ElasticityExplicitTet4* const integrator,
+ topology::SolutionFields* fields)
+{ // _initialize
+ CPPUNIT_ASSERT(0 != mesh);
+ CPPUNIT_ASSERT(0 != integrator);
+ CPPUNIT_ASSERT(0 != _data);
+ CPPUNIT_ASSERT(0 != _quadrature);
+ CPPUNIT_ASSERT(0 != _material);
+
+ // Setup mesh
+ mesh->createSieveMesh(_data->cellDim);
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ ALE::Obj<SieveMesh::sieve_type> sieve =
+ new SieveMesh::sieve_type(mesh->comm());
+ CPPUNIT_ASSERT(!sieve.isNull());
+
+ // Cells and vertices
+ const bool interpolate = false;
+ ALE::Obj<ALE::Mesh::sieve_type> s =
+ new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
+
+ ALE::SieveBuilder<ALE::Mesh>::buildTopology(s,
+ _data->cellDim, _data->numCells,
+ const_cast<int*>(_data->cells),
+ _data->numVertices,
+ interpolate, _data->numBasis);
+ std::map<ALE::Mesh::point_type,ALE::Mesh::point_type> renumbering;
+ ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
+ sieveMesh->setSieve(sieve);
+ sieveMesh->stratify();
+ ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, _data->spaceDim,
+ _data->vertices);
+
+ // Material ids
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->heightStratum(0);
+ CPPUNIT_ASSERT(!cells.isNull());
+ const ALE::Obj<SieveMesh::label_type>& labelMaterials =
+ sieveMesh->createLabel("material-id");
+ CPPUNIT_ASSERT(!labelMaterials.isNull());
+ int i = 0;
+ for(SieveMesh::label_sequence::iterator e_iter=cells->begin();
+ e_iter != cells->end();
+ ++e_iter)
+ sieveMesh->setValue(labelMaterials, *e_iter, _data->matId);
+
+ // Setup quadrature
+ _quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
+ _data->basisDerivRef, _data->numQuadPts,
+ _data->numBasis, _data->cellDim,
+ _data->quadPts, _data->numQuadPts, _data->cellDim,
+ _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 material
+ spatialdata::spatialdb::SimpleIOAscii iohandler;
+ iohandler.filename(_data->matDBFilename);
+ spatialdata::spatialdb::SimpleDB dbProperties;
+ dbProperties.ioHandler(&iohandler);
+
+ _material->id(_data->matId);
+ _material->label(_data->matLabel);
+ _material->dbProperties(&dbProperties);
+ _material->normalizer(normalizer);
+
+ integrator->quadrature(_quadrature);
+ integrator->gravityField(_gravityField);
+ integrator->timeStep(_data->dt);
+ integrator->material(_material);
+ integrator->initialize(*mesh);
+
+ // Setup fields
+ CPPUNIT_ASSERT(0 != fields);
+ fields->add("residual", "residual");
+ fields->add("dispIncr(t->t+dt)", "displacement_increment");
+ fields->add("disp(t)", "displacement");
+ fields->add("disp(t-dt)", "displacement");
+ fields->solutionName("dispIncr(t->t+dt)");
+
+ topology::Field<topology::Mesh>& residual = fields->get("residual");
+ residual.newSection(topology::FieldBase::VERTICES_FIELD, _data->spaceDim);
+ residual.allocate();
+ residual.zero();
+ fields->copyLayout("residual");
+
+ const int fieldSize = _data->spaceDim * _data->numVertices;
+ topology::Field<topology::Mesh>& dispIncr = fields->get("dispIncr(t->t+dt)");
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ topology::Field<topology::Mesh>& dispTmdt = fields->get("disp(t-dt)");
+ const ALE::Obj<RealSection>& dispIncrSection = dispIncr.section();
+ const ALE::Obj<RealSection>& dispTSection = dispT.section();
+ const ALE::Obj<RealSection>& dispTmdtSection = dispTmdt.section();
+ CPPUNIT_ASSERT(!dispIncrSection.isNull());
+ CPPUNIT_ASSERT(!dispTSection.isNull());
+ CPPUNIT_ASSERT(!dispTmdtSection.isNull());
+ const int offset = _data->numCells;
+ for (int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
+ dispIncrSection->updatePoint(iVertex+offset,
+ &_data->fieldTIncr[iVertex*_data->spaceDim]);
+ dispTSection->updatePoint(iVertex+offset,
+ &_data->fieldT[iVertex*_data->spaceDim]);
+ dispTmdtSection->updatePoint(iVertex+offset,
+ &_data->fieldTmdt[iVertex*_data->spaceDim]);
+ } // for
+} // _initialize
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.hh 2010-02-10 00:09:40 UTC (rev 16247)
@@ -0,0 +1,129 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/feassemble/TestElasticityExplicitTet4.hh
+ *
+ * @brief C++ TestElasticityExplicitTet4 object
+ *
+ * C++ unit testing for ElasticityExplicitTet4.
+ */
+
+#if !defined(pylith_feassemble_testelasticityexplicittet4_hh)
+#define pylith_feassemble_testelasticityexplicittet4_hh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "pylith/feassemble/feassemblefwd.hh" // forward declarations
+#include "pylith/topology/topologyfwd.hh" // USES Mesh, SolutionFields
+#include "pylith/materials/materialsfwd.hh" // USES ElasticMaterial
+
+#include "spatialdata/spatialdb/spatialdbfwd.hh" // USES GravityField
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace feassemble {
+ class TestElasticityExplicitTet4;
+ class IntegratorData;
+ } // feassemble
+} // pylith
+
+/// C++ unit testing for ElasticityExplicitTet4
+class pylith::feassemble::TestElasticityExplicitTet4 : public CppUnit::TestFixture
+{ // class TestElasticityExplicitTet4
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestElasticityExplicitTet4 );
+
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testTimeStep );
+ CPPUNIT_TEST( testMaterial );
+ CPPUNIT_TEST( testNeedNewJacobian );
+ CPPUNIT_TEST( testUseSolnIncr );
+ CPPUNIT_TEST( testInitialize );
+ CPPUNIT_TEST( testIntegrateResidual );
+ CPPUNIT_TEST( testIntegrateJacobian );
+ CPPUNIT_TEST( testIntegrateJacobianLumped );
+ CPPUNIT_TEST( testUpdateStateVars );
+ CPPUNIT_TEST( testStableTimeStep );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Setup testing data.
+ void setUp(void);
+
+ /// Tear down testing data.
+ void tearDown(void);
+
+ /// Test constructor.
+ void testConstructor(void);
+
+ /// Test timeStep().
+ void testTimeStep(void);
+
+ /// Test material().
+ void testMaterial(void);
+
+ /// Test needNewJacobian().
+ void testNeedNewJacobian(void);
+
+ /// Test useSolnIncr().
+ void testUseSolnIncr(void);
+
+ /// Test initialize().
+ void testInitialize(void);
+
+ /// Test integrateResidual().
+ void testIntegrateResidual(void);
+
+ /// Test integrateJacobian().
+ void testIntegrateJacobian(void);
+
+ /// Test integrateJacobianLumped().
+ void testIntegrateJacobianLumped(void);
+
+ /// Test updateStateVars().
+ void testUpdateStateVars(void);
+
+ /// Test StableTimeStep().
+ void testStableTimeStep(void);
+
+ // PROTECTED MEMBERS //////////////////////////////////////////////////
+protected :
+
+ IntegratorData* _data; ///< Data for testing.
+ materials::ElasticMaterial* _material; ///< Elastic material.
+ Quadrature<topology::Mesh>* _quadrature; ///< Quadrature information.
+ spatialdata::spatialdb::GravityField* _gravityField; ///< Gravity field.
+
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /** Initialize elasticity integrator.
+ *
+ * @param mesh Finite-element mesh to initialize.
+ * @param integrator ElasticityIntegrator to initialize.
+ * @param fields Solution fields.
+ */
+ void _initialize(topology::Mesh* mesh,
+ ElasticityExplicitTet4* const integrator,
+ topology::SolutionFields* const fields);
+
+}; // class TestElasticityExplicitTet4
+
+#endif // pylith_feassemble_testelasticityexplicittet4_hh
+
+
+// End of file
More information about the CIG-COMMITS
mailing list