[cig-commits] r16227 - in short/3D/PyLith/trunk: . libsrc/bc libsrc/faults libsrc/feassemble libsrc/materials libsrc/problems modulesrc/problems pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Wed Feb 3 16:31:00 PST 2010
Author: brad
Date: 2010-02-03 16:30:58 -0800 (Wed, 03 Feb 2010)
New Revision: 16227
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/materials/Makefile.am
short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
short/3D/PyLith/trunk/libsrc/problems/Formulation.hh
short/3D/PyLith/trunk/modulesrc/problems/Formulation.i
short/3D/PyLith/trunk/pylith/problems/ExplicitLumped.py
Log:
More optimization for lumped solve. Also added defines to control level of logging detail.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/TODO 2010-02-04 00:30:58 UTC (rev 16227)
@@ -65,6 +65,12 @@
Unit tests
+ setting up cohesive info
+
+ ElasticityExplicit
+ integrateResidualLumped()
+
+
FaultCohesiveDynL
integrateJacobian (lumped)
adjustSolnLumped
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2010-02-04 00:30:58 UTC (rev 16227)
@@ -30,6 +30,7 @@
#include <sstream> // USES std::ostringstream
//#define PRECOMPUTE_GEOMETRY
+//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
@@ -304,18 +305,19 @@
topology::Mesh::RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
numBasis*spaceDim);
- _logger->eventEnd(setupEvent);
-
#if !defined(PRECOMPUTE_GEOMETRY)
- _logger->eventBegin(geometryEvent);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveSubMesh->getRealSection("coordinates");
RestrictVisitor coordsVisitor(*coordinates,
coordinatesCell.size(), &coordinatesCell[0]);
- _logger->eventEnd(geometryEvent);
#endif
+ _logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
+
// Get parameters used in integration.
const double dt = _dt;
assert(dt > 0);
@@ -324,7 +326,9 @@
c_iter != cellsEnd;
++c_iter) {
// Get geometry information for current cell
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
+#endif
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
@@ -332,13 +336,16 @@
sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
_quadrature->computeGeometry(coordinatesCell, *c_iter);
#endif
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
+ _logger->eventBegin(restrictEvent);
+#endif
// Reset element vector to zero
_resetCellVector();
// Restrict input fields to cell
- _logger->eventBegin(restrictEvent);
dispTVisitor.clear();
sieveSubMesh->restrictClosure(*c_iter, dispTVisitor);
const double* dispTCell = dispTVisitor.getValues();
@@ -349,8 +356,10 @@
dampersSection->restrictPoint(*c_iter, &dampersCell[0], dampersCell.size());
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
+#endif
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -373,16 +382,25 @@
} // for
} // for
} // for
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
// Assemble cell contribution into field
- _logger->eventBegin(updateEvent);
residualVisitor.clear();
sieveSubMesh->updateClosure(*c_iter, residualVisitor);
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
+#endif
} // for
+ PetscLogFlops(cells->size()*numQuadPts*(3+numBasis*(1+numBasis*(5*spaceDim))));
- PetscLogFlops(cells->size()*numQuadPts*(3+numBasis*(1+numBasis*(5*spaceDim))));
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+#endif
} // integrateResidual
// ----------------------------------------------------------------------
@@ -441,8 +459,6 @@
(int) pow(sieveMesh->getSieve()->getMaxConeSize(),
sieveMesh->depth())*spaceDim);
- _logger->eventEnd(setupEvent);
-
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
assert(0 != jacobianMat);
@@ -455,20 +471,25 @@
_initCellMatrix();
#if !defined(PRECOMPUTE_GEOMETRY)
- _logger->eventBegin(geometryEvent);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveSubMesh->getRealSection("coordinates");
RestrictVisitor coordsVisitor(*coordinates,
coordinatesCell.size(), &coordinatesCell[0]);
- _logger->eventEnd(geometryEvent);
#endif
+ _logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
+
for (SieveSubMesh::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
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
@@ -476,15 +497,19 @@
sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
_quadrature->computeGeometry(coordinatesCell, *c_iter);
#endif
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
+ _logger->eventBegin(restrictEvent);
+#endif
// Get damping constants
- _logger->eventBegin(restrictEvent);
assert(numQuadPts*spaceDim == dampersSection->getFiberDimension(*c_iter));
const double* dampingConstsCell = dampersSection->restrictPoint(*c_iter);
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
-
_logger->eventBegin(computeEvent);
+#endif
// Reset element vector to zero
_resetCellMatrix();
@@ -510,20 +535,28 @@
} // for
} // for
} // for
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
// Assemble cell contribution into PETSc Matrix
- _logger->eventBegin(updateEvent);
jacobianVisitor.clear();
PetscErrorCode err = updateOperator(jacobianMat, *sieveSubMesh->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
-
PetscLogFlops(cells->size()*numQuadPts*(3+numBasis*(1+numBasis*(1+2*spaceDim))));
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+#endif
+
_needNewJacobian = false;
} // integrateJacobian
@@ -589,23 +622,26 @@
topology::Mesh::UpdateAddVisitor jacobianVisitor(*jacobianSection,
&_cellVector[0]);
- _logger->eventEnd(setupEvent);
-
#if !defined(PRECOMPUTE_GEOMETRY)
- _logger->eventBegin(geometryEvent);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveSubMesh->getRealSection("coordinates");
RestrictVisitor coordsVisitor(*coordinates,
coordinatesCell.size(), &coordinatesCell[0]);
- _logger->eventEnd(geometryEvent);
#endif
+ _logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
+
for (SieveSubMesh::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
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
@@ -613,15 +649,19 @@
sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
_quadrature->computeGeometry(coordinatesCell, *c_iter);
#endif
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
+ _logger->eventBegin(restrictEvent);
+#endif
// Get damping constants
- _logger->eventBegin(restrictEvent);
assert(numQuadPts*spaceDim == dampersSection->getFiberDimension(*c_iter));
const double* dampingConstsCell = dampersSection->restrictPoint(*c_iter);
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
-
_logger->eventBegin(computeEvent);
+#endif
// Reset element vector to zero
_resetCellMatrix();
@@ -649,16 +689,24 @@
} // for
_lumpCellMatrix();
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
// Assemble cell contribution into lumped matrix.
- _logger->eventBegin(updateEvent);
jacobianVisitor.clear();
sieveSubMesh->updateClosure(*c_iter, jacobianVisitor);
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
+#endif
} // for
PetscLogFlops(cells->size()*numQuadPts*(3+numBasis*(1+numBasis*(1+2*spaceDim))));
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+#endif
+
_needNewJacobian = false;
} // integrateJacobian
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2010-02-04 00:30:58 UTC (rev 16227)
@@ -34,6 +34,7 @@
#include <stdexcept> // USES std::runtime_error
//#define PRECOMPUTE_GEOMETRY
+//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
@@ -246,6 +247,9 @@
assert(!dispTIncrSection.isNull());
_logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
@@ -254,7 +258,9 @@
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(restrictEvent);
+#endif
slipSection->restrictPoint(v_fault, &slipVertex[0], slipVertex.size());
@@ -271,8 +277,10 @@
dispTIncrSection->restrictPoint(v_positive, &dispTIncrVertexP[0], dispTIncrVertexP.size());
dispTIncrSection->restrictPoint(v_lagrange, &dispTIncrVertexL[0], dispTIncrVertexL.size());
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
+#endif
// Compute current estimate of displacement at time t+dt using
// solution increment.
@@ -297,8 +305,10 @@
residualVertexL[kDim] -= (dispTpdtVertexP[iDim] - dispTpdtVertexN[iDim])
* orientationVertex[kDim*spaceDim+iDim];
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
+#endif
assert(residualVertexN.size() == residualSection->getFiberDimension(v_negative));
residualSection->updateAddPoint(v_negative, &residualVertexN[0]);
@@ -309,11 +319,15 @@
assert(residualVertexL.size() == residualSection->getFiberDimension(v_lagrange));
residualSection->updateAddPoint(v_lagrange, &residualVertexL[0]);
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
-
+#endif
} // for
+ PetscLogFlops(numVertices*spaceDim*spaceDim*8);
- PetscLogFlops(numVertices*spaceDim*spaceDim*8);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+#endif
} // integrateResidualAssembled
// ----------------------------------------------------------------------
@@ -412,20 +426,27 @@
sieveMesh->depth()) * spaceDim);
_logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
for (SieveMesh::label_sequence::iterator c_iter = cellsCohesiveBegin; c_iter
!= cellsCohesiveEnd; ++c_iter) {
const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(restrictEvent);
+#endif
matrixCell = 0.0;
// Get orientations at fault cell's vertices.
orientationVisitor.clear();
faultSieveMesh->restrictClosure(c_fault, orientationVisitor);
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
+#endif
for (int iConstraint = 0; iConstraint < numConstraintVert; ++iConstraint) {
// Blocks in cell matrix associated with normal cohesive
@@ -462,18 +483,27 @@
} // for
} // for
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
// Insert cell contribution into PETSc Matrix
- _logger->eventBegin(updateEvent);
jacobianVisitor.clear();
PetscErrorCode err = updateOperator(jacobianMatrix, *sieveMesh->getSieve(),
jacobianVisitor, *c_iter, &matrixCell[0], INSERT_VALUES);
CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
-
+#endif
} // for
PetscLogFlops(cellsCohesiveSize*numConstraintVert*spaceDim*spaceDim*4);
+
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+#endif
+
_needNewJacobian = false;
} // integrateJacobianAssembled
@@ -509,6 +539,9 @@
assert(!jacobianSection.isNull());
_logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
@@ -516,13 +549,20 @@
assert(jacobianSection->getFiberDimension(v_lagrange) == spaceDim);
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(updateEvent);
+#endif
jacobianSection->updatePoint(v_lagrange, &jacobianVertex[0]);
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
+#endif
} // for
-
PetscLogFlops(0);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+#endif
+
_needNewJacobian = false;
} // integrateJacobianAssembled
@@ -599,6 +639,10 @@
_logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(computeEvent);
+#endif
+
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
@@ -606,7 +650,9 @@
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(restrictEvent);
+#endif
// Get orientations at fault cell's vertices.
orientationSection->restrictPoint(v_fault, &orientationVertex[0], orientationVertex.size());
@@ -626,8 +672,10 @@
dispTSection->restrictPoint(v_negative, &dispTVertexN[0], dispTVertexN.size());
dispTSection->restrictPoint(v_positive, &dispTVertexP[0], dispTVertexP.size());
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
+#endif
switch (spaceDim) { // switch
case 1: {
@@ -783,8 +831,10 @@
throw std::logic_error("Unknown spatial dimension.");
} // switch
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
+#endif
assert(solutionVertexN.size() == solutionSection->getFiberDimension(v_negative));
solutionSection->updateAddPoint(v_negative, &solutionVertexN[0]);
@@ -795,7 +845,9 @@
assert(solutionVertexL.size() == solutionSection->getFiberDimension(v_lagrange));
solutionSection->updateAddPoint(v_lagrange, &solutionVertexL[0]);
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
+#endif
} // for
switch(spaceDim) {
@@ -812,6 +864,10 @@
assert(0);
throw std::logic_error("Unknown spatial dimension.");
} // switch
+
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+#endif
} // adjustSolnLumped
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2010-02-04 00:30:58 UTC (rev 16227)
@@ -35,6 +35,7 @@
#include <stdexcept> // USES std::runtime_error
//#define PRECOMPUTE_GEOMETRY
+//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
@@ -211,13 +212,18 @@
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) {
// Compute geometry information for current cell
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
+#endif
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
@@ -225,23 +231,31 @@
sieveMesh->restrictClosure(*c_iter, coordsVisitor);
_quadrature->computeGeometry(coordinatesCell, *c_iter);
#endif
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
+ _logger->eventBegin(stateVarsEvent);
+#endif
// Get state variables for cell.
- _logger->eventBegin(stateVarsEvent);
_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
- _logger->eventBegin(restrictEvent);
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
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -251,7 +265,6 @@
// Compute body force vector if gravity is being used.
if (0 != _gravityField) {
- _logger->eventBegin(computeEvent);
const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();
assert(0 != cs);
@@ -279,49 +292,321 @@
} // for
} // for
PetscLogFlops(numQuadPts * (2 + numBasis * (1 + 2 * spaceDim)));
- _logger->eventEnd(computeEvent);
} // if
// Compute action for inertial terms
- _logger->eventBegin(computeEvent);
const double_array& density = _material->calcDensity();
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt =
- quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- const double valI = wt*basis[iQuad*numBasis+iBasis];
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const double valIJ = valI * basis[iQuad*numBasis+jBasis];
- for (int iDim=0; iDim < spaceDim; ++iDim)
- _cellVector[iBasis*spaceDim+iDim] +=
- valIJ * (dispTCell[jBasis*spaceDim+iDim]
- - dispTmdtCell[jBasis*spaceDim+iDim]);
+ for (int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad]
+ / dt2;
+ for (int iBasis = 0; iBasis < numBasis; ++iBasis) {
+ const double valI = wt * basis[iQuad * numBasis + iBasis];
+ for (int jBasis = 0; jBasis < numBasis; ++jBasis) {
+ const double valIJ = valI * basis[iQuad * numBasis + jBasis];
+ for (int iDim = 0; iDim < spaceDim; ++iDim)
+ _cellVector[iBasis * spaceDim + iDim] += valIJ * (dispTCell[jBasis
+ * spaceDim + iDim] - dispTmdtCell[jBasis * spaceDim + iDim]);
} // for
} // for
} // for
PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(6*spaceDim))));
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
+ _logger->eventBegin(stressEvent);
+#endif
// Compute B(transpose) * sigma, first computing strains
- _logger->eventBegin(stressEvent);
calcTotalStrainFn(&strainCell, basisDeriv, dispTCell,
numBasis, numQuadPts);
const double_array& stressCell = _material->calcStress(strainCell, true);
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stressEvent);
+ _logger->eventBegin(computeEvent);
+#endif
- _logger->eventBegin(computeEvent);
CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell);
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
// Assemble cell contribution into field
- _logger->eventBegin(updateEvent);
residualVisitor.clear();
sieveMesh->updateClosure(*c_iter, residualVisitor);
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
+#endif
} // for
+
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+#endif
} // integrateResidual
// ----------------------------------------------------------------------
+// Integrate constributions to residual term (r) for operator.
+void
+pylith::feassemble::ElasticityExplicit::integrateResidualLumped(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateResidualLumped
+ /// Member prototype for _elasticityResidualXD()
+ typedef void (pylith::feassemble::ElasticityExplicit::*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
+ 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();
+ /** :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.");
+
+ // Set variables dependent on dimension of cell
+ totalStrain_fn_type calcTotalStrainFn;
+ elasticityResidual_fn_type elasticityResidualFn;
+ if (1 == cellDim) {
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityExplicit::_elasticityResidual1D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
+ } else if (2 == cellDim) {
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityExplicit::_elasticityResidual2D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
+ } else if (3 == cellDim) {
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityExplicit::_elasticityResidual3D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D;
+ } else
+ assert(0);
+
+ // 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) {
+ // Compute geometry information for current cell
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventBegin(geometryEvent);
+#endif
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
+#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
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+ const double_array& quadPtsNondim = _quadrature->quadPts();
+
+ // Compute body force vector if gravity is being used.
+ if (0 != _gravityField) {
+ const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();
+ assert(0 != cs);
+
+ // Get density at quadrature points for this cell
+ const double_array& density = _material->calcDensity();
+
+ quadPtsGlobal = quadPtsNondim;
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
+ lengthScale);
+
+ // Compute action for element body forces
+ for (int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
+ 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 wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
+ for (int iBasis = 0, iQ = iQuad * numBasis; iBasis < numBasis; ++iBasis) {
+ const double valI = wt * basis[iQ + iBasis];
+ for (int iDim = 0; iDim < spaceDim; ++iDim) {
+ _cellVector[iBasis * spaceDim + iDim] += valI * gravVec[iDim];
+ } // for
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts * (2 + numBasis * (1 + 2 * spaceDim)));
+ } // if
+
+ // Compute action for inertial terms
+ const double_array& density = _material->calcDensity();
+ 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) {
+ const double valIJ = basis[iQ + iBasis] * valJ;
+ for (int iDim = 0; iDim < spaceDim; ++iDim) {
+ _cellVector[iBasis*spaceDim+iDim] += valIJ *
+ (dispTCell[iBasis*spaceDim+iDim] - dispTmdtCell[iBasis*spaceDim+iDim]);
+ } // for
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(4+numBasis+numBasis*(1+spaceDim*2)));
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+ _logger->eventBegin(stressEvent);
+#endif
+
+ // Compute B(transpose) * sigma, first computing strains
+ calcTotalStrainFn(&strainCell, basisDeriv, dispTCell,
+ numBasis, numQuadPts);
+ const double_array& stressCell = _material->calcStress(strainCell, true);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(stressEvent);
+ _logger->eventBegin(computeEvent);
+#endif
+
+ CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell);
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _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)
+ _logger->eventEnd(computeEvent);
+#endif
+} // integrateResidualLumped
+
+// ----------------------------------------------------------------------
// Compute matrix associated with operator.
void
pylith::feassemble::ElasticityExplicit::integrateJacobian(
@@ -400,13 +685,18 @@
&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
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
@@ -414,12 +704,19 @@
sieveMesh->restrictClosure(*c_iter, coordsVisitor);
_quadrature->computeGeometry(coordinatesCell, *c_iter);
#endif
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
+ _logger->eventBegin(stateVarsEvent);
+#endif
// Get state variables for cell.
- _logger->eventBegin(stateVarsEvent);
_material->retrievePropsAndVars(*c_iter);
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
+ _logger->eventBegin(computeEvent);
+#endif
// Reset element matrix to zero
_resetCellMatrix();
@@ -432,10 +729,9 @@
const double_array& density = _material->calcDensity();
// Compute Jacobian for inertial terms
- _logger->eventBegin(computeEvent);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt =
- quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
+ quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const double valI = wt*basis[iQ+iBasis];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
@@ -449,18 +745,28 @@
} // for
} // for
PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+spaceDim))));
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
// Assemble cell contribution into PETSc matrix.
- _logger->eventBegin(updateEvent);
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)
+ _logger->eventEnd(computeEvent);
+#endif
+
_needNewJacobian = false;
_material->resetNeedNewJacobian();
} // integrateJacobian
@@ -535,13 +841,17 @@
&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
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
@@ -549,15 +859,22 @@
sieveMesh->restrictClosure(*c_iter, coordsVisitor);
_quadrature->computeGeometry(coordinatesCell, *c_iter);
#endif
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
+ _logger->eventBegin(stateVarsEvent);
+#endif
// Get state variables for cell.
- _logger->eventBegin(stateVarsEvent);
_material->retrievePropsAndVars(*c_iter);
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
+ _logger->eventBegin(computeEvent);
+#endif
// Reset element matrix to zero
- _resetCellMatrix();
+ _resetCellVector();
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -567,33 +884,40 @@
const double_array& density = _material->calcDensity();
// Compute Jacobian for inertial terms
- _logger->eventBegin(computeEvent);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt =
- quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- const double valI = wt*basis[iQ+iBasis];
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const double valIJ = valI * basis[iQ+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] += valIJ;
- } // for
- } // for
+ 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) {
+ const double valIJ = basis[iQ + iBasis] * valJ;
+ for (int iDim = 0; iDim < spaceDim; ++iDim)
+ _cellVector[iBasis*spaceDim+iDim] += valIJ;
} // for
} // for
- PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+spaceDim))));
- _lumpCellMatrix();
+ PetscLogFlops(numQuadPts*(4+numBasis+numBasis*(1+spaceDim)));
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
// Assemble cell contribution into lumped matrix.
- _logger->eventBegin(updateEvent);
jacobianVisitor.clear();
sieveMesh->updateClosure(*c_iter, jacobianVisitor);
+
+#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
+#endif
} // for
+#if !defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+#endif
+
_needNewJacobian = false;
_material->resetNeedNewJacobian();
} // integrateJacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2010-02-04 00:30:58 UTC (rev 16227)
@@ -101,6 +101,16 @@
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.
*
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2010-02-04 00:30:58 UTC (rev 16227)
@@ -148,11 +148,34 @@
* @param t Current time
* @param fields Solution fields
*/
- virtual
+ virtual
void integrateResidualAssembled(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
+ */
+ virtual
+ void integrateResidualLumped(const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields);
+
+ /** Integrate contributions to residual term (r) for operator that
+ * do not require assembly over cells, vertices, or processors.
+ *
+ * @param residual Field containing values for residual
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ virtual
+ void integrateResidualLumpedAssembled(const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields);
+
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
*
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2010-02-04 00:30:58 UTC (rev 16227)
@@ -71,9 +71,31 @@
const topology::Field<topology::Mesh>& residual,
const double t,
topology::SolutionFields* const fields) {
- _needNewJacobian = false;
} // integrateResidualAssembled
+// Integrate contributions to residual term (r) for operator.
+template<typename quadrature_type>
+inline
+void
+pylith::feassemble::Integrator<quadrature_type>::integrateResidualLumped(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields) {
+ integrateResidual(residual, t, fields);
+} // integrateResidual
+
+// Integrate contributions to residual term (r) for operator that
+// do not require assembly over cells, vertices, or processors.
+template<typename quadrature_type>
+inline
+void
+pylith::feassemble::Integrator<quadrature_type>::integrateResidualLumpedAssembled(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields) {
+ integrateResidualAssembled(residual, t, fields);
+} // integrateResidualAssembled
+
// Integrate contributions to Jacobian matrix (A) associated with
// operator.
template<typename quadrature_type>
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2010-02-04 00:30:58 UTC (rev 16227)
@@ -189,23 +189,6 @@
} // _dimProperties
// ----------------------------------------------------------------------
-// Compute density at location from properties.
-void
-pylith::materials::ElasticIsotropic3D::_calcDensity(double* const density,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars)
-{ // _calcDensity
- assert(0 != density);
- assert(0 != properties);
- assert(_numPropsQuadPt == numProperties);
- assert(0 == numStateVars);
-
- density[0] = properties[p_density];
-} // _calcDensity
-
-// ----------------------------------------------------------------------
// Compute stress tensor at location from properties.
void
pylith::materials::ElasticIsotropic3D::_calcStress(double* const stress,
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh 2010-02-04 00:30:58 UTC (rev 16227)
@@ -196,6 +196,8 @@
}; // class ElasticIsotropic3D
+#include "ElasticIsotropic3D.icc" // inline methods
+
#endif // pylith_materials_elasticisotropic3d_hh
Modified: short/3D/PyLith/trunk/libsrc/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Makefile.am 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/libsrc/materials/Makefile.am 2010-02-04 00:30:58 UTC (rev 16227)
@@ -17,6 +17,7 @@
ElasticStrain1D.hh \
ElasticStress1D.hh \
ElasticIsotropic3D.hh \
+ ElasticIsotropic3D.icc \
ElasticMaterial.hh \
ElasticMaterial.icc \
ElasticPlaneStrain.hh \
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.cc 2010-02-04 00:30:58 UTC (rev 16227)
@@ -178,6 +178,59 @@
} // reformResidual
// ----------------------------------------------------------------------
+// Reform system residual.
+void
+pylith::problems::Formulation::reformResidualLumped(const PetscVec* tmpResidualVec,
+ const PetscVec* tmpSolutionVec)
+{ // reformResidualLumped
+ assert(0 != _fields);
+
+ // Update section view of field.
+ if (0 != tmpSolutionVec) {
+ topology::Field<topology::Mesh>& solution = _fields->solution();
+ solution.scatterVectorToSection(*tmpSolutionVec);
+ } // if
+
+ // Set residual to zero.
+ topology::Field<topology::Mesh>& residual = _fields->get("residual");
+ residual.zero();
+
+ // Add in contributions that require assembly.
+ int numIntegrators = _meshIntegrators.size();
+ assert(numIntegrators > 0); // must have at least 1 bulk integrator
+ for (int i=0; i < numIntegrators; ++i) {
+ _meshIntegrators[i]->timeStep(_dt);
+ _meshIntegrators[i]->integrateResidualLumped(residual, _t, _fields);
+ } // for
+ numIntegrators = _submeshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i) {
+ _submeshIntegrators[i]->timeStep(_dt);
+ _submeshIntegrators[i]->integrateResidualLumped(residual, _t, _fields);
+ } // for
+
+ // Assemble residual.
+ residual.complete();
+
+ // Add in contributions that do not require assembly.
+ numIntegrators = _meshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _meshIntegrators[i]->integrateResidualLumpedAssembled(residual, _t, _fields);
+ numIntegrators = _submeshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _submeshIntegrators[i]->integrateResidualLumpedAssembled(residual, _t, _fields);
+
+ // Update PETSc view of residual
+ if (0 != tmpResidualVec)
+ residual.scatterSectionToVector(*tmpResidualVec);
+ else
+ residual.scatterSectionToVector();
+
+ // TODO: Move this to SolverLinear
+ if (0 != tmpResidualVec)
+ VecScale(*tmpResidualVec, -1.0);
+} // reformResidualLumped
+
+// ----------------------------------------------------------------------
// Reform system Jacobian.
void
pylith::problems::Formulation::reformJacobian(const PetscVec* tmpSolutionVec)
Modified: short/3D/PyLith/trunk/libsrc/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/libsrc/problems/Formulation.hh 2010-02-04 00:30:58 UTC (rev 16227)
@@ -111,6 +111,14 @@
void reformResidual(const PetscVec* tmpResidualVec =0,
const PetscVec* tmpSolutionVec =0);
+ /** Reform system residual.
+ *
+ * @param tmpResidualVec Temporary PETSc vector for residual.
+ * @param tmpSolutionVec Temporary PETSc vector for solution.
+ */
+ void reformResidualLumped(const PetscVec* tmpResidualVec =0,
+ const PetscVec* tmpSolutionVec =0);
+
/* Reform system Jacobian.
*
* @param tmpSolveSolnVec Temporary PETSc vector for solution.
Modified: short/3D/PyLith/trunk/modulesrc/problems/Formulation.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/problems/Formulation.i 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/modulesrc/problems/Formulation.i 2010-02-04 00:30:58 UTC (rev 16227)
@@ -85,6 +85,14 @@
void reformResidual(const PetscVec* tmpResidualVec =0,
const PetscVec* tmpSolveSolnVec =0);
+ /** Reform system residual for case with lumped Jacobian.
+ *
+ * @param tmpResidualVec Temporary PETSc vector for residual.
+ * @param tmpSolveSolnVec Temporary PETSc vector for solution.
+ */
+ void reformResidualLumped(const PetscVec* tmpResidualVec =0,
+ const PetscVec* tmpSolveSolnVec =0);
+
/* Reform system Jacobian.
*
* @param tmpSolveSolnVec Temporary PETSc vector for solution.
Modified: short/3D/PyLith/trunk/pylith/problems/ExplicitLumped.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/ExplicitLumped.py 2010-02-03 22:08:47 UTC (rev 16226)
+++ short/3D/PyLith/trunk/pylith/problems/ExplicitLumped.py 2010-02-04 00:30:58 UTC (rev 16227)
@@ -215,6 +215,21 @@
return
+ def _reformResidual(self, t, dt):
+ """
+ Reform residual vector for operator.
+ """
+ self._info.log("Integrating residual term in operator.")
+ self._eventLogger.stagePush("Reform Residual")
+
+ self.updateSettings(self.jacobian, self.fields, t, dt)
+ self.reformResidualLumped()
+
+ self._eventLogger.stagePop()
+ self._debug.log(resourceUsageString())
+ return
+
+
def _reformJacobian(self, t, dt):
"""
Reform Jacobian matrix for operator.
More information about the CIG-COMMITS
mailing list