[cig-commits] r15880 - in short/3D/PyLith/branches/pylith-friction: . libsrc/bc libsrc/feassemble
brad at geodynamics.org
brad at geodynamics.org
Mon Oct 26 15:56:12 PDT 2009
Author: brad
Date: 2009-10-26 15:56:12 -0700 (Mon, 26 Oct 2009)
New Revision: 15880
Modified:
short/3D/PyLith/branches/pylith-friction/TODO
short/3D/PyLith/branches/pylith-friction/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/branches/pylith-friction/libsrc/bc/AbsorbingDampers.hh
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicit.hh
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.cc
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.hh
Log:
Added routine to lump cell matrix. Added integrateJacobian() for lumped Jacobian to ElasticityExplicit and AbsorbingDampers.
Modified: short/3D/PyLith/branches/pylith-friction/TODO
===================================================================
--- short/3D/PyLith/branches/pylith-friction/TODO 2009-10-26 21:21:11 UTC (rev 15879)
+++ short/3D/PyLith/branches/pylith-friction/TODO 2009-10-26 22:56:12 UTC (rev 15880)
@@ -5,6 +5,13 @@
Data structure to hold fault Lagrange vertices and conventional vertices
array = numVertics x 3 (i, j, k)
+Status?
+
+ Field split
+ uniform global refinement for tets with faults
+ customized SNES line search
+
+
----------------------------------------------------------------------
FRICTION
----------------------------------------------------------------------
@@ -29,6 +36,10 @@
----------------------------------------------------------------------
Unit tests
+
+ Integrator
+ _lumpCellMatrix
+
FaultCohesiveKin
integrateJacobian (lumped)
adjustSolnLumped
@@ -44,16 +55,6 @@
integrateJacobian (lumped)
-FaultCohesiveKin (C++)
- integrateJacobian (lumped)
- adjustSolnLumped()
-
-ElasticityExplicit
- integrateJacobian (lumped)
-
-AbsorbingDampers
- integrateJacobian (lumped)
-
FaultCohesiveDynL (C++)
integrateJacobian (lumped)
adjustSolnLumped()
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/bc/AbsorbingDampers.cc 2009-10-26 21:21:11 UTC (rev 15879)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/bc/AbsorbingDampers.cc 2009-10-26 22:56:12 UTC (rev 15880)
@@ -480,6 +480,114 @@
} // integrateJacobian
// ----------------------------------------------------------------------
+// Integrate contributions to Jacobian matrix (A) associated with
+void
+pylith::bc::AbsorbingDampers::integrateJacobian(
+ const topology::Field<topology::Mesh>& jacobian,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateJacobian
+ assert(0 != _quadrature);
+ assert(0 != _boundaryMesh);
+ assert(0 != fields);
+
+ // 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();
+
+ // Get cell information
+ const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
+ assert(!sieveSubMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ sieveSubMesh->heightStratum(1);
+ assert(!cells.isNull());
+ const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ assert(dt > 0);
+
+ // Get sections
+ const ALE::Obj<SubRealSection>& dampersSection =
+ _parameters->get("damping constants").section();
+ assert(!dampersSection.isNull());
+
+ const topology::Field<topology::Mesh>& solution = fields->solution();
+ const ALE::Obj<SieveMesh>& sieveMesh = solution.mesh().sieveMesh();
+ const ALE::Obj<RealSection>& solutionSection = solution.section();
+ assert(!solutionSection.isNull());
+
+ const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
+ assert(!jacobianSection.isNull());
+ topology::Mesh::UpdateAddVisitor jacobianVisitor(*jacobianSection,
+ &_cellVector[0]);
+
+#if !defined(PRECOMPUTE_GEOMETRY)
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveSubMesh->getRealSection("coordinates");
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
+#endif
+
+ // Allocate matrix for cell values.
+ _initCellMatrix();
+
+ for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
+
+ // Reset element vector to zero
+ _resetCellMatrix();
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ assert(numQuadPts*spaceDim == dampersSection->getFiberDimension(*c_iter));
+ const double* dampingConstsCell = dampersSection->restrictPoint(*c_iter);
+
+ // Compute Jacobian for absorbing bc terms
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt =
+ quadWts[iQuad] * jacobianDet[iQuad] / (2.0 * dt);
+ 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 * dampingConstsCell[iQuad*spaceDim+iDim];
+ } // for
+ } // for
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+2*spaceDim))));
+ _lumpCellMatrix();
+
+ // Assemble cell contribution into lumped matrix.
+ jacobianVisitor.clear();
+ sieveMesh->updateAdd(*c_iter, jacobianVisitor);
+ } // for
+
+ _needNewJacobian = false;
+} // integrateJacobian
+
+// ----------------------------------------------------------------------
// Verify configuration is acceptable.
void
pylith::bc::AbsorbingDampers::verifyConfiguration(const topology::Mesh& mesh) const
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/bc/AbsorbingDampers.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/bc/AbsorbingDampers.hh 2009-10-26 21:21:11 UTC (rev 15879)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/bc/AbsorbingDampers.hh 2009-10-26 22:56:12 UTC (rev 15880)
@@ -107,6 +107,17 @@
const double t,
topology::SolutionFields* const fields);
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param jacobian Jacobian of system.
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateJacobian(const topology::Field<topology::Mesh>& jacobian,
+ const double t,
+ topology::SolutionFields* const fields);
+
/** Verify configuration is acceptable.
*
* @param mesh Finite-element mesh
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicit.cc 2009-10-26 21:21:11 UTC (rev 15879)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicit.cc 2009-10-26 22:56:12 UTC (rev 15880)
@@ -431,5 +431,137 @@
_material->resetNeedNewJacobian();
} // integrateJacobian
+// ----------------------------------------------------------------------
+// Compute matrix associated with operator.
+void
+pylith::feassemble::ElasticityExplicit::integrateJacobian(
+ const topology::Field<topology::Mesh>& jacobian,
+ const double t,
+ topology::SolutionFields* fields)
+{ // integrateJacobian
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ 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 parameters used in integration.
+ const double dt = _dt;
+ const double dt2 = dt*dt;
+ assert(dt > 0);
+
+ // Get sections
+ double_array dispTCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& dispTSection =
+ fields->get("disp(t)").section();
+ assert(!dispTSection.isNull());
+
+ 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);
+
+ // Loop over cells
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ _logger->eventBegin(geometryEvent);
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
+ _logger->eventEnd(geometryEvent);
+
+ // Get state variables for cell.
+ _logger->eventBegin(stateVarsEvent);
+ _material->retrievePropsAndVars(*c_iter);
+ _logger->eventEnd(stateVarsEvent);
+
+ // Reset element matrix to zero
+ _resetCellMatrix();
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ // Get material physical properties at quadrature points for this cell
+ 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
+ } // for
+ PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+spaceDim))));
+ _lumpCellMatrix();
+ _logger->eventEnd(computeEvent);
+
+ // Assemble cell contribution into lumped matrix.
+ _logger->eventBegin(updateEvent);
+ jacobianVisitor.clear();
+ sieveMesh->updateAdd(*c_iter, jacobianVisitor);
+ _logger->eventEnd(updateEvent);
+ } // for
+
+ _needNewJacobian = false;
+ _material->resetNeedNewJacobian();
+} // integrateJacobian
+
+
// End of file
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicit.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicit.hh 2009-10-26 21:21:11 UTC (rev 15879)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicit.hh 2009-10-26 22:56:12 UTC (rev 15880)
@@ -112,6 +112,17 @@
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(const topology::Field<topology::Mesh>& jacobian,
+ const double t,
+ topology::SolutionFields* const fields);
+
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.cc 2009-10-26 21:21:11 UTC (rev 15879)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.cc 2009-10-26 22:56:12 UTC (rev 15880)
@@ -153,5 +153,29 @@
_cellMatrix = 0.0;
} // _resetCellMatrix
+// ----------------------------------------------------------------------
+// Lump cell matrix, putting the result in the cell vector using
+// equivalent forces for rigid body motion.
+template<typename quadrature_type>
+void
+pylith::feassemble::Integrator<quadrature_type>::_lumpCellMatrix(void)
+{ // _lumpCellMatrix
+ assert(0 != _quadrature);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ _cellVector = 0.0;
+ double value = 0.0;
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ value = 0.0;
+ for (int jBasis=0, indexI=iBasis*spaceDim; jBasis < numBasis; ++jBasis)
+ value += _cellMatrix[indexI+jBasis*spaceDim+iDim];
+ _cellVector[iBasis*spaceDim+iDim] = value;
+ } // for
+
+ PetscLogFlops(numBasis*numBasis*spaceDim);
+} // _lumpCellMatrix
+
+
// End of file
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.hh 2009-10-26 21:21:11 UTC (rev 15879)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.hh 2009-10-26 22:56:12 UTC (rev 15880)
@@ -256,6 +256,10 @@
/// Zero out matrix containing result of integration for cell.
void _resetCellMatrix(void);
+ /// Lump cell matrix, putting the result in the cell vector using
+ /// equivalent forces for rigid body motion.
+ void _lumpCellMatrix(void);
+
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
More information about the CIG-COMMITS
mailing list