[cig-commits] r15943 - in short/3D/PyLith/branches/pylith-friction/libsrc: . feassemble
brad at geodynamics.org
brad at geodynamics.org
Mon Nov 9 16:12:12 PST 2009
Author: brad
Date: 2009-11-09 16:12:11 -0800 (Mon, 09 Nov 2009)
New Revision: 15943
Added:
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicitLgDeform.cc
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicitLgDeform.hh
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.hh
Modified:
short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Makefile.am
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/feassemblefwd.hh
Log:
Started implementing large deformations.
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am 2009-11-09 21:27:57 UTC (rev 15942)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am 2009-11-10 00:12:11 UTC (rev 15943)
@@ -71,6 +71,7 @@
feassemble/Quadrature2Din3D.cc \
feassemble/Quadrature3D.cc \
feassemble/IntegratorElasticity.cc \
+ feassemble/IntegratorElasticityLgDeform.cc \
feassemble/ElasticityImplicit.cc \
feassemble/ElasticityExplicit.cc \
materials/Metadata.cc \
Added: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicitLgDeform.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicitLgDeform.cc 2009-11-10 00:12:11 UTC (rev 15943)
@@ -0,0 +1,532 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "ElasticityImplicit.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/EventLogger.hh" // USES EventLogger
+#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 "spatialdata/spatialdb/GravityField.hh" // USES GravityField
+
+#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::ElasticityImplicit::ElasticityImplicit(void) :
+ _dtm1(-1.0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::ElasticityImplicit::~ElasticityImplicit(void)
+{ // destructor
+ deallocate();
+} // destructor
+
+// ----------------------------------------------------------------------
+// Deallocate PETSc and local data structures.
+void
+pylith::feassemble::ElasticityImplicit::deallocate(void)
+{ // deallocate
+ IntegratorElasticity::deallocate();
+} // deallocate
+
+// ----------------------------------------------------------------------
+// Set time step for advancing from time t to time t+dt.
+void
+pylith::feassemble::ElasticityImplicit::timeStep(const double dt)
+{ // timeStep
+ if (_dt != -1.0)
+ _dtm1 = _dt;
+ else
+ _dtm1 = dt;
+ _dt = dt;
+ if (0 != _material)
+ _material->timeStep(_dt);
+} // timeStep
+
+// ----------------------------------------------------------------------
+// Get stable time step for advancing from time t to time t+dt.
+double
+pylith::feassemble::ElasticityImplicit::stableTimeStep(const topology::Mesh& mesh) const
+{ // stableTimeStep
+ assert(0 != _material);
+ return _material->stableTimeStepImplicit(mesh);
+} // stableTimeStep
+
+// ----------------------------------------------------------------------
+// Set flag for setting constraints for total field solution or
+// incremental field solution.
+void
+pylith::feassemble::ElasticityImplicit::useSolnIncr(const bool flag)
+{ // useSolnIncr
+ assert(0 != _material);
+ _useSolnIncr = flag;
+ _material->useElasticBehavior(!_useSolnIncr);
+} // useSolnIncr
+
+// ----------------------------------------------------------------------
+// Integrate constributions to residual term (r) for operator.
+void
+pylith::feassemble::ElasticityImplicit::integrateResidual(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateResidual
+ /// Member prototype for _elasticityResidualXD()
+ typedef void (pylith::feassemble::ElasticityImplicit::*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();
+ 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::ElasticityImplicit::_elasticityResidual1D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
+ } else if (2 == cellDim) {
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityImplicit::_elasticityResidual2D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
+ } else if (3 == cellDim) {
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityImplicit::_elasticityResidual3D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D;
+ } else
+ assert(0);
+
+ // Allocate vectors for cell values.
+ double_array dispTCell(numBasis*spaceDim);
+ double_array dispTIncrCell(numBasis*spaceDim);
+ double_array dispTpdtCell(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>& dispTIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispTIncrSection.isNull());
+ topology::Mesh::RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
+ numBasis*spaceDim,
+ &dispTIncrCell[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());
+
+ _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();
+
+ // Restrict input fields to cell
+ _logger->eventBegin(restrictEvent);
+ dispTVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+ dispTIncrVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+ _logger->eventBegin(restrictEvent);
+
+ // 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 current estimate of displacement at time t+dt using
+ // solution increment.
+ dispTpdtCell = dispTCell + dispTIncrCell;
+
+ // 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 B(transpose) * sigma, first computing strains
+ _logger->eventBegin(stressEvent);
+ calcTotalStrainFn(&strainCell, basisDeriv, dispTpdtCell,
+ numBasis, numQuadPts);
+ const double_array& stressCell = _material->calcStress(strainCell, true);
+ _logger->eventEnd(stressEvent);
+
+ _logger->eventBegin(computeEvent);
+ CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell);
+ _logger->eventEnd(computeEvent);
+
+#if 0 // DEBUGGING
+ std::cout << "Updating residual for cell " << *c_iter << std::endl;
+ for(int i = 0; i < _quadrature->spaceDim() * _quadrature->numBasis(); ++i) {
+ std::cout << " v["<<i<<"]: " << _cellVector[i] << std::endl;
+ }
+#endif
+ // Assemble cell contribution into field
+ _logger->eventBegin(updateEvent);
+ residualVisitor.clear();
+ sieveMesh->updateClosure(*c_iter, residualVisitor);
+ _logger->eventEnd(updateEvent);
+ } // for
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Compute stiffness matrix.
+void
+pylith::feassemble::ElasticityImplicit::integrateJacobian(
+ topology::Jacobian* jacobian,
+ const double t,
+ topology::SolutionFields* fields)
+{ // integrateJacobian
+ /// Member prototype for _elasticityJacobianXD()
+ typedef void (pylith::feassemble::ElasticityImplicit::*elasticityJacobian_fn_type)
+ (const double_array&);
+
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(0 != _logger);
+ 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.");
+
+ // Set variables dependent on dimension of cell
+ totalStrain_fn_type calcTotalStrainFn;
+ elasticityJacobian_fn_type elasticityJacobianFn;
+ if (1 == cellDim) {
+ elasticityJacobianFn =
+ &pylith::feassemble::ElasticityImplicit::_elasticityJacobian1D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
+ } else if (2 == cellDim) {
+ elasticityJacobianFn =
+ &pylith::feassemble::ElasticityImplicit::_elasticityJacobian2D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
+ } else if (3 == cellDim) {
+ elasticityJacobianFn =
+ &pylith::feassemble::ElasticityImplicit::_elasticityJacobian3D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D;
+ } else
+ assert(0);
+
+ // Allocate vector for total strain
+ double_array dispTCell(numBasis*spaceDim);
+ double_array dispTIncrCell(numBasis*spaceDim);
+ double_array dispTpdtCell(numBasis*spaceDim);
+ double_array strainCell(numQuadPts*tensorSize);
+ strainCell = 0.0;
+
+ // 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>& dispTIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispTIncrSection.isNull());
+ topology::Mesh::RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
+ numBasis*spaceDim,
+ &dispTIncrCell[0]);
+
+ // Get sparse matrix
+ const PetscMat jacobianMat = jacobian->matrix();
+ assert(0 != jacobianMat);
+
+ // Get parameters used in integration.
+ const double 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();
+
+ // Restrict input fields to cell
+ _logger->eventBegin(restrictEvent);
+ dispTVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+ dispTIncrVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+ _logger->eventBegin(restrictEvent);
+
+ // 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();
+
+ // Compute current estimate of displacement at time t+dt using
+ // solution increment.
+ dispTpdtCell = dispTCell + dispTIncrCell;
+
+ _logger->eventBegin(computeEvent);
+ // Compute strains
+ calcTotalStrainFn(&strainCell, basisDeriv, dispTpdtCell,
+ numBasis, numQuadPts);
+
+ // Get "elasticity" matrix at quadrature points for this cell
+ const double_array& elasticConsts =
+ _material->calcDerivElastic(strainCell);
+
+ CALL_MEMBER_FN(*this, elasticityJacobianFn)(elasticConsts);
+ _logger->eventEnd(computeEvent);
+
+ if (_quadrature->checkConditioning()) {
+ int n = numBasis*spaceDim;
+ int lwork = 5*n;
+ int idummy = 0;
+ int lierr = 0;
+ double *elemMat = new double[n*n];
+ double *svalues = new double[n];
+ double *work = new double[lwork];
+ double minSV = 0;
+ double maxSV = 0;
+ double sdummy = 0;
+
+ const int n2 = n*n;
+ for (int i = 0; i < n2; ++i)
+ elemMat[i] = _cellMatrix[i];
+ lapack_dgesvd("N", "N", &n, &n, elemMat, &n, svalues,
+ &sdummy, &idummy, &sdummy, &idummy, work,
+ &lwork, &lierr);
+ if (lierr)
+ throw std::runtime_error("Lapack SVD failed");
+ minSV = svalues[n-7];
+ maxSV = svalues[0];
+ std::cout << "Element " << *c_iter << std::endl;
+ for(int i = 0; i < n; ++i)
+ std::cout << " sV["<<i<<"] = " << svalues[i] << std::endl;
+ std::cout << " kappa(elemMat) = " << maxSV/minSV << std::endl;
+ delete [] elemMat;
+ delete [] svalues;
+ delete [] work;
+ } // if
+
+ // 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
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicitLgDeform.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicitLgDeform.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/ElasticityImplicitLgDeform.hh 2009-11-10 00:12:11 UTC (rev 15943)
@@ -0,0 +1,140 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/feassemble/ElasticityImplicit.hh
+ *
+ * @brief Implicit time integration of quasi-static elasticity equation
+ * 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 = 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(t+dt) is a vector that depends on the
+ * field at time t+dt and t, and du is zero at unknown DOF and set to
+ * the constrained values at known DOF.
+ *
+ * Contributions to the RHS (b) include body forces, which are either
+ * independent of u (small strain case) or are computed based on the
+ * displacements at time t. The RHS also includes the internal force
+ * vector, which is either constant for a time step (small strain,
+ * elastic rheology) or changes with each iteration (large strain or
+ * non-elastic rheology). The internal force vector is subtracted from the
+ * existing force vector to get the residual. This object also computes
+ * the entire stiffness matrix (A).
+ *
+ * See governing equations section of user manual for more
+ * information.
+ */
+
+#if !defined(pylith_feassemble_elasticityimplicit_hh)
+#define pylith_feassemble_elasticityimplicit_hh
+
+// Include directives ---------------------------------------------------
+#include "IntegratorElasticity.hh" // ISA IntegratorElasticity
+
+// ElasticityImplicit ---------------------------------------------------
+class pylith::feassemble::ElasticityImplicit : public IntegratorElasticity
+{ // ElasticityImplicit
+ friend class TestElasticityImplicit; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ ElasticityImplicit(void);
+
+ /// Destructor
+ ~ElasticityImplicit(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);
+
+ /** Get stable time step for advancing from time t to time t+dt.
+ *
+ * Default is current time step.
+ *
+ * @param mesh Finite-element mesh.
+ * @returns Time step
+ */
+ double stableTimeStep(const topology::Mesh& mesh) const;
+
+ /** 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 residual part of RHS for 3-D finite elements.
+ * Includes gravity and element internal force contribution.
+ *
+ * We assume that the effects of boundary conditions are already
+ * included in the residual (tractions, concentrated nodal forces,
+ * and contributions to internal force vector due to
+ * displacement/velocity BC). This routine computes the additional
+ * external loads due to body forces plus the
+ * element internal forces for the current stress state.
+ *
+ * @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);
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ /// Not implemented.
+ ElasticityImplicit(const ElasticityImplicit&);
+
+ /// Not implemented
+ const ElasticityImplicit& operator=(const ElasticityImplicit&);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ double _dtm1; ///< Time step for t-dt1 -> t
+
+}; // ElasticityImplicit
+
+#endif // pylith_feassemble_elasticityimplicit_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc 2009-11-10 00:12:11 UTC (rev 15943)
@@ -0,0 +1,751 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "IntegratorElasticityLgDeform.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/Fields.hh" // USES Fields
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+#include "pylith/utils/EventLogger.hh" // USES EventLogger
+
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
+#include "pylith/utils/array.hh" // USES double_array
+
+#include <cstring> // USES memcpy()
+#include <strings.h> // USES strcasecmp()
+#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::IntegratorElasticityLgDeform::IntegratorElasticityLgDeform(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::IntegratorElasticityLgDeform::~IntegratorElasticityLgDeform(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Determine whether we need to recompute the Jacobian.
+bool
+pylith::feassemble::IntegratorElasticityLgDeform::needNewJacobian(void)
+{ // needNewJacobian
+ assert(0 != _material);
+ if (!_needNewJacobian)
+ _needNewJacobian = _material->needNewJacobian();
+ return _needNewJacobian;
+} // needNewJacobian
+
+// ----------------------------------------------------------------------
+// Update state variables as needed.
+void
+pylith::feassemble::IntegratorElasticityLgDeform::updateStateVars(
+ const double t,
+ topology::SolutionFields* const fields)
+{ // updateStateVars
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(0 != fields);
+
+ // No need to update state vars if material doesn't have any.
+ if (!_material->hasStateVars())
+ return;
+
+ // Get cell information that doesn't depend on particular cell
+ const int cellDim = _quadrature->cellDim();
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int tensorSize = _material->tensorSize();
+ totalStrain_fn_type calcTotalStrainFn;
+ if (1 == cellDim) {
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain1D;
+ } else if (2 == cellDim) {
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain2D;
+ } else if (3 == cellDim) {
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain3D;
+ } else {
+ std::cerr << "Bad cell dimension '" << cellDim << "'." << std::endl;
+ assert(0);
+ throw std::logic_error("Bad cell dimension in IntegratorElasticityLgDeform.");
+ } // else
+
+ // Allocate arrays for cell data.
+ double_array dispCell(numBasis*spaceDim);
+ double_array strainCell(numQuadPts*tensorSize);
+ strainCell = 0.0;
+
+ // 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 fields
+ const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
+ const ALE::Obj<RealSection>& dispSection = disp.section();
+ assert(!dispSection.isNull());
+ topology::Mesh::RestrictVisitor dispVisitor(*dispSection,
+ dispCell.size(), &dispCell[0]);
+
+#if !defined(PRECOMPUTE_GEOMETRY)
+ 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]);
+#endif
+
+ // Loop over cells
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Retrieve geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
+
+ // Restrict input fields to cell
+ dispVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispVisitor);
+
+ // Get cell geometry information that depends on cell
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ // Compute strains
+ calcTotalStrainFn(&strainCell, basisDeriv, dispCell,
+ numBasis, numQuadPts);
+
+ // Update material state
+ _material->updateStateVars(strainCell, *c_iter);
+ } // for
+} // updateStateVars
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_calcStrainStressField(
+ topology::Field<topology::Mesh>* field,
+ const char* name,
+ topology::SolutionFields* const fields)
+{ // _calcStrainStressField
+ assert(0 != field);
+ assert(0 != _quadrature);
+ assert(0 != _material);
+
+ const bool calcStress = (0 == strcasecmp(name, "stress")) ? true : false;
+
+ // Get cell information that doesn't depend on particular cell
+ const int cellDim = _quadrature->cellDim();
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int tensorSize = _material->tensorSize();
+ totalStrain_fn_type calcTotalStrainFn;
+ if (1 == cellDim) {
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain1D;
+ } else if (2 == cellDim) {
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain2D;
+ } else if (3 == cellDim) {
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain3D;
+ } else {
+ std::cerr << "Bad cell dimension '" << cellDim << "'." << std::endl;
+ assert(0);
+ throw std::logic_error("Bad cell dimension in IntegratorElasticityLgDeform.");
+ } // else
+
+ // Allocate arrays for cell data.
+ double_array dispCell(numBasis*spaceDim);
+ double_array strainCell(numQuadPts*tensorSize);
+ strainCell = 0.0;
+ double_array stressCell(numQuadPts*tensorSize);
+ stressCell = 0.0;
+
+ // Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = field->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 field
+ const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
+ const ALE::Obj<RealSection>& dispSection = disp.section();
+ assert(!dispSection.isNull());
+ topology::Mesh::RestrictVisitor dispVisitor(*dispSection,
+ dispCell.size(), &dispCell[0]);
+
+#if !defined(PRECOMPUTE_GEOMETRY)
+ 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]);
+#endif
+
+ const ALE::Obj<RealSection>& fieldSection = field->section();
+ assert(!fieldSection.isNull());
+
+ // Loop over cells
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Retrieve geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
+
+ // Restrict input fields to cell
+ dispVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispVisitor);
+
+ // Get cell geometry information that depends on cell
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ // Compute strains
+ calcTotalStrainFn(&strainCell, basisDeriv, dispCell,
+ numBasis, numQuadPts);
+
+ if (!calcStress)
+ fieldSection->updatePoint(*c_iter, &strainCell[0]);
+ else {
+ _material->retrievePropsAndVars(*c_iter);
+ stressCell = _material->calcStress(strainCell);
+ fieldSection->updatePoint(*c_iter, &stressCell[0]);
+ } // else
+ } // for
+} // _calcStrainStressField
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_calcStressFromStrain(
+ topology::Field<topology::Mesh>* field)
+{ // _calcStressFromStrain
+ assert(0 != field);
+ assert(0 != _quadrature);
+ assert(0 != _material);
+
+ const int cellDim = _quadrature->cellDim();
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int tensorSize = _material->tensorSize();
+
+ // Allocate arrays for cell data.
+ double_array strainCell(numQuadPts*tensorSize);
+ strainCell = 0.0;
+ double_array stressCell(numQuadPts*tensorSize);
+ stressCell = 0.0;
+
+ // Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = field->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 field
+ const ALE::Obj<RealSection>& fieldSection = field->section();
+ assert(!fieldSection.isNull());
+
+ // Loop over cells
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ fieldSection->restrictPoint(*c_iter, &strainCell[0], strainCell.size());
+ _material->retrievePropsAndVars(*c_iter);
+ stressCell = _material->calcStress(strainCell);
+ fieldSection->updatePoint(*c_iter, &stressCell[0]);
+ } // for
+} // _calcStressFromStrain
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in residual for 1-D cells.
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_elasticityResidual1D(
+ const double_array& stress)
+{ // _elasticityResidual1D
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const double_array& quadWts = _quadrature->quadWts();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ assert(1 == cellDim);
+ assert(quadWts.size() == numQuadPts);
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double s11 = stress[iQuad];
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const double N1 = wt*basisDeriv[iQuad*numBasis+iBasis ];
+ _cellVector[iBasis*spaceDim ] -= N1*s11;
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(1+numBasis*5));
+} // _elasticityResidual1D
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in residual for 2-D cells.
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_elasticityResidual2D(
+ const double_array& stress)
+{ // _elasticityResidual2D
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const double_array& quadWts = _quadrature->quadWts();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ assert(2 == cellDim);
+ assert(quadWts.size() == numQuadPts);
+ const int stressSize = 3;
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double s11 = stress[iQuad*stressSize+0];
+ const double s22 = stress[iQuad*stressSize+1];
+ const double s12 = stress[iQuad*stressSize+2];
+ for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
+ iBasis < numBasis;
+ ++iBasis) {
+ const double N1 = wt*basisDeriv[iQ+iBasis*spaceDim ];
+ const double N2 = wt*basisDeriv[iQ+iBasis*spaceDim+1];
+ _cellVector[iBasis*spaceDim ] -= N1*s11 + N2*s12;
+ _cellVector[iBasis*spaceDim+1] -= N1*s12 + N2*s22;
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(1+numBasis*(8+2+9)));
+} // _elasticityResidual2D
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in residual for 3-D cells.
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_elasticityResidual3D(
+ const double_array& stress)
+{ // _elasticityResidual3D
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const double_array& quadWts = _quadrature->quadWts();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ assert(3 == cellDim);
+ assert(quadWts.size() == numQuadPts);
+ const int stressSize = 6;
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double s11 = stress[iQuad*stressSize+0];
+ const double s22 = stress[iQuad*stressSize+1];
+ const double s33 = stress[iQuad*stressSize+2];
+ const double s12 = stress[iQuad*stressSize+3];
+ const double s23 = stress[iQuad*stressSize+4];
+ const double s13 = stress[iQuad*stressSize+5];
+
+ for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
+ iBasis < numBasis;
+ ++iBasis) {
+ const int iBlock = iBasis*spaceDim;
+ const double N1 = wt*basisDeriv[iQ+iBasis*spaceDim+0];
+ const double N2 = wt*basisDeriv[iQ+iBasis*spaceDim+1];
+ const double N3 = wt*basisDeriv[iQ+iBasis*spaceDim+2];
+
+ _cellVector[iBlock ] -= N1*s11 + N2*s12 + N3*s13;
+ _cellVector[iBlock+1] -= N1*s12 + N2*s22 + N3*s23;
+ _cellVector[iBlock+2] -= N1*s13 + N2*s23 + N3*s33;
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(1+numBasis*(3+12)));
+} // _elasticityResidual3D
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in Jacobian for 1-D cells.
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_elasticityJacobian1D(
+ const double_array& elasticConsts)
+{ // _elasticityJacobian1D
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const double_array& quadWts = _quadrature->quadWts();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ assert(1 == cellDim);
+ assert(quadWts.size() == numQuadPts);
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double C1111 = elasticConsts[iQuad];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const double valI = wt*basisDeriv[iQ+iBasis]*C1111;
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double valIJ = valI * basisDeriv[iQ+jBasis];
+ const int iBlock = iBasis*spaceDim * (numBasis*spaceDim);
+ const int jBlock = jBasis*spaceDim;
+ _cellMatrix[iBlock+jBlock] += valIJ;
+ } // for
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*3)));
+} // _elasticityJacobian1D
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in Jacobian for 2-D cells.
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_elasticityJacobian2D(
+ const double_array& elasticConsts)
+{ // _elasticityJacobian2D
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const double_array& quadWts = _quadrature->quadWts();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ assert(2 == cellDim);
+ assert(quadWts.size() == numQuadPts);
+ const int numConsts = 6;
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ // tau_ij = C_ijkl * e_kl
+ // = C_ijlk * 0.5 (u_k,l + u_l,k)
+ // = 0.5 * C_ijkl * (u_k,l + u_l,k)
+ // divide C_ijkl by 2 if k != l
+ const double C1111 = elasticConsts[iQuad*numConsts+0];
+ const double C1122 = elasticConsts[iQuad*numConsts+1];
+ const double C1112 = elasticConsts[iQuad*numConsts+2]/2.0;
+ const double C2222 = elasticConsts[iQuad*numConsts+3];
+ const double C2212 = elasticConsts[iQuad*numConsts+4]/2.0;
+ const double C1212 = elasticConsts[iQuad*numConsts+5]/2.0;
+ for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
+ iBasis < numBasis;
+ ++iBasis) {
+ const double Nip = wt*basisDeriv[iQ+iBasis*spaceDim ];
+ const double Niq = wt*basisDeriv[iQ+iBasis*spaceDim+1];
+ const int iBlock = (iBasis*spaceDim ) * (numBasis*spaceDim);
+ const int iBlock1 = (iBasis*spaceDim+1) * (numBasis*spaceDim);
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double Njp = basisDeriv[iQ+jBasis*spaceDim ];
+ const double Njq = basisDeriv[iQ+jBasis*spaceDim+1];
+ const double ki0j0 =
+ C1111 * Nip * Njp + C1112 * Niq * Njp +
+ C1112 * Nip * Njq + C1212 * Niq * Njq;
+ const double ki0j1 =
+ C1122 * Nip * Njq + C2212 * Niq * Njq +
+ C1112 * Nip * Njp + C1212 * Niq * Njp;
+ const double ki1j0 =
+ C1122 * Niq * Njp + C2212 * Niq * Njq +
+ C1112 * Nip * Njp + C1212 * Nip * Njq;
+ const double ki1j1 =
+ C2222 * Niq * Njq + C2212 * Nip * Njq +
+ C2212 * Niq * Njp + C1212 * Nip * Njp;
+ const int jBlock = (jBasis*spaceDim );
+ const int jBlock1 = (jBasis*spaceDim+1);
+ _cellMatrix[iBlock +jBlock ] += ki0j0;
+ _cellMatrix[iBlock +jBlock1] += ki0j1;
+ _cellMatrix[iBlock1+jBlock ] += ki1j0;
+ _cellMatrix[iBlock1+jBlock1] += ki1j1;
+ } // for
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*(3*11+4))));
+} // _elasticityJacobian2D
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in Jacobian for 3-D cells.
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_elasticityJacobian3D(
+ const double_array& elasticConsts)
+{ // _elasticityJacobian3D
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const double_array& quadWts = _quadrature->quadWts();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ assert(3 == cellDim);
+ assert(quadWts.size() == numQuadPts);
+ const int numConsts = 21;
+
+ // Compute Jacobian for consistent tangent matrix
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ // tau_ij = C_ijkl * e_kl
+ // = C_ijlk * 0.5 (u_k,l + u_l,k)
+ // = 0.5 * C_ijkl * (u_k,l + u_l,k)
+ // divide C_ijkl by 2 if k != l
+ const double C1111 = elasticConsts[iQuad*numConsts+ 0];
+ const double C1122 = elasticConsts[iQuad*numConsts+ 1];
+ const double C1133 = elasticConsts[iQuad*numConsts+ 2];
+ const double C1112 = elasticConsts[iQuad*numConsts+ 3] / 2.0;
+ const double C1123 = elasticConsts[iQuad*numConsts+ 4] / 2.0;
+ const double C1113 = elasticConsts[iQuad*numConsts+ 5] / 2.0;
+ const double C2222 = elasticConsts[iQuad*numConsts+ 6];
+ const double C2233 = elasticConsts[iQuad*numConsts+ 7];
+ const double C2212 = elasticConsts[iQuad*numConsts+ 8] / 2.0;
+ const double C2223 = elasticConsts[iQuad*numConsts+ 9] / 2.0;
+ const double C2213 = elasticConsts[iQuad*numConsts+10] / 2.0;
+ const double C3333 = elasticConsts[iQuad*numConsts+11];
+ const double C3312 = elasticConsts[iQuad*numConsts+12] / 2.0;
+ const double C3323 = elasticConsts[iQuad*numConsts+13] / 2.0;
+ const double C3313 = elasticConsts[iQuad*numConsts+14] / 2.0;
+ const double C1212 = elasticConsts[iQuad*numConsts+15] / 2.0;
+ const double C1223 = elasticConsts[iQuad*numConsts+16] / 2.0;
+ const double C1213 = elasticConsts[iQuad*numConsts+17] / 2.0;
+ const double C2323 = elasticConsts[iQuad*numConsts+18] / 2.0;
+ const double C2313 = elasticConsts[iQuad*numConsts+19] / 2.0;
+ const double C1313 = elasticConsts[iQuad*numConsts+20] / 2.0;
+ for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
+ iBasis < numBasis;
+ ++iBasis) {
+ const double Nip = wt*basisDeriv[iQ+iBasis*spaceDim+0];
+ const double Niq = wt*basisDeriv[iQ+iBasis*spaceDim+1];
+ const double Nir = wt*basisDeriv[iQ+iBasis*spaceDim+2];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const double Njp = basisDeriv[iQ+jBasis*spaceDim+0];
+ const double Njq = basisDeriv[iQ+jBasis*spaceDim+1];
+ const double Njr = basisDeriv[iQ+jBasis*spaceDim+2];
+ const double ki0j0 =
+ C1111 * Nip * Njp + C1112 * Niq * Njp + C1113 * Nir * Njp +
+ C1112 * Nip * Njq + C1212 * Niq * Njq + C1213 * Nir * Njq +
+ C1113 * Nip * Njr + C1213 * Niq * Njr + C1313 * Nir * Njr;
+ const double ki0j1 =
+ C1122 * Nip * Njq + C2212 * Niq * Njq + C2213 * Nir * Njq +
+ C1112 * Nip * Njp + C1212 * Niq * Njp + C1213 * Nir * Njp +
+ C1123 * Nip * Njr + C1223 * Niq * Njr + C2313 * Nir * Njr;
+ const double ki0j2 =
+ C1133 * Nip * Njr + C3312 * Niq * Njr + C3313 * Nir * Njr +
+ C1123 * Nip * Njq + C1223 * Niq * Njq + C2313 * Nir * Njq +
+ C1113 * Nip * Njp + C1213 * Niq * Njp + C1313 * Nir * Njp;
+ const double ki1j0 =
+ C1122 * Niq * Njp + C1112 * Nip * Njp + C1123 * Nir * Njp +
+ C2212 * Niq * Njq + C1212 * Nip * Njq + C1223 * Nir * Njq +
+ C2213 * Niq * Njr + C1213 * Nip * Njr + C2313 * Nir * Njr;
+ const double ki1j1 =
+ C2222 * Niq * Njq + C2212 * Nip * Njq + C2223 * Nir * Njq +
+ C2212 * Niq * Njp + C1212 * Nip * Njp + C1223 * Nir * Njp +
+ C2223 * Niq * Njr + C1223 * Nip * Njr + C2323 * Nir * Njr;
+ const double ki1j2 =
+ C2233 * Niq * Njr + C3312 * Nip * Njr + C3323 * Nir * Njr +
+ C2223 * Niq * Njq + C1223 * Nip * Njq + C2323 * Nir * Njq +
+ C2213 * Niq * Njp + C1213 * Nip * Njp + C2313 * Nir * Njp;
+ const double ki2j0 =
+ C1133 * Nir * Njp + C1123 * Niq * Njp + C1113 * Nip * Njp +
+ C3312 * Nir * Njq + C1223 * Niq * Njq + C1213 * Nip * Njq +
+ C3313 * Nir * Njr + C2313 * Niq * Njr + C1313 * Nip * Njr;
+ const double ki2j1 =
+ C2233 * Nir * Njq + C2223 * Niq * Njq + C2213 * Nip * Njq +
+ C3312 * Nir * Njp + C1223 * Niq * Njp + C1213 * Nip * Njp +
+ C3323 * Nir * Njr + C2323 * Niq * Njr + C2313 * Nip * Njr;
+ const double ki2j2 =
+ C3333 * Nir * Njr + C3323 * Niq * Njr + C3313 * Nip * Njr +
+ C3323 * Nir * Njq + C2323 * Niq * Njq + C2313 * Nip * Njq +
+ C3313 * Nir * Njp + C2313 * Niq * Njp + C1313 * Nip * Njp;
+ const int iBlock = iBasis*spaceDim * (numBasis*spaceDim);
+ const int iBlock1 = (iBasis*spaceDim+1) * (numBasis*spaceDim);
+ const int iBlock2 = (iBasis*spaceDim+2) * (numBasis*spaceDim);
+ const int jBlock = jBasis*spaceDim;
+ const int jBlock1 = jBasis*spaceDim+1;
+ const int jBlock2 = jBasis*spaceDim+2;
+ _cellMatrix[iBlock +jBlock ] += ki0j0;
+ _cellMatrix[iBlock +jBlock1] += ki0j1;
+ _cellMatrix[iBlock +jBlock2] += ki0j2;
+ _cellMatrix[iBlock1+jBlock ] += ki1j0;
+ _cellMatrix[iBlock1+jBlock1] += ki1j1;
+ _cellMatrix[iBlock1+jBlock2] += ki1j2;
+ _cellMatrix[iBlock2+jBlock ] += ki2j0;
+ _cellMatrix[iBlock2+jBlock1] += ki2j1;
+ _cellMatrix[iBlock2+jBlock2] += ki2j2;
+ } // for
+ } // for
+ } // for
+ PetscLogFlops(numQuadPts*(1+numBasis*(3+numBasis*(6*26+9))));
+} // _elasticityJacobian3D
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain1D(
+ double_array* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis,
+ const int numQuadPts)
+{ // calcTotalStrain1D
+ assert(0 != strain);
+
+ const int dim = 1;
+
+ assert(basisDeriv.size() == numQuadPts*numBasis*dim);
+ assert(disp.size() == numBasis*dim);
+
+ (*strain) = 0.0;
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ (*strain)[iQuad] += basisDeriv[iQuad*numBasis+iBasis] * disp[iBasis];
+ } // for
+} // calcTotalStrain1D
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain2D(
+ double_array* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis,
+ const int numQuadPts)
+{ // calcTotalStrain2D
+ assert(0 != strain);
+
+ const int dim = 2;
+
+ assert(basisDeriv.size() == numQuadPts*numBasis*dim);
+ assert(disp.size() == numBasis*dim);
+
+ (*strain) = 0.0;
+ const int strainSize = 3;
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
+ (*strain)[iQuad*strainSize+0] +=
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
+ (*strain)[iQuad*strainSize+1] +=
+ basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
+ (*strain)[iQuad*strainSize+2] +=
+ 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
+ } // for
+} // calcTotalStrain2D
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain3D(
+ double_array* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis,
+ const int numQuadPts)
+{ // calcTotalStrain3D
+ assert(0 != strain);
+
+ const int dim = 3;
+
+ assert(basisDeriv.size() == numQuadPts*numBasis*dim);
+ assert(disp.size() == numBasis*dim);
+
+ (*strain) = 0.0;
+ const int strainSize = 6;
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
+ (*strain)[iQuad*strainSize+0] +=
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
+ (*strain)[iQuad*strainSize+1] +=
+ basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
+ (*strain)[iQuad*strainSize+2] +=
+ basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+2];
+ (*strain)[iQuad*strainSize+3] +=
+ 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
+ (*strain)[iQuad*strainSize+4] +=
+ 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+1] +
+ basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+2]);
+ (*strain)[iQuad*strainSize+5] +=
+ 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim ] +
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+2]);
+ } // for
+} // calcTotalStrain3D
+
+// ----------------------------------------------------------------------
+// Calculate deformation tensor.
+void
+pylith::feassemble::IntegratorElasticityLgDeform::_calcDeformation(
+ double_array* deform,
+ const double_array& basisDeriv,
+ const double_array& coords,
+ const double_array& disp,
+ const int numBasis,
+ const int numQuadPts,
+ const int dim)
+{ // _calcDeformation
+ assert(0 != deform);
+
+ assert(basisDeriv.size() == numQuadPts*numBasis*dim);
+ assert(disp.size() == numBasis*dim);
+
+ const int deformSize = dim*dim;
+
+ (*deform) = 0.0;
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis)
+ for (int iDim=0, indexD=0; iDim < dim; ++iDim)
+ for (int jDim=0; jDim < dim; ++jDim, ++indexD)
+ (*deform)[iQuad*deformSize+indexD] +=
+ basisDeriv[iQ+iBasis*dim+jDim] *
+ (coords[iBasis*dim+iDim] + disp[iBasis*dim+iDim]);
+} // _calcDeformation
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.hh 2009-11-10 00:12:11 UTC (rev 15943)
@@ -0,0 +1,207 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/feassemble/IntegratorElasticityLgDeform.hh
+ *
+ * @brief Object containing general elasticity operations for implicit
+ * and explicit time integration of the elasticity equation for large
+ * rigid body rotations and small strains.
+ */
+
+#if !defined(pylith_feassemble_integratorelasticitylgdeform_hh)
+#define pylith_feassemble_integratorelasticitylgdeform_hh
+
+// Include directives ---------------------------------------------------
+#include "IntegratorElasticity.hh" // ISA IntegratorElasticity
+
+// IntegratorElasticity -------------------------------------------------
+/** @brief General elasticity operations for implicit and explicit
+ * time integration of the elasticity equation for large rigid body
+ * rotations and small strains.
+ */
+class pylith::feassemble::IntegratorElasticityLgDeform :
+ public IntegratorElasticity
+{ // IntegratorElasticityLgDeform
+ friend class TestIntegratorElasticityLgDeform; // unit testing
+
+// PUBLIC TYPEDEFS //////////////////////////////////////////////////////
+public :
+
+ typedef void (*totalStrain_fn_type)(double_array*,
+ const double_array&,
+ const double_array&,
+ const int,
+ const int);
+
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ IntegratorElasticityLgDeform(void);
+
+ /// Destructor
+ virtual
+ ~IntegratorElasticityLgDeform(void);
+
+ /** Determine whether we need to recompute the Jacobian.
+ *
+ * @returns True if Jacobian needs to be recomputed, false otherwise.
+ */
+ bool needNewJacobian(void);
+
+ /** Update state variables as needed.
+ *
+ * @param t Current time
+ * @param fields Solution fields
+ * @param mesh Finite-element mesh
+ */
+ void updateStateVars(const double t,
+ topology::SolutionFields* const fields);
+
+// PROTECTED METHODS ////////////////////////////////////////////////////
+protected :
+
+ /** Calculate stress or strain field from solution field.
+ *
+ * @param field Field in which to store stress or strain.
+ * @param name Name of field to compute ('total-strain' or 'stress')
+ * @param fields Manager for solution fields.
+ */
+ void _calcStrainStressField(topology::Field<topology::Mesh>* field,
+ const char* name,
+ topology::SolutionFields* const fields);
+
+ /** Calculate stress field from total strain field. Stress field
+ * replaces strain field in section.
+ *
+ * @param field Field in which to store stress.
+ */
+ void _calcStressFromStrain(topology::Field<topology::Mesh>* field);
+
+
+ /** Integrate elasticity term in residual for 1-D cells.
+ *
+ * @param stress Stress tensor for cell at quadrature points.
+ */
+ void _elasticityResidual1D(const double_array& stress);
+
+ /** Integrate elasticity term in residual for 2-D cells.
+ *
+ * @param stress Stress tensor for cell at quadrature points.
+ */
+ void _elasticityResidual2D(const double_array& stress);
+
+ /** Integrate elasticity term in residual for 3-D cells.
+ *
+ * @param stress Stress tensor for cell at quadrature points.
+ */
+ void _elasticityResidual3D(const double_array& stress);
+
+ /** Integrate elasticity term in Jacobian for 1-D cells.
+ *
+ * @param elasticConsts Matrix of elasticity constants at quadrature points.
+ */
+ void _elasticityJacobian1D(const double_array& elasticConsts);
+
+ /** Integrate elasticity term in Jacobian for 2-D cells.
+ *
+ * @param elasticConsts Matrix of elasticity constants at quadrature points.
+ */
+ void _elasticityJacobian2D(const double_array& elasticConsts);
+
+ /** Integrate elasticity term in Jacobian for 3-D cells.
+ *
+ * @param elasticConsts Matrix of elasticity constants at quadrature points.
+ */
+ void _elasticityJacobian3D(const double_array& elasticConsts);
+
+ /** Compute total strain in at quadrature points of a cell.
+ *
+ * @param strain Strain tensor at quadrature points.
+ * @param basisDeriv Derivatives of basis functions at quadrature points.
+ * @param disp Displacement at vertices of cell.
+ * @param numBasis Number of basis functions for cell.
+ * @param numQuadPts Number of quadrature points.
+ */
+ static
+ void _calcTotalStrain1D(double_array* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis,
+ const int numQuadPts);
+
+ /** Compute total strain in at quadrature points of a cell.
+ *
+ * @param strain Strain tensor at quadrature points.
+ * @param basisDeriv Derivatives of basis functions at quadrature points.
+ * @param disp Displacement at vertices of cell.
+ * @param numBasis Number of basis functions for cell.
+ * @param numQuadPts Number of quadrature points.
+ */
+ static
+ void _calcTotalStrain2D(double_array* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis,
+ const int numQuadPts);
+
+ /** Compute total strain in at quadrature points of a cell.
+ *
+ * @param strain Strain tensor at quadrature points.
+ * @param basisDeriv Derivatives of basis functions at quadrature points.
+ * @param disp Displacement at vertices of cell.
+ * @param numBasis Number of basis functions for cell.
+ * @param numQuadPts Number of quadrature points.
+ */
+ static
+ void _calcTotalStrain3D(double_array* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis,
+ const int numQuadPts);
+
+ /** Calculate deformation tensor.
+ *
+ * @param deform Deformation tensor for cell at quadrature points.
+ * @param basisDeriv Derivatives of basis functions at quadrature points.
+ * @param coords Coordinates of DOF of cell.
+ * @param disp Displacements of DOF of cell.
+ * @param numBasis Number of basis functions for cell.
+ * @param numQuadPts Number of quadrature points.
+ * @param dim Dimension of cell.
+ */
+ static
+ void _calcDeformation(double_array* deform,
+ const double_array& basisDeriv,
+ const double_array& coords,
+ const double_array& disp,
+ const int numBasis,
+ const int numQuadPts,
+ const int dim);
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ /// Not implemented.
+ IntegratorElasticityLgDeform(const IntegratorElasticityLgDeform&);
+
+ /// Not implemented
+ const IntegratorElasticityLgDeform& operator=(const IntegratorElasticityLgDeform&);
+
+}; // IntegratorElasticityLgDeform
+
+#endif // pylith_feassemble_integratorelasticitylgdeform_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-09 21:27:57 UTC (rev 15942)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Makefile.am 2009-11-10 00:12:11 UTC (rev 15943)
@@ -24,6 +24,7 @@
Integrator.icc \
Integrator.cc \
IntegratorElasticity.hh \
+ IntegratorElasticityLgDeform.hh \
GeometryPoint1D.hh \
GeometryPoint2D.hh \
GeometryPoint3D.hh \
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/feassemblefwd.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/feassemblefwd.hh 2009-11-09 21:27:57 UTC (rev 15942)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/feassemblefwd.hh 2009-11-10 00:12:11 UTC (rev 15943)
@@ -56,6 +56,10 @@
class ElasticityImplicit;
class ElasticityExplicit;
+ class IntegratorElasticityLgDeform;
+ class ElasticityImplicitLgDeform;
+ class ElasticityExplicitLgDeform;
+
} // feassemble
} // pylith
More information about the CIG-COMMITS
mailing list