[cig-commits] r15957 - in short/3D/PyLith/branches/pylith-friction/libsrc: . feassemble
brad at geodynamics.org
brad at geodynamics.org
Wed Nov 11 18:06:05 PST 2009
Author: brad
Date: 2009-11-11 18:06:05 -0800 (Wed, 11 Nov 2009)
New Revision: 15957
Added:
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicitLgDeform.cc
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicitLgDeform.hh
Modified:
short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Makefile.am
Log:
Added explicit elasticity integrator for large deformations.
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am 2009-11-12 01:26:26 UTC (rev 15956)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am 2009-11-12 02:06:05 UTC (rev 15957)
@@ -75,6 +75,7 @@
feassemble/ElasticityExplicit.cc \
feassemble/IntegratorElasticityLgDeform.cc \
feassemble/ElasticityImplicitLgDeform.cc \
+ feassemble/ElasticityExplicitLgDeform.cc \
materials/Metadata.cc \
materials/Material.cc \
materials/ElasticMaterial.cc \
Added: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicitLgDeform.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicitLgDeform.cc 2009-11-12 02:06:05 UTC (rev 15957)
@@ -0,0 +1,604 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "ElasticityExplicitLgDeform.hh" // implementation of class methods
+
+#include "Quadrature.hh" // USES Quadrature
+#include "CellGeometry.hh" // USES CellGeometry
+
+#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 PRECOMPUTE_GEOMETRY
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::ElasticityExplicitLgDeform::ElasticityExplicitLgDeform(void) :
+ _dtm1(-1.0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::ElasticityExplicitLgDeform::~ElasticityExplicitLgDeform(void)
+{ // destructor
+ deallocate();
+} // destructor
+
+// ----------------------------------------------------------------------
+// Deallocate PETSc and local data structures.
+void
+pylith::feassemble::ElasticityExplicitLgDeform::deallocate(void)
+{ // deallocate
+ IntegratorElasticityLgDeform::deallocate();
+} // deallocate
+
+// ----------------------------------------------------------------------
+// Set time step for advancing from time t to time t+dt.
+void
+pylith::feassemble::ElasticityExplicitLgDeform::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::ElasticityExplicitLgDeform::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::ElasticityExplicitLgDeform::integrateResidual(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateResidual
+ /// Member prototype for _elasticityResidualXD()
+ typedef void (pylith::feassemble::ElasticityExplicitLgDeform::*elasticityResidual_fn_type)
+ (const double_array&, 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& vertices = _quadrature->refGeometry().vertices();
+ 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::ElasticityExplicitLgDeform::_elasticityResidual1D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain1D;
+ } else if (2 == cellDim) {
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityExplicitLgDeform::_elasticityResidual2D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain2D;
+ } else if (3 == cellDim) {
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityExplicitLgDeform::_elasticityResidual3D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain3D;
+ } else
+ assert(0);
+
+ // Allocate vectors for cell values.
+ double_array dispTCell(numBasis*spaceDim);
+ double_array dispTmdtCell(numBasis*spaceDim);
+ double_array deformCell(numQuadPts*spaceDim*spaceDim);
+ double_array strainCell(numQuadPts*tensorSize);
+ 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);
+
+ // 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 vector to zero
+ _resetCellVector();
+
+ // 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();
+
+ // Restrict input fields to cell
+ _logger->eventBegin(restrictEvent);
+ dispTVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+ dispTmdtVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTmdtVisitor);
+ _logger->eventEnd(restrictEvent);
+
+ // 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);
+
+ // 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)));
+ _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
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(6*spaceDim))));
+ _logger->eventEnd(computeEvent);
+
+ // Compute B(transpose) * sigma, first computing strains
+ _logger->eventBegin(stressEvent);
+ _calcDeformation(&deformCell, basisDeriv, vertices, dispTCell,
+ numBasis, numQuadPts, spaceDim);
+ calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
+ const double_array& stressCell = _material->calcStress(strainCell, true);
+ _logger->eventEnd(stressEvent);
+
+ _logger->eventBegin(computeEvent);
+ CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell, dispTCell);
+ _logger->eventEnd(computeEvent);
+
+ // Assemble cell contribution into field
+ _logger->eventBegin(updateEvent);
+ residualVisitor.clear();
+ sieveMesh->updateClosure(*c_iter, residualVisitor);
+ _logger->eventEnd(updateEvent);
+ } // for
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Compute matrix associated with operator.
+void
+pylith::feassemble::ElasticityExplicitLgDeform::integrateJacobian(
+ topology::Jacobian* jacobian,
+ const double t,
+ topology::SolutionFields* fields)
+{ // integrateJacobian
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(0 != jacobian);
+ assert(0 != fields);
+
+ const int setupEvent = _logger->eventId("ElIJ setup");
+ const int geometryEvent = _logger->eventId("ElIJ geometry");
+ const int computeEvent = _logger->eventId("ElIJ compute");
+ const int restrictEvent = _logger->eventId("ElIJ restrict");
+ const int stateVarsEvent = _logger->eventId("ElIJ stateVars");
+ const int updateEvent = _logger->eventId("ElIJ update");
+
+ _logger->eventBegin(setupEvent);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const int tensorSize = _material->tensorSize();
+ if (cellDim != spaceDim)
+ throw std::logic_error("Don't know how to integrate elasticity " \
+ "contribution to Jacobian matrix for cells with " \
+ "different dimensions than the spatial dimension.");
+
+ // Allocate vectors for cell data.
+ double_array dispTCell(numBasis*spaceDim);
+
+ // Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const int materialId = _material->id();
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<RealSection>& dispTSection =
+ fields->get("disp(t)").section();
+ assert(!dispTSection.isNull());
+
+ // Get sparse matrix
+ const PetscMat jacobianMat = jacobian->matrix();
+ assert(0 != jacobianMat);
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ const double dt2 = dt*dt;
+ assert(dt > 0);
+
+ const ALE::Obj<SieveMesh::order_type>& globalOrder =
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", dispTSection);
+ assert(!globalOrder.isNull());
+ // We would need to request unique points here if we had an interpolated mesh
+ topology::Mesh::IndicesVisitor jacobianVisitor(*dispTSection, *globalOrder,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
+
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
+ _logger->eventEnd(setupEvent);
+
+ // 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))));
+ _logger->eventEnd(computeEvent);
+
+ // 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.");
+ _logger->eventEnd(updateEvent);
+ } // for
+
+ _needNewJacobian = false;
+ _material->resetNeedNewJacobian();
+} // integrateJacobian
+
+// ----------------------------------------------------------------------
+// Compute matrix associated with operator.
+void
+pylith::feassemble::ElasticityExplicitLgDeform::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->updateClosure(*c_iter, jacobianVisitor);
+ _logger->eventEnd(updateEvent);
+ } // for
+
+ _needNewJacobian = false;
+ _material->resetNeedNewJacobian();
+} // integrateJacobian
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicitLgDeform.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicitLgDeform.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityExplicitLgDeform.hh 2009-11-12 02:06:05 UTC (rev 15957)
@@ -0,0 +1,147 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/feassemble/ElasticityExplicitLgDeform.hh
+ *
+ * @brief Explicit time integration of dynamic elasticity equation
+ * with large rigid body motion and small strains using
+ * finite-elements.
+ */
+
+#if !defined(pylith_feassemble_elasticityexplicitlgdeform_hh)
+#define pylith_feassemble_elasticityexplicitlgdeform_hh
+
+// Include directives ---------------------------------------------------
+#include "IntegratorElasticityLgDeform.hh" // ISA IntegratorElasticityLgDeform
+
+// ElasticityExplicitLgDeform -------------------------------------------
+/**@brief Explicit time integration of the dynamic elasticity equation
+ * with large rigid body motion and small strains using 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::ElasticityExplicitLgDeform :
+ public IntegratorElasticityLgDeform
+{ // ElasticityExplicitLgDeform
+ friend class TestElasticityExplicitLgDeform; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ ElasticityExplicitLgDeform(void);
+
+ /// Destructor
+ ~ElasticityExplicitLgDeform(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 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(const topology::Field<topology::Mesh>& jacobian,
+ const double t,
+ topology::SolutionFields* const fields);
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ /// Not implemented.
+ ElasticityExplicitLgDeform(const ElasticityExplicitLgDeform&);
+
+ /// Not implemented
+ const ElasticityExplicitLgDeform& operator=(const ElasticityExplicitLgDeform&);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ double _dtm1; ///< Time step for t-dt1 -> t
+
+}; // ElasticityExplicitLgDeform
+
+#endif // pylith_feassemble_elasticityexplicitlgdeform_hh
+
+
+// End of file
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Makefile.am 2009-11-12 01:26:26 UTC (rev 15956)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Makefile.am 2009-11-12 02:06:05 UTC (rev 15957)
@@ -19,6 +19,7 @@
Constraint.hh \
Constraint.icc \
ElasticityExplicit.hh \
+ ElasticityExplicitLgDeform.hh \
ElasticityImplicit.hh \
ElasticityImplicitLgDeform.hh \
Integrator.hh \
More information about the CIG-COMMITS
mailing list