[cig-commits] r16826 - in short/3D/PyLith/trunk: examples/twocells/twoquad4 libsrc libsrc/feassemble modulesrc/feassemble pylith pylith/feassemble pylith/problems unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Fri May 28 19:23:03 PDT 2010
Author: brad
Date: 2010-05-28 19:23:03 -0700 (Fri, 28 May 2010)
New Revision: 16826
Added:
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.hh
short/3D/PyLith/trunk/libsrc/feassemble/tet4_elasticity.wxm
short/3D/PyLith/trunk/libsrc/feassemble/tri3_elasticity.wxm
short/3D/PyLith/trunk/modulesrc/feassemble/ElasticityExplicitTri3.i
short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicitTri3.py
short/3D/PyLith/trunk/pylith/problems/ExplicitLumpedTri3.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.hh
Modified:
short/3D/PyLith/trunk/examples/twocells/twoquad4/pylithapp.cfg
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.hh
short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
short/3D/PyLith/trunk/libsrc/feassemble/feassemblefwd.hh
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/pylith/feassemble/ElasticityExplicitTet4.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc
Log:
Added optimized elasticity integrator for tri3 cells. Still needs some debugging.
Modified: short/3D/PyLith/trunk/examples/twocells/twoquad4/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/twocells/twoquad4/pylithapp.cfg 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/examples/twocells/twoquad4/pylithapp.cfg 2010-05-29 02:23:03 UTC (rev 16826)
@@ -96,23 +96,3 @@
# start_in_debugger = true
# debugger_timeout = 100
-
-# TESTING OF FAULT PRECONDITIONER
-# Field split
-#[pylithapp.timedependent.formulation]
-#split_fields = True
-#use_custom_constraint_pc = True
-#matrix_type = aij
-
-#[pylithapp.petsc]
-#fs_pc_type = fieldsplit
-#fs_pc_fieldsplit_real_diagonal =
-#fs_pc_fieldsplit_type = multiplicative
-#fs_fieldsplit_0_pc_type = ml
-#fs_fieldsplit_1_pc_type = ml
-#fs_fieldsplit_2_pc_type = ilu
-#fs_fieldsplit_2_pc_type = jacobi
-#fs_fieldsplit_0_ksp_type = gmres
-#fs_fieldsplit_1_ksp_type = gmres
-#fs_fieldsplit_2_ksp_type = richardson
-
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2010-05-29 02:23:03 UTC (rev 16826)
@@ -75,6 +75,7 @@
feassemble/IntegratorElasticity.cc \
feassemble/ElasticityImplicit.cc \
feassemble/ElasticityExplicit.cc \
+ feassemble/ElasticityExplicitTri3.cc \
feassemble/ElasticityExplicitTet4.cc \
feassemble/IntegratorElasticityLgDeform.cc \
feassemble/ElasticityImplicitLgDeform.cc \
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2010-05-29 02:23:03 UTC (rev 16826)
@@ -319,6 +319,7 @@
// Compute B(transpose) * sigma, first computing strains
calcTotalStrainFn(&strainCell, basisDeriv, dispCell,
numBasis, numQuadPts);
+
const double_array& stressCell = _material->calcStress(strainCell, true);
#if defined(DETAILED_EVENT_LOGGING)
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc 2010-05-29 02:23:03 UTC (rev 16826)
@@ -249,7 +249,8 @@
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];
+ quadPtsGlobal[iDim] +=
+ coordinatesCell[iBasis*spaceDim+iDim] / numBasis;
_normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
lengthScale);
@@ -545,7 +546,8 @@
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];
+ quadPtsGlobal[iDim] +=
+ coordinatesCell[iBasis*spaceDim+iDim] / numBasis;
_normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
lengthScale);
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.hh 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.hh 2010-05-29 02:23:03 UTC (rev 16826)
@@ -143,7 +143,7 @@
*/
double _volume(const double_array& coordinatesCell) const;
- /** Compute volume of tetrahedral cell.
+ /** Compute derivatives of basis functions of tetrahedral cell.
*
* @param coordinatesCell Coordinates of vertices of cell.
* @returns Derivatives of basis functions.
Added: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.cc 2010-05-29 02:23:03 UTC (rev 16826)
@@ -0,0 +1,915 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "ElasticityExplicitTri3.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;
+typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Mesh::UpdateAddVisitor UpdateAddVisitor;
+
+// ----------------------------------------------------------------------
+const int pylith::feassemble::ElasticityExplicitTri3::_spaceDim = 2;
+const int pylith::feassemble::ElasticityExplicitTri3::_cellDim = 2;
+const int pylith::feassemble::ElasticityExplicitTri3::_tensorSize = 3;
+const int pylith::feassemble::ElasticityExplicitTri3::_numBasis = 3;
+const int pylith::feassemble::ElasticityExplicitTri3::_numQuadPts = 1;
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::ElasticityExplicitTri3::ElasticityExplicitTri3(void) :
+ _dtm1(-1.0)
+{ // constructor
+ _basisDerivArray.resize(_numQuadPts*_numBasis*_spaceDim);
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::ElasticityExplicitTri3::~ElasticityExplicitTri3(void)
+{ // destructor
+ deallocate();
+} // destructor
+
+// ----------------------------------------------------------------------
+// Deallocate PETSc and local data structures.
+void
+pylith::feassemble::ElasticityExplicitTri3::deallocate(void)
+{ // deallocate
+ IntegratorElasticity::deallocate();
+} // deallocate
+
+// ----------------------------------------------------------------------
+// Set time step for advancing from time t to time t+dt.
+void
+pylith::feassemble::ElasticityExplicitTri3::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::ElasticityExplicitTri3::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::ElasticityExplicitTri3::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 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
+ double_array accCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& accSection =
+ fields->get("acceleration(t)").section();
+ assert(!accSection.isNull());
+ RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
+
+ double_array dispCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispSection =
+ fields->get("disp(t)").section();
+ assert(!dispSection.isNull());
+ RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
+
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->lengthScale();
+ const double gravityScale =
+ _normalizer->pressureScale() / (_normalizer->lengthScale() *
+ _normalizer->densityScale());
+
+ _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 area = _area(coordinatesCell);
+ assert(area > 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
+ accVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, accVisitor);
+
+ dispVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispVisitor);
+
+#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] +=
+ coordinatesCell[iBasis*spaceDim+iDim] / numBasis;
+ _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] * area / 3.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] * area / 36.0;
+ 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 * accCell[jBasis*spaceDim+iDim];
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(3 + numBasis*numBasis*spaceDim*2);
+ _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 x1 = coordinatesCell[2];
+ const double y1 = coordinatesCell[3];
+
+ const double x2 = coordinatesCell[4];
+ const double y2 = coordinatesCell[5];
+
+ const double scaleB = 2.0 * area;
+ const double b0 = (y1 - y2) / scaleB;
+ const double c0 = (x2 - x1) / scaleB;
+
+ const double b1 = (y2 - y0) / scaleB;
+ const double c1 = (x0 - x2) / scaleB;
+
+ const double b2 = (y0 - y1) / scaleB;
+ const double c2 = (x1 - x0) / scaleB;
+
+ assert(strainCell.size() == 3);
+ strainCell[0] = b2*dispCell[4] + b1*dispCell[2] + b0*dispCell[0];
+
+ strainCell[1] = c2*dispCell[5] + c1*dispCell[3] + c0*dispCell[1];
+
+ strainCell[2] = (b2*dispCell[5] + c2*dispCell[4] + b1*dispCell[3] +
+ c1*dispCell[2] + b0*dispCell[1] + c0*dispCell[0]) / 2.0;
+
+
+ const double_array& stressCell = _material->calcStress(strainCell, true);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(34);
+ _logger->eventEnd(stressEvent);
+ _logger->eventBegin(computeEvent);
+#endif
+
+ assert(_cellVector.size() == 6);
+ assert(stressCell.size() == 3);
+ _cellVector[0] -= (c0*stressCell[2] + b0*stressCell[0]) * area;
+ _cellVector[1] -= (b0*stressCell[2] + c0*stressCell[1]) * area;
+ _cellVector[2] -= (c1*stressCell[2] + b1*stressCell[0]) * area;
+ _cellVector[3] -= (b1*stressCell[2] + c1*stressCell[1]) * area;
+ _cellVector[4] -= (c2*stressCell[2] + b2*stressCell[0]) * area;
+ _cellVector[5] -= (b2*stressCell[2] + c2*stressCell[1]) * area;
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(30);
+ _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()*(3 + numBasis*numBasis*spaceDim*2 + 34+30));
+ _logger->eventEnd(computeEvent);
+#endif
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Integrate constributions to residual term (r) for operator.
+void
+pylith::feassemble::ElasticityExplicitTri3::integrateResidualLumped(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateResidualLumped
+ /// Member prototype for _elasticityResidualXD()
+ typedef void (pylith::feassemble::ElasticityExplicitTri3::*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 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();
+
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ assert(!sieve.isNull());
+
+ // Get sections
+ double_array accCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& accSection =
+ fields->get("acceleration(t)").section();
+ assert(!accSection.isNull());
+ RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
+
+ double_array dispCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispSection =
+ fields->get("disp(t)").section();
+ assert(!dispSection.isNull());
+ RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ 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.
+ 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
+#if 1
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+
+ accVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, accVisitor);
+
+ dispVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispVisitor);
+#else
+ coordsVisitor.clear();
+ sieve->orientedConeOpt(*c_iter, coordsVisitor, numBasis, spaceDim);
+
+ accVisitor.clear();
+ sieve->orientedConeOpt(*c_iter, accVisitor, numBasis, spaceDim);
+
+ dispVisitor.clear();
+ sieve->orientedConeOpt(*c_iter, dispVisitor, numBasis, spaceDim);
+#endif
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(restrictEvent);
+ _logger->eventBegin(geometryEvent);
+#endif
+
+ // Compute geometry information for current cell
+ const double area = _area(coordinatesCell);
+ assert(area > 0.0);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(geometryEvent);
+ _logger->eventBegin(stateVarsEvent);
+#endif
+
+ // Get state variables for cell.
+ _material->retrievePropsAndVars(*c_iter);
+
+ // Get density at quadrature points for this cell
+ const double_array& density = _material->calcDensity();
+
+#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] +=
+ coordinatesCell[iBasis*spaceDim+iDim] / numBasis;
+ _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] * area / 3.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] * area / 3.0;
+ _cellVector -= wtVertex * accCell;
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(2 + numBasis*spaceDim*2);
+ _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 x1 = coordinatesCell[2];
+ const double y1 = coordinatesCell[3];
+
+ const double x2 = coordinatesCell[4];
+ const double y2 = coordinatesCell[5];
+
+ const double scaleB = 2.0 * area;
+ const double b0 = (y1 - y2) / scaleB;
+ const double c0 = (x2 - x1) / scaleB;
+
+ const double b1 = (y2 - y0) / scaleB;
+ const double c1 = (x0 - x2) / scaleB;
+
+ const double b2 = (y0 - y1) / scaleB;
+ const double c2 = (x1 - x0) / scaleB;
+
+ assert(strainCell.size() == 3);
+ strainCell[0] = b2*dispCell[4] + b1*dispCell[2] + b0*dispCell[0];
+
+ strainCell[1] = c2*dispCell[5] + c1*dispCell[3] + c0*dispCell[1];
+
+ strainCell[2] = (b2*dispCell[5] + c2*dispCell[4] + b1*dispCell[3] +
+ c1*dispCell[2] + b0*dispCell[1] + c0*dispCell[0]) / 2.0;
+
+ const double_array& stressCell = _material->calcStress(strainCell, true);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(34);
+ _logger->eventEnd(stressEvent);
+ _logger->eventBegin(computeEvent);
+#endif
+
+ assert(_cellVector.size() == 6);
+ assert(stressCell.size() == 3);
+ _cellVector[0] -= (c0*stressCell[2] + b0*stressCell[0]) * area;
+ _cellVector[1] -= (b0*stressCell[2] + c0*stressCell[1]) * area;
+ _cellVector[2] -= (c1*stressCell[2] + b1*stressCell[0]) * area;
+ _cellVector[3] -= (b1*stressCell[2] + c1*stressCell[1]) * area;
+ _cellVector[4] -= (c2*stressCell[2] + b2*stressCell[0]) * area;
+ _cellVector[5] -= (b2*stressCell[2] + c2*stressCell[1]) * area;
+
+#if defined(DETAILED_EVENT_LOGGING)
+ PetscLogFlops(30);
+ _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*2 + 34+30));
+ _logger->eventEnd(computeEvent);
+#endif
+} // integrateResidualLumped
+
+// ----------------------------------------------------------------------
+// Compute matrix associated with operator.
+void
+pylith::feassemble::ElasticityExplicitTri3::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.");
+
+ // 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
+ double_array dispCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispSection =
+ fields->get("disp(t)").section();
+ assert(!dispSection.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", dispSection);
+ assert(!globalOrder.isNull());
+ // We would need to request unique points here if we had an interpolated mesh
+ topology::Mesh::IndicesVisitor jacobianVisitor(*dispSection, *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 area = _area(coordinatesCell);
+ assert(area > 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] * area / (36.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::ElasticityExplicitTri3::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 area = _area(coordinatesCell);
+ assert(area > 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
+
+ // Compute Jacobian for inertial terms
+ const double_array& density = _material->calcDensity();
+ _cellVector = density[0] * area / (3.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 area of triangular cell.
+double
+pylith::feassemble::ElasticityExplicitTri3::_area(
+ const double_array& coordinatesCell) const
+{ // __area
+ assert(6 == coordinatesCell.size());
+
+ const double x0 = coordinatesCell[0];
+ const double y0 = coordinatesCell[1];
+
+ const double x1 = coordinatesCell[2];
+ const double y1 = coordinatesCell[3];
+
+ const double x2 = coordinatesCell[4];
+ const double y2 = coordinatesCell[5];
+
+ const double area = 0.5*((x1-x0)*(y2-y0) - (x2-x0)*(y1-y0));
+ PetscLogFlops(8);
+
+ return area;
+} // _area
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.hh 2010-05-29 02:23:03 UTC (rev 16826)
@@ -0,0 +1,180 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/feassemble/ElasticityExplicitTri3.hh
+ *
+ * @brief Explicit time integration of dynamic elasticity equation
+ * using linear triangular finite-elements.
+ */
+
+#if !defined(pylith_feassemble_elasticityexplicittri3_hh)
+#define pylith_feassemble_elasticityexplicittri3_hh
+
+// Include directives ---------------------------------------------------
+#include "IntegratorElasticity.hh" // ISA IntegratorElasticity
+
+// ElasticityExplicitTri3 ---------------------------------------------------
+/**@brief Explicit time integration of the dynamic elasticity equation
+ * using linear triangular 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::ElasticityExplicitTri3 : public IntegratorElasticity
+{ // ElasticityExplicitTri3
+ friend class TestElasticityExplicitTri3; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ ElasticityExplicitTri3(void);
+
+ /// Destructor
+ ~ElasticityExplicitTri3(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 area of triangular cell.
+ *
+ * @param coordinatesCell Coordinates of vertices of cell.
+ * @returns Area of cell.
+ */
+ double _area(const double_array& coordinatesCell) const;
+
+ /** Compute derivatives of basis functions of triangular 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.
+ ElasticityExplicitTri3(const ElasticityExplicitTri3&);
+
+ /// Not implemented
+ const ElasticityExplicitTri3& operator=(const ElasticityExplicitTri3&);
+
+}; // ElasticityExplicitTri3
+
+#endif // pylith_feassemble_elasticityexplicittri3_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2010-05-29 02:23:03 UTC (rev 16826)
@@ -19,6 +19,7 @@
Constraint.hh \
Constraint.icc \
ElasticityExplicit.hh \
+ ElasticityExplicitTri3.hh \
ElasticityExplicitTet4.hh \
ElasticityExplicitLgDeform.hh \
ElasticityImplicit.hh \
Modified: short/3D/PyLith/trunk/libsrc/feassemble/feassemblefwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/feassemblefwd.hh 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/libsrc/feassemble/feassemblefwd.hh 2010-05-29 02:23:03 UTC (rev 16826)
@@ -57,6 +57,7 @@
class ElasticityExplicit;
class ElasticityExplicitTet4;
+ class ElasticityExplicitTri3;
class IntegratorElasticityLgDeform;
class ElasticityImplicitLgDeform;
Added: short/3D/PyLith/trunk/libsrc/feassemble/tet4_elasticity.wxm
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/tet4_elasticity.wxm (rev 0)
+++ short/3D/PyLith/trunk/libsrc/feassemble/tet4_elasticity.wxm 2010-05-29 02:23:03 UTC (rev 16826)
@@ -0,0 +1,54 @@
+/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
+/* [ Created with wxMaxima version 0.8.2 ] */
+
+/* [wxMaxima: input start ] */
+b : matrix([1,y1,z1],[1,y2,z2],[1,y3,z3]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+-determinant(b);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+c : matrix([x1, 1, z1],[x2,1,z2],[x3,1,z3]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+-determinant(c);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+d : matrix([x0, y0,1],[x1,y1,1],[x2,y2,1]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+determinant(d);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+B : matrix([b1,0,0,b2,0,0,b3,0,0,b4,0,0],
+[0,c1,0,0,c2,0,0,c3,0,0,c4,0],
+[0,0,d1,0,0,d2,0,0,d3,0,0,d4],
+[c1,b1,0,c2,b2,0,c3,b3,0,c4,b4,0],
+[0,d1,c1,0,d2,c2,0,d3,c3,0,d4,c4],
+[d1,0,b1,d2,0,b2,d3,0,b3,d4,0,b4]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+disp : matrix([dispTCell0,dispTCell1,dispTCell2,dispTCell3,dispTCell4,dispTCell5,dispTCell6,dispTCell7,dispTCell8,dispTCell9,dispTCell10,dispTCell11]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+B . transpose(disp);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+stress : matrix([stressCell0, stressCell1, stressCell2, stressCell3, stressCell4, stressCell5]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+transpose(B) . transpose(stress);
+/* [wxMaxima: input end ] */
+
+/* Maxima can't load/batch files which end with a comment! */
+"Created with wxMaxima"$
Added: short/3D/PyLith/trunk/libsrc/feassemble/tri3_elasticity.wxm
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/tri3_elasticity.wxm (rev 0)
+++ short/3D/PyLith/trunk/libsrc/feassemble/tri3_elasticity.wxm 2010-05-29 02:23:03 UTC (rev 16826)
@@ -0,0 +1,27 @@
+/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
+/* [ Created with wxMaxima version 0.8.2 ] */
+
+/* [wxMaxima: input start ] */
+B : matrix([b0,0,b1,0,b2,0],
+[0,c0,0,c1,0,c2],
+[c0,b0,c1,b1,c2,b2]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+disp : matrix([dispCell0,dispCell1,dispCell2,dispCell3,dispCell4,dispCell5]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+B . transpose(disp);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+stress : matrix([stressCell0, stressCell1, stressCell2]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+transpose(B) . transpose(stress);
+/* [wxMaxima: input end ] */
+
+/* Maxima can't load/batch files which end with a comment! */
+"Created with wxMaxima"$
Added: short/3D/PyLith/trunk/modulesrc/feassemble/ElasticityExplicitTri3.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/ElasticityExplicitTri3.i (rev 0)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/ElasticityExplicitTri3.i 2010-05-29 02:23:03 UTC (rev 16826)
@@ -0,0 +1,98 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file modulesrc/feassemble/ElasticityExplicitTri3.i
+ *
+ * @brief Python interface to C++ ElasticityExplicitTri3 object.
+ */
+
+namespace pylith {
+ namespace feassemble {
+
+ class ElasticityExplicitTri3 : public IntegratorElasticity
+ { // ElasticityExplicitTri3
+
+ // PUBLIC MEMBERS /////////////////////////////////////////////////
+ public :
+
+ /// Constructor
+ ElasticityExplicitTri3(void);
+
+ /// Destructor
+ ~ElasticityExplicitTri3(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 residual term (r) for operator.
+ *
+ * @param residual Field containing values for residual
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateResidualLumped(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);
+
+ }; // ElasticityExplicitTri3
+
+ } // feassemble
+} // pylith
+
+
+// End of file
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/Makefile.am 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/Makefile.am 2010-05-29 02:23:03 UTC (rev 16826)
@@ -38,6 +38,7 @@
IntegratorElasticity.i \
ElasticityImplicit.i \
ElasticityExplicit.i \
+ ElasticityExplicitTri3.i \
ElasticityExplicitTet4.i \
IntegratorElasticityLgDeform.i \
ElasticityImplicitLgDeform.i \
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.i 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.i 2010-05-29 02:23:03 UTC (rev 16826)
@@ -35,6 +35,7 @@
#include "pylith/feassemble/Quadrature.hh"
#include "pylith/feassemble/ElasticityImplicit.hh"
#include "pylith/feassemble/ElasticityExplicit.hh"
+#include "pylith/feassemble/ElasticityExplicitTri3.hh"
#include "pylith/feassemble/ElasticityExplicitTet4.hh"
#include "pylith/feassemble/ElasticityImplicitLgDeform.hh"
#include "pylith/feassemble/ElasticityExplicitLgDeform.hh"
@@ -85,6 +86,7 @@
%include "ElasticityImplicit.i"
%include "ElasticityExplicit.i"
%include "ElasticityExplicitTet4.i"
+%include "ElasticityExplicitTri3.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-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/pylith/Makefile.am 2010-05-29 02:23:03 UTC (rev 16826)
@@ -44,6 +44,7 @@
feassemble/Constraint.py \
feassemble/ElasticityExplicit.py \
feassemble/ElasticityExplicitTet4.py \
+ feassemble/ElasticityExplicitTri3.py \
feassemble/ElasticityExplicitLgDeform.py \
feassemble/ElasticityImplicit.py \
feassemble/ElasticityImplicitLgDeform.py \
@@ -118,6 +119,7 @@
problems/__init__.py \
problems/Explicit.py \
problems/ExplicitLumped.py \
+ problems/ExplicitLumpedTri3.py \
problems/ExplicitLumpedTet4.py \
problems/ExplicitLgDeform.py \
problems/Formulation.py \
Modified: short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicitTet4.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicitTet4.py 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicitTet4.py 2010-05-29 02:23:03 UTC (rev 16826)
@@ -29,7 +29,7 @@
# PUBLIC METHODS /////////////////////////////////////////////////////
- def __init__(self, name="elasticityexplicit"):
+ def __init__(self, name="elasticityexplicittet4"):
"""
Constructor.
"""
Added: short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicitTri3.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicitTri3.py (rev 0)
+++ short/3D/PyLith/trunk/pylith/feassemble/ElasticityExplicitTri3.py 2010-05-29 02:23:03 UTC (rev 16826)
@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/feassemble/ElasticityExplicitTri3.py
+##
+## @brief Python object for explicit time integration of dynamic
+## elasticity equation using finite-elements.
+##
+## Factory: integrator
+
+from IntegratorElasticity import IntegratorElasticity
+from feassemble import ElasticityExplicitTri3 as ModuleElasticityExplicitTri3
+
+# ElasticityExplicitTri3 class
+class ElasticityExplicitTri3(IntegratorElasticity, ModuleElasticityExplicitTri3):
+ """
+ Python object for explicit time integration of dynamic elasticity
+ equation using finite-elements.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="elasticityexplicittri3"):
+ """
+ Constructor.
+ """
+ IntegratorElasticity.__init__(self, name)
+ ModuleElasticityExplicitTri3.__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)
+ ModuleElasticityExplicitTri3.initialize(self, self.mesh)
+ self._initializeOutput(totalTime, numTimeSteps, normalizer)
+
+ self._eventLogger.eventEnd(logEvent)
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def integrator():
+ """
+ Factory associated with ElasticityExplicitTri3.
+ """
+ return ElasticityExplicitTri3()
+
+
+# End of file
Added: short/3D/PyLith/trunk/pylith/problems/ExplicitLumpedTri3.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/ExplicitLumpedTri3.py (rev 0)
+++ short/3D/PyLith/trunk/pylith/problems/ExplicitLumpedTri3.py 2010-05-29 02:23:03 UTC (rev 16826)
@@ -0,0 +1,80 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/problems/ExplicitLumpedTri3.py
+##
+## @brief Python ExplicitLumpedTri3 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
+
+# ExplicitLumpedTri3 class
+class ExplicitLumpedTri3(ExplicitLumped):
+ """
+ Python ExplicitLumpedTri3 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="explicitlumpedtri3"):
+ """
+ Constructor.
+ """
+ ExplicitLumped.__init__(self, name)
+ return
+
+
+ def elasticityIntegrator(self):
+ """
+ Get integrator for elastic material.
+ """
+ from pylith.feassemble.ElasticityExplicitTri3 import ElasticityExplicitTri3
+ return ElasticityExplicitTri3()
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based using inventory.
+ """
+ ExplicitLumped._configure(self)
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def pde_formulation():
+ """
+ Factory associated with ExplicitLumped.
+ """
+ return ExplicitLumpedTri3()
+
+
+# End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2010-05-29 02:23:03 UTC (rev 16826)
@@ -59,6 +59,7 @@
TestElasticityExplicitGrav2DQuadratic.cc \
TestElasticityExplicitGrav3DLinear.cc \
TestElasticityExplicitGrav3DQuadratic.cc \
+ TestElasticityExplicitTri3.cc \
TestElasticityExplicitTet4.cc \
TestElasticityImplicit.cc \
TestElasticityImplicit1DLinear.cc \
@@ -119,6 +120,7 @@
TestElasticityExplicitGrav2DQuadratic.hh \
TestElasticityExplicitGrav3DLinear.hh \
TestElasticityExplicitGrav3DQuadratic.hh \
+ TestElasticityExplicitTri3.hh \
TestElasticityExplicitTet4.hh \
TestElasticityImplicit.hh \
TestElasticityImplicit1DLinear.hh \
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc 2010-05-29 01:10:33 UTC (rev 16825)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc 2010-05-29 02:23:03 UTC (rev 16826)
@@ -231,7 +231,7 @@
const int size = residualSection->sizeWithBC();
CPPUNIT_ASSERT_EQUAL(sizeE, size);
-#if 1 // DEBUGGING
+#if 0 // DEBUGGING
residual.view("RESIDUAL");
std::cout << "EXPECTED RESIDUAL" << std::endl;
for (int i=0; i < size; ++i)
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc 2010-05-29 02:23:03 UTC (rev 16826)
@@ -0,0 +1,536 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestElasticityExplicitTri3.hh" // Implementation of class methods
+
+#include "pylith/feassemble/ElasticityExplicitTri3.hh" // USES ElasticityExplicitTri3
+#include "data/ElasticityExplicitData2DLinear.hh"
+#include "pylith/feassemble/GeometryTri2D.hh" // USES GeometryTri2D
+
+#include "pylith/materials/ElasticPlaneStrain.hh" // USES ElasticPlaneStrain
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#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::TestElasticityExplicitTri3 );
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::feassemble::TestElasticityExplicitTri3::setUp(void)
+{ // setUp
+ _quadrature = new Quadrature<topology::Mesh>();
+ CPPUNIT_ASSERT(0 != _quadrature);
+ GeometryTri2D geometry;
+ _quadrature->refGeometry(&geometry);
+
+ _data = new ElasticityExplicitData2DLinear();
+ CPPUNIT_ASSERT(0 != _data);
+ _material = new materials::ElasticPlaneStrain;
+ CPPUNIT_ASSERT(0 != _material);
+ CPPUNIT_ASSERT_EQUAL(std::string("ElasticPlaneStrain"),
+ std::string(_data->matType));
+ _gravityField = 0;
+} // setUp
+
+// ----------------------------------------------------------------------
+// Tear down testing data.
+void
+pylith::feassemble::TestElasticityExplicitTri3::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::TestElasticityExplicitTri3::testConstructor(void)
+{ // testConstructor
+ ElasticityExplicitTri3 integrator;
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test timeStep().
+void
+pylith::feassemble::TestElasticityExplicitTri3::testTimeStep(void)
+{ // testTimeStep
+ ElasticityExplicitTri3 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::TestElasticityExplicitTri3::testMaterial(void)
+{ // testMaterial
+ ElasticityExplicitTri3 integrator;
+
+ materials::ElasticPlaneStrain 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::TestElasticityExplicitTri3::testNeedNewJacobian(void)
+{ // testNeedNewJacobian
+ ElasticityExplicitTri3 integrator;
+
+ materials::ElasticPlaneStrain 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::TestElasticityExplicitTri3::testUseSolnIncr(void)
+{ // testUseSolnIncr
+ ElasticityExplicitTri3 integrator;
+
+ materials::ElasticPlaneStrain 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::TestElasticityExplicitTri3::testInitialize(void)
+{ // testInitialize
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ ElasticityExplicitTri3 integrator;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+
+} // testInitialize
+
+// ----------------------------------------------------------------------
+// Test integrateResidual().
+void
+pylith::feassemble::TestElasticityExplicitTri3::testIntegrateResidual(void)
+{ // testIntegrateResidual
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ ElasticityExplicitTri3 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 1 // 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 integrateResidual().
+void
+pylith::feassemble::TestElasticityExplicitTri3::testIntegrateResidualLumped(void)
+{ // testIntegrateResidualLumped
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ ElasticityExplicitTri3 integrator;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+
+ topology::Field<topology::Mesh>& residual = fields.get("residual");
+ const double t = 1.0;
+ integrator.integrateResidualLumped(residual, t, &fields);
+
+ const double* valsE = _data->valsResidualLumped;
+ 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 1 // 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);
+} // testIntegrateResidualLumped
+
+// ----------------------------------------------------------------------
+// Test integrateJacobian().
+void
+pylith::feassemble::TestElasticityExplicitTri3::testIntegrateJacobian(void)
+{ // testIntegrateJacobian
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ ElasticityExplicitTri3 integrator;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+ integrator._needNewJacobian = true;
+
+ topology::Jacobian jacobian(fields.solution());
+
+ 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::TestElasticityExplicitTri3::testIntegrateJacobianLumped(void)
+{ // testIntegrateJacobian
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ ElasticityExplicitTri3 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::TestElasticityExplicitTri3::testUpdateStateVars(void)
+{ // testUpdateStateVars
+ CPPUNIT_ASSERT(0 != _data);
+
+ topology::Mesh mesh;
+ ElasticityExplicitTri3 integrator;
+ topology::SolutionFields fields(mesh);
+ _initialize(&mesh, &integrator, &fields);
+
+ const double t = 1.0;
+ integrator.updateStateVars(t, &fields);
+} // testUpdateStateVars
+
+// ----------------------------------------------------------------------
+// Test StableTimeStep().
+void
+pylith::feassemble::TestElasticityExplicitTri3::testStableTimeStep(void)
+{ // testStableTimeStep
+ topology::Mesh mesh;
+ ElasticityExplicitTri3 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::TestElasticityExplicitTri3::_initialize(
+ topology::Mesh* mesh,
+ ElasticityExplicitTri3* 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);
+
+ const int spaceDim = _data->spaceDim;
+ const double dt = _data->dt;
+
+ // 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, 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,
+ spaceDim);
+
+ spatialdata::units::Nondimensional normalizer;
+ spatialdata::geocoords::CSCart cs;
+ cs.setSpaceDim(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->add("acceleration(t)", "acceleration");
+ fields->solutionName("dispIncr(t->t+dt)");
+
+ topology::Field<topology::Mesh>& residual = fields->get("residual");
+ residual.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
+ residual.allocate();
+ residual.zero();
+ fields->copyLayout("residual");
+
+ const int fieldSize = 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 =
+ fields->get("dispIncr(t->t+dt)").section();
+ const ALE::Obj<RealSection>& dispTSection =
+ fields->get("disp(t)").section();
+ const ALE::Obj<RealSection>& dispTmdtSection =
+ fields->get("disp(t-dt)").section();
+ const ALE::Obj<RealSection>& accSection =
+ fields->get("acceleration(t)").section();
+ CPPUNIT_ASSERT(!dispIncrSection.isNull());
+ CPPUNIT_ASSERT(!dispTSection.isNull());
+ CPPUNIT_ASSERT(!dispTmdtSection.isNull());
+ CPPUNIT_ASSERT(!accSection.isNull());
+
+ double_array accVertex(spaceDim);
+
+ const int offset = _data->numCells;
+ for (int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
+ dispIncrSection->updatePoint(iVertex+offset,
+ &_data->fieldTIncr[iVertex*spaceDim]);
+ dispTSection->updatePoint(iVertex+offset,
+ &_data->fieldT[iVertex*spaceDim]);
+ dispTmdtSection->updatePoint(iVertex+offset,
+ &_data->fieldTmdt[iVertex*spaceDim]);
+
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ accVertex[iDim] = (_data->fieldTIncr[iVertex*spaceDim+iDim] -
+ _data->fieldT[iVertex*spaceDim+iDim] +
+ _data->fieldTmdt[iVertex*spaceDim+iDim]) / (dt*dt);
+ accSection->updatePoint(iVertex+offset, &accVertex[0]);
+ } // for
+} // _initialize
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.hh 2010-05-29 02:23:03 UTC (rev 16826)
@@ -0,0 +1,133 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/feassemble/TestElasticityExplicitTri3.hh
+ *
+ * @brief C++ TestElasticityExplicitTri3 object
+ *
+ * C++ unit testing for ElasticityExplicitTri3.
+ */
+
+#if !defined(pylith_feassemble_testelasticityexplicittri3_hh)
+#define pylith_feassemble_testelasticityexplicittri3_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 TestElasticityExplicitTri3;
+ class ElasticityExplicitData;
+ } // feassemble
+} // pylith
+
+/// C++ unit testing for ElasticityExplicitTri3
+class pylith::feassemble::TestElasticityExplicitTri3 : public CppUnit::TestFixture
+{ // class TestElasticityExplicitTri3
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestElasticityExplicitTri3 );
+
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testTimeStep );
+ CPPUNIT_TEST( testMaterial );
+ CPPUNIT_TEST( testNeedNewJacobian );
+ CPPUNIT_TEST( testUseSolnIncr );
+ CPPUNIT_TEST( testInitialize );
+ CPPUNIT_TEST( testIntegrateResidual );
+ CPPUNIT_TEST( testIntegrateResidualLumped );
+ 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 integrateResidualLumped().
+ void testIntegrateResidualLumped(void);
+
+ /// Test integrateJacobian().
+ void testIntegrateJacobian(void);
+
+ /// Test integrateJacobianLumped().
+ void testIntegrateJacobianLumped(void);
+
+ /// Test updateStateVars().
+ void testUpdateStateVars(void);
+
+ /// Test StableTimeStep().
+ void testStableTimeStep(void);
+
+ // PROTECTED MEMBERS //////////////////////////////////////////////////
+protected :
+
+ ElasticityExplicitData* _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,
+ ElasticityExplicitTri3* const integrator,
+ topology::SolutionFields* const fields);
+
+}; // class TestElasticityExplicitTri3
+
+#endif // pylith_feassemble_testelasticityexplicittri3_hh
+
+
+// End of file
More information about the CIG-COMMITS
mailing list