[cig-commits] r8299 - in short/3D/PyLith/trunk: libsrc
libsrc/feassemble unittests/libtests/feassemble
unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Fri Nov 16 15:40:39 PST 2007
Author: brad
Date: 2007-11-16 15:40:38 -0800 (Fri, 16 Nov 2007)
New Revision: 8299
Added:
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py
Log:
Factored out common routines from ElasticityImplicit and ElasticityExplicit into IntegratorElasticity.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2007-11-16 23:40:38 UTC (rev 8299)
@@ -52,6 +52,7 @@
feassemble/GeometryQuad3D.cc \
feassemble/GeometryHex3D.cc \
feassemble/Integrator.cc \
+ feassemble/IntegratorElasticity.cc \
feassemble/Quadrature.cc \
feassemble/Quadrature0D.cc \
feassemble/Quadrature1D.cc \
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2007-11-16 23:40:38 UTC (rev 8299)
@@ -34,8 +34,7 @@
// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityExplicit::ElasticityExplicit(void) :
- _dtm1(-1.0),
- _material(0)
+ _dtm1(-1.0)
{ // constructor
} // constructor
@@ -43,7 +42,6 @@
// Destructor
pylith::feassemble::ElasticityExplicit::~ElasticityExplicit(void)
{ // destructor
- _material = 0; // Don't manage memory for material
} // destructor
// ----------------------------------------------------------------------
@@ -62,35 +60,14 @@
} // timeStep
// ----------------------------------------------------------------------
-// Set material.
-void
-pylith::feassemble::ElasticityExplicit::material(materials::ElasticMaterial* m)
-{ // material
- _material = m;
- if (0 != _material)
- _material->timeStep(_dt);
-} // material
-
-// ----------------------------------------------------------------------
-// Determine whether we need to recompute the Jacobian.
-bool
-pylith::feassemble::ElasticityExplicit::needNewJacobian(void)
-{ // needNewJacobian
- assert(0 != _material);
- if (!_needNewJacobian)
- _needNewJacobian = _material->needNewJacobian();
- return _needNewJacobian;
-} // needNewJacobian
-
-// ----------------------------------------------------------------------
// Set flag for setting constraints for total field solution or
// incremental field solution.
void
pylith::feassemble::ElasticityExplicit::useSolnIncr(const bool flag)
{ // useSolnIncr
- assert(0 != _material);
- _useSolnIncr = flag;
- _material->useElasticBehavior(!_useSolnIncr);
+ throw std::logic_error("Incremental solution not supported for "
+ "explicit time integration of elasticity "
+ "equation.");
} // useSolnIncr
// ----------------------------------------------------------------------
@@ -371,459 +348,5 @@
_material->resetNeedNewJacobian();
} // integrateJacobian
-// ----------------------------------------------------------------------
-// Update state variables as needed.
-void
-pylith::feassemble::ElasticityExplicit::updateState(
- const double t,
- const ALE::Obj<real_section_type>& disp,
- const ALE::Obj<Mesh>& mesh)
-{ // updateState
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(!disp.isNull());
- // No need to update state if using elastic behavior
- if (!_material->usesUpdateState())
- return;
-
- // Set variables dependent on dimension of cell
- const int cellDim = _quadrature->cellDim();
- int tensorSize = 0;
- Elasticity::totalStrain_fn_type calcTotalStrainFn;
- if (1 == cellDim) {
- tensorSize = 1;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain1D;
- } else if (2 == cellDim) {
- tensorSize = 3;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain2D;
- } else if (3 == cellDim) {
- tensorSize = 6;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain3D;
- } else
- assert(0);
-
- // Get cell information
- const ALE::Obj<ALE::Mesh::label_sequence>& cells =
- mesh->getLabelStratum("material-id", _material->id());
- assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
-
- // Get sections
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
-
- // Get cell geometry information that doesn't depend on cell
- const int numQuadPts = _quadrature->numQuadPts();
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
-
- const int cellVecSize = numBasis*spaceDim;
- double_array dispCell(cellVecSize);
-
- // Allocate vector for total strain
- std::vector<double_array> totalStrain(numQuadPts);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- totalStrain[iQuad].resize(tensorSize);
- totalStrain[iQuad] = 0.0;
- } // for
-
-#ifdef FASTER
- const int dispTTag = _dispTTags[_material->id()];
-#endif
-
- int c_index = 0;
- for (Mesh::label_sequence::iterator c_iter=cells->begin();
- c_iter != cellsEnd;
- ++c_iter, ++c_index) {
- // Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
-
- // Restrict input fields to cell
-#ifdef FASTER
- mesh->restrict(disp, dispTTag, c_index, &dispCell[0], cellVecSize);
-#else
- mesh->restrict(disp, *c_iter, &dispCell[0], cellVecSize);
-#endif
-
- // Get cell geometry information that depends on cell
- const double_array& basisDeriv = _quadrature->basisDeriv();
-
- // Compute strains
- calcTotalStrainFn(&totalStrain, basisDeriv, dispCell, numBasis);
-
- // Update material state
- _material->updateState(totalStrain, *c_iter);
- } // for
-} // updateState
-
-// ----------------------------------------------------------------------
-// Verify configuration is acceptable.
-void
-pylith::feassemble::ElasticityExplicit::verifyConfiguration(
- const ALE::Obj<Mesh>& mesh)
-{ // verifyConfiguration
- assert(0 != _quadrature);
- assert(0 != _material);
-
- const int dimension = mesh->getDimension();
-
- // check compatibility of mesh and material
- if (_material->dimension() != dimension) {
- std::ostringstream msg;
- msg << "Material '" << _material->label()
- << "' is incompatible with mesh.\n"
- << "Dimension of mesh: " << dimension
- << ", dimension of material: " << _material->dimension()
- << ".";
- throw std::runtime_error(msg.str());
- } // if
-
- // check compatibility of mesh and quadrature scheme
- if (_quadrature->cellDim() != dimension) {
- std::ostringstream msg;
- msg << "Quadrature is incompatible with cells for material '"
- << _material->label() << "'.\n"
- << "Dimension of mesh: " << dimension
- << ", dimension of quadrature: " << _quadrature->cellDim()
- << ".";
- throw std::runtime_error(msg.str());
- } // if
- const int numCorners = _quadrature->refGeometry().numCorners();
- const ALE::Obj<ALE::Mesh::label_sequence>& cells =
- mesh->getLabelStratum("material-id", _material->id());
- assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
- const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
- assert(!sieve.isNull());
- for (Mesh::label_sequence::iterator c_iter=cells->begin();
- c_iter != cellsEnd;
- ++c_iter) {
- const int cellNumCorners = sieve->nCone(*c_iter, mesh->depth())->size();
- if (numCorners != cellNumCorners) {
- std::ostringstream msg;
- msg << "Quadrature is incompatible with cell in material '"
- << _material->label() << "'.\n"
- << "Cell " << *c_iter << " has " << cellNumCorners
- << " vertices but quadrature reference cell has "
- << numCorners << " vertices.";
- throw std::runtime_error(msg.str());
- } // if
- } // for
-} // verifyConfiguration
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in residual for 1-D cells.
-void
-pylith::feassemble::ElasticityExplicit::_elasticityResidual1D(
- const std::vector<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][0];
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- const double N1 = wt*basisDeriv[iQuad*numBasis+iBasis ];
- _cellVector[iBasis*spaceDim ] -= N1*s11;
- } // for
- } // for
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*5));
-} // _elasticityResidual1D
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in residual for 2-D cells.
-void
-pylith::feassemble::ElasticityExplicit::_elasticityResidual2D(
- const std::vector<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);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- const double s22 = stress[iQuad][1];
- const double s12 = stress[iQuad][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
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(8+2+9)));
-} // _elasticityResidual2D
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in residual for 3-D cells.
-void
-pylith::feassemble::ElasticityExplicit::_elasticityResidual3D(
- const std::vector<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);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- const double s22 = stress[iQuad][1];
- const double s33 = stress[iQuad][2];
- const double s12 = stress[iQuad][3];
- const double s23 = stress[iQuad][4];
- const double s13 = stress[iQuad][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
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(3+12)));
-} // _elasticityResidual3D
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in Jacobian for 1-D cells.
-void
-pylith::feassemble::ElasticityExplicit::_elasticityJacobian1D(
- const std::vector<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][0];
- 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
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(2+numBasis*3)));
-} // _elasticityJacobian1D
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in Jacobian for 2-D cells.
-void
-pylith::feassemble::ElasticityExplicit::_elasticityJacobian2D(
- const std::vector<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);
-
- 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][0];
- const double C1122 = elasticConsts[iQuad][1];
- const double C1112 = elasticConsts[iQuad][2]/2.0;
- const double C2222 = elasticConsts[iQuad][3];
- const double C2212 = elasticConsts[iQuad][4]/2.0;
- const double C1212 = elasticConsts[iQuad][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
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(2+numBasis*(3*11+4))));
-} // _elasticityJacobian2D
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in Jacobian for 3-D cells.
-void
-pylith::feassemble::ElasticityExplicit::_elasticityJacobian3D(
- const std::vector<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);
-
- // 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][ 0];
- const double C1122 = elasticConsts[iQuad][ 1];
- const double C1133 = elasticConsts[iQuad][ 2];
- const double C1112 = elasticConsts[iQuad][ 3]/2.0;
- const double C1123 = elasticConsts[iQuad][ 4]/2.0;
- const double C1113 = elasticConsts[iQuad][ 5]/2.0;
- const double C2222 = elasticConsts[iQuad][ 6];
- const double C2233 = elasticConsts[iQuad][ 7];
- const double C2212 = elasticConsts[iQuad][ 8]/2.0;
- const double C2223 = elasticConsts[iQuad][ 9]/2.0;
- const double C2213 = elasticConsts[iQuad][10]/2.0;
- const double C3333 = elasticConsts[iQuad][11];
- const double C3312 = elasticConsts[iQuad][12]/2.0;
- const double C3323 = elasticConsts[iQuad][13]/2.0;
- const double C3313 = elasticConsts[iQuad][14]/2.0;
- const double C1212 = elasticConsts[iQuad][15]/2.0;
- const double C1223 = elasticConsts[iQuad][16]/2.0;
- const double C1213 = elasticConsts[iQuad][17]/2.0;
- const double C2323 = elasticConsts[iQuad][18]/2.0;
- const double C2313 = elasticConsts[iQuad][19]/2.0;
- const double C1313 = elasticConsts[iQuad][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
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(3+numBasis*(6*26+9))));
-} // _elasticityJacobian3D
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2007-11-16 23:40:38 UTC (rev 8299)
@@ -55,7 +55,7 @@
#if !defined(pylith_feassemble_elasticityexplicit_hh)
#define pylith_feassemble_elasticityexplicit_hh
-#include "Integrator.hh" // ISA Integrator
+#include "IntegratorElasticity.hh" // ISA IntegratorElasticity
#include "pylith/utils/array.hh" // USES std::vector, double_array
namespace pylith {
@@ -63,10 +63,6 @@
class ElasticityExplicit;
class TestElasticityExplicit;
} // feassemble
-
- namespace materials {
- class ElasticMaterial;
- } // feassemble
} // pylith
namespace spatialdata {
@@ -78,7 +74,7 @@
} // geocoords
} // spatialdata
-class pylith::feassemble::ElasticityExplicit : public Integrator
+class pylith::feassemble::ElasticityExplicit : public IntegratorElasticity
{ // ElasticityExplicit
friend class TestElasticityExplicit; // unit testing
@@ -91,24 +87,12 @@
/// Destructor
~ElasticityExplicit(void);
- /** Set material.
- *
- * @param m Elastic material.
- */
- void material(materials::ElasticMaterial* m);
-
/** Set time step for advancing from time t to time t+dt.
*
* @param dt Time step
*/
void timeStep(const double dt);
- /** Determine whether we need to recompute the Jacobian.
- *
- * @returns True if Jacobian needs to be recomputed, false otherwise.
- */
- bool needNewJacobian(void);
-
/** Set flag for setting constraints for total field solution or
* incremental field solution.
*
@@ -141,64 +125,9 @@
topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh);
- /** Update state variables as needed.
- *
- * @param t Current time
- * @param field Current solution field.
- * @param mesh Finite-element mesh
- */
- void updateState(const double t,
- const ALE::Obj<real_section_type>& field,
- const ALE::Obj<Mesh>& mesh);
-
- /** Verify configuration is acceptable.
- *
- * @param mesh Finite-element mesh
- */
- void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
-
// PRIVATE METHODS //////////////////////////////////////////////////////
private :
- /** Integrate elasticity term in residual for 1-D cells.
- *
- * @param stress Stress tensor for cell at quadrature points.
- */
- void _elasticityResidual1D(const std::vector<double_array>& stress);
-
- /** Integrate elasticity term in residual for 2-D cells.
- *
- * @param stress Stress tensor for cell at quadrature points.
- */
- void _elasticityResidual2D(const std::vector<double_array>& stress);
-
- /** Integrate elasticity term in residual for 3-D cells.
- *
- * @param stress Stress tensor for cell at quadrature points.
- */
- void _elasticityResidual3D(const std::vector<double_array>& stress);
-
- /** Integrate elasticity term in Jacobian for 1-D cells.
- *
- * @param elasticConsts Matrix of elasticity constants at quadrature points.
- */
- void _elasticityJacobian1D(const std::vector<double_array>& elasticConsts);
-
- /** Integrate elasticity term in Jacobian for 2-D cells.
- *
- * @param elasticConsts Matrix of elasticity constants at quadrature points.
- */
- void _elasticityJacobian2D(const std::vector<double_array>& elasticConsts);
-
- /** Integrate elasticity term in Jacobian for 3-D cells.
- *
- * @param elasticConsts Matrix of elasticity constants at quadrature points.
- */
- void _elasticityJacobian3D(const std::vector<double_array>& elasticConsts);
-
-// PRIVATE METHODS //////////////////////////////////////////////////////
-private :
-
/// Not implemented.
ElasticityExplicit(const ElasticityExplicit& i);
@@ -210,9 +139,6 @@
double _dtm1; ///< Time step for t-dt1 -> t
- /// Elastic material associated with integrator
- materials::ElasticMaterial* _material;
-
// Optimization
std::map<int, int> _dispTTags; ///< Tags indexing dispT field
std::map<int, int> _dispTmdtTags; ///< Tags indexing dispTmdt field
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-11-16 23:40:38 UTC (rev 8299)
@@ -34,8 +34,7 @@
// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityImplicit::ElasticityImplicit(void) :
- _dtm1(-1.0),
- _material(0)
+ _dtm1(-1.0)
{ // constructor
} // constructor
@@ -43,7 +42,6 @@
// Destructor
pylith::feassemble::ElasticityImplicit::~ElasticityImplicit(void)
{ // destructor
- _material = 0; // Don't manage memory for material
} // destructor
// ----------------------------------------------------------------------
@@ -62,27 +60,6 @@
} // timeStep
// ----------------------------------------------------------------------
-// Set material.
-void
-pylith::feassemble::ElasticityImplicit::material(materials::ElasticMaterial* m)
-{ // material
- _material = m;
- if (0 != _material)
- _material->timeStep(_dt);
-} // material
-
-// ----------------------------------------------------------------------
-// Determine whether we need to recompute the Jacobian.
-bool
-pylith::feassemble::ElasticityImplicit::needNewJacobian(void)
-{ // needNewJacobian
- assert(0 != _material);
- if (!_needNewJacobian)
- _needNewJacobian = _material->needNewJacobian();
- return _needNewJacobian;
-} // needNewJacobian
-
-// ----------------------------------------------------------------------
// Set flag for setting constraints for total field solution or
// incremental field solution.
void
@@ -435,464 +412,5 @@
_material->resetNeedNewJacobian();
} // integrateJacobian
-// ----------------------------------------------------------------------
-// Update state variables as needed.
-void
-pylith::feassemble::ElasticityImplicit::updateState(
- const double t,
- const ALE::Obj<real_section_type>& disp,
- const ALE::Obj<Mesh>& mesh)
-{ // updateState
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(!disp.isNull());
- // No need to update state if using elastic behavior
- if (!_material->usesUpdateState())
- return;
-
- // Set variables dependent on dimension of cell
- const int cellDim = _quadrature->cellDim();
- int tensorSize = 0;
- Elasticity::totalStrain_fn_type calcTotalStrainFn;
- if (1 == cellDim) {
- tensorSize = 1;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain1D;
- } else if (2 == cellDim) {
- tensorSize = 3;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain2D;
- } else if (3 == cellDim) {
- tensorSize = 6;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain3D;
- } else
- assert(0);
-
- // Get cell information
- const ALE::Obj<ALE::Mesh::label_sequence>& cells =
- mesh->getLabelStratum("material-id", _material->id());
- assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
-
- // Get sections
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
-
- // Get cell geometry information that doesn't depend on cell
- const int numQuadPts = _quadrature->numQuadPts();
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
-
- // Precompute the geometric and function space information
- _quadrature->precomputeGeometry(mesh, coordinates, cells);
-
- const int cellVecSize = numBasis*spaceDim;
- double_array dispCell(cellVecSize);
-
- // Allocate vector for total strain
- std::vector<double_array> totalStrain(numQuadPts);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- totalStrain[iQuad].resize(tensorSize);
- totalStrain[iQuad] = 0.0;
- } // for
-
-#ifdef FASTER
- const int dTag = _dTags[_material->id()];
-#endif
-
- int c_index = 0;
- for (Mesh::label_sequence::iterator c_iter=cells->begin();
- c_iter != cellsEnd;
- ++c_iter, ++c_index) {
- // Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
-
- // Restrict input fields to cell
-#ifdef FASTER
- mesh->restrict(disp, dTag, c_index, &dispCell[0], cellVecSize);
-#else
- mesh->restrict(disp, *c_iter, &dispCell[0], cellVecSize);
-#endif
-
- // Get cell geometry information that depends on cell
- const double_array& basisDeriv = _quadrature->basisDeriv();
-
- // Compute strains
- calcTotalStrainFn(&totalStrain, basisDeriv, dispCell, numBasis);
-
- // Update material state
- _material->updateState(totalStrain, *c_iter);
- } // for
-
- _material->useElasticBehavior(false);
-} // updateState
-
-// ----------------------------------------------------------------------
-// Verify configuration is acceptable.
-void
-pylith::feassemble::ElasticityImplicit::verifyConfiguration(
- const ALE::Obj<Mesh>& mesh)
-{ // verifyConfiguration
- assert(0 != _quadrature);
- assert(0 != _material);
-
- const int dimension = mesh->getDimension();
-
- // check compatibility of mesh and material
- if (_material->dimension() != dimension) {
- std::ostringstream msg;
- msg << "Material '" << _material->label()
- << "' is incompatible with mesh.\n"
- << "Dimension of mesh: " << dimension
- << ", dimension of material: " << _material->dimension()
- << ".";
- throw std::runtime_error(msg.str());
- } // if
-
- // check compatibility of mesh and quadrature scheme
- if (_quadrature->cellDim() != dimension) {
- std::ostringstream msg;
- msg << "Quadrature is incompatible with cells for material '"
- << _material->label() << "'.\n"
- << "Dimension of mesh: " << dimension
- << ", dimension of quadrature: " << _quadrature->cellDim()
- << ".";
- throw std::runtime_error(msg.str());
- } // if
- const int numCorners = _quadrature->refGeometry().numCorners();
- const ALE::Obj<ALE::Mesh::label_sequence>& cells =
- mesh->getLabelStratum("material-id", _material->id());
- assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
- const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
- assert(!sieve.isNull());
- for (Mesh::label_sequence::iterator c_iter=cells->begin();
- c_iter != cellsEnd;
- ++c_iter) {
- const int cellNumCorners = sieve->nCone(*c_iter, mesh->depth())->size();
- if (numCorners != cellNumCorners) {
- std::ostringstream msg;
- msg << "Quadrature is incompatible with cell in material '"
- << _material->label() << "'.\n"
- << "Cell " << *c_iter << " has " << cellNumCorners
- << " vertices but quadrature reference cell has "
- << numCorners << " vertices.";
- throw std::runtime_error(msg.str());
- } // if
- } // for
-} // verifyConfiguration
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in residual for 1-D cells.
-void
-pylith::feassemble::ElasticityImplicit::_elasticityResidual1D(
- const std::vector<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][0];
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- const double N1 = wt*basisDeriv[iQuad*numBasis+iBasis ];
- _cellVector[iBasis*spaceDim ] -= N1*s11;
- } // for
- } // for
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*5));
-} // _elasticityResidual1D
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in residual for 2-D cells.
-void
-pylith::feassemble::ElasticityImplicit::_elasticityResidual2D(
- const std::vector<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);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- const double s22 = stress[iQuad][1];
- const double s12 = stress[iQuad][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
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(8+2+9)));
-} // _elasticityResidual2D
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in residual for 3-D cells.
-void
-pylith::feassemble::ElasticityImplicit::_elasticityResidual3D(
- const std::vector<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);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- const double s22 = stress[iQuad][1];
- const double s33 = stress[iQuad][2];
- const double s12 = stress[iQuad][3];
- const double s23 = stress[iQuad][4];
- const double s13 = stress[iQuad][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
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(3+12)));
-} // _elasticityResidual3D
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in Jacobian for 1-D cells.
-void
-pylith::feassemble::ElasticityImplicit::_elasticityJacobian1D(
- const std::vector<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][0];
- 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
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(2+numBasis*3)));
-} // _elasticityJacobian1D
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in Jacobian for 2-D cells.
-void
-pylith::feassemble::ElasticityImplicit::_elasticityJacobian2D(
- const std::vector<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);
-
- 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][0];
- const double C1122 = elasticConsts[iQuad][1];
- const double C1112 = elasticConsts[iQuad][2]/2.0;
- const double C2222 = elasticConsts[iQuad][3];
- const double C2212 = elasticConsts[iQuad][4]/2.0;
- const double C1212 = elasticConsts[iQuad][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
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(2+numBasis*(3*11+4))));
-} // _elasticityJacobian2D
-
-// ----------------------------------------------------------------------
-// Integrate elasticity term in Jacobian for 3-D cells.
-void
-pylith::feassemble::ElasticityImplicit::_elasticityJacobian3D(
- const std::vector<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);
-
- // 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][ 0];
- const double C1122 = elasticConsts[iQuad][ 1];
- const double C1133 = elasticConsts[iQuad][ 2];
- const double C1112 = elasticConsts[iQuad][ 3]/2.0;
- const double C1123 = elasticConsts[iQuad][ 4]/2.0;
- const double C1113 = elasticConsts[iQuad][ 5]/2.0;
- const double C2222 = elasticConsts[iQuad][ 6];
- const double C2233 = elasticConsts[iQuad][ 7];
- const double C2212 = elasticConsts[iQuad][ 8]/2.0;
- const double C2223 = elasticConsts[iQuad][ 9]/2.0;
- const double C2213 = elasticConsts[iQuad][10]/2.0;
- const double C3333 = elasticConsts[iQuad][11];
- const double C3312 = elasticConsts[iQuad][12]/2.0;
- const double C3323 = elasticConsts[iQuad][13]/2.0;
- const double C3313 = elasticConsts[iQuad][14]/2.0;
- const double C1212 = elasticConsts[iQuad][15]/2.0;
- const double C1223 = elasticConsts[iQuad][16]/2.0;
- const double C1213 = elasticConsts[iQuad][17]/2.0;
- const double C2323 = elasticConsts[iQuad][18]/2.0;
- const double C2313 = elasticConsts[iQuad][19]/2.0;
- const double C1313 = elasticConsts[iQuad][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
- PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(3+numBasis*(6*26+9))));
-} // _elasticityJacobian3D
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh 2007-11-16 23:40:38 UTC (rev 8299)
@@ -47,7 +47,7 @@
#if !defined(pylith_feassemble_elasticityimplicit_hh)
#define pylith_feassemble_elasticityimplicit_hh
-#include "Integrator.hh" // ISA Integrator
+#include "IntegratorElasticity.hh" // ISA IntegratorElasticity
#include "pylith/utils/array.hh" // USES std::vector, double_array
namespace pylith {
@@ -55,10 +55,6 @@
class ElasticityImplicit;
class TestElasticityImplicit;
} // feassemble
-
- namespace materials {
- class ElasticMaterial;
- } // feassemble
} // pylith
namespace spatialdata {
@@ -70,7 +66,7 @@
} // geocoords
} // spatialdata
-class pylith::feassemble::ElasticityImplicit : public Integrator
+class pylith::feassemble::ElasticityImplicit : public IntegratorElasticity
{ // ElasticityImplicit
friend class TestElasticityImplicit; // unit testing
@@ -83,24 +79,12 @@
/// Destructor
~ElasticityImplicit(void);
- /** Set material.
- *
- * @param m Elastic material.
- */
- void material(materials::ElasticMaterial* m);
-
/** Set time step for advancing from time t to time t+dt.
*
* @param dt Time step
*/
void timeStep(const double dt);
- /** Determine whether we need to recompute the Jacobian.
- *
- * @returns True if Jacobian needs to be recomputed, false otherwise.
- */
- bool needNewJacobian(void);
-
/** Set flag for setting constraints for total field solution or
* incremental field solution.
*
@@ -140,61 +124,6 @@
topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh);
- /** Update state variables as needed.
- *
- * @param t Current time
- * @param field Current solution field.
- * @param mesh Finite-element mesh
- */
- void updateState(const double t,
- const ALE::Obj<real_section_type>& field,
- const ALE::Obj<Mesh>& mesh);
-
- /** Verify configuration is acceptable.
- *
- * @param mesh Finite-element mesh
- */
- void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
-
-// PRIVATE METHODS //////////////////////////////////////////////////////
-private :
-
- /** Integrate elasticity term in residual for 1-D cells.
- *
- * @param stress Stress tensor for cell at quadrature points.
- */
- void _elasticityResidual1D(const std::vector<double_array>& stress);
-
- /** Integrate elasticity term in residual for 2-D cells.
- *
- * @param stress Stress tensor for cell at quadrature points.
- */
- void _elasticityResidual2D(const std::vector<double_array>& stress);
-
- /** Integrate elasticity term in residual for 3-D cells.
- *
- * @param stress Stress tensor for cell at quadrature points.
- */
- void _elasticityResidual3D(const std::vector<double_array>& stress);
-
- /** Integrate elasticity term in Jacobian for 1-D cells.
- *
- * @param elasticConsts Matrix of elasticity constants at quadrature points.
- */
- void _elasticityJacobian1D(const std::vector<double_array>& elasticConsts);
-
- /** Integrate elasticity term in Jacobian for 2-D cells.
- *
- * @param elasticConsts Matrix of elasticity constants at quadrature points.
- */
- void _elasticityJacobian2D(const std::vector<double_array>& elasticConsts);
-
- /** Integrate elasticity term in Jacobian for 3-D cells.
- *
- * @param elasticConsts Matrix of elasticity constants at quadrature points.
- */
- void _elasticityJacobian3D(const std::vector<double_array>& elasticConsts);
-
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
@@ -209,9 +138,6 @@
double _dtm1; ///< Time step for t-dt1 -> t
- /// Elastic material associated with integrator
- materials::ElasticMaterial* _material;
-
// Optimization
std::map<int, int> _dTags; ///< Tags indexing dispTBctpdt field
std::map<int, int> _rTags; ///< tags indexing residual field
Added: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2007-11-16 23:40:38 UTC (rev 8299)
@@ -0,0 +1,521 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "IntegratorElasticity.hh" // implementation of class methods
+
+#include "Quadrature.hh" // USES Quadrature
+#include "Elasticity.hh" // USES Elasticity
+#include "CellGeometry.hh" // USES CellGeometry
+
+#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
+#include "pylith/utils/array.hh" // USES double_array
+
+#include <assert.h> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+
+//#define FASTER
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::IntegratorElasticity::IntegratorElasticity(void) :
+ _material(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::IntegratorElasticity::~IntegratorElasticity(void)
+{ // destructor
+ _material = 0; // Don't manage memory for material
+} // destructor
+
+// ----------------------------------------------------------------------
+// Set material.
+void
+pylith::feassemble::IntegratorElasticity::material(materials::ElasticMaterial* m)
+{ // material
+ _material = m;
+ if (0 != _material)
+ _material->timeStep(_dt);
+} // material
+
+// ----------------------------------------------------------------------
+// Determine whether we need to recompute the Jacobian.
+bool
+pylith::feassemble::IntegratorElasticity::needNewJacobian(void)
+{ // needNewJacobian
+ assert(0 != _material);
+ if (!_needNewJacobian)
+ _needNewJacobian = _material->needNewJacobian();
+ return _needNewJacobian;
+} // needNewJacobian
+
+// ----------------------------------------------------------------------
+// Update state variables as needed.
+void
+pylith::feassemble::IntegratorElasticity::updateState(
+ const double t,
+ const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh)
+{ // updateState
+ assert(0 != _quadrature);
+ assert(0 != _material);
+ assert(!disp.isNull());
+
+ // No need to update state if using elastic behavior
+ if (!_material->usesUpdateState())
+ return;
+
+ // Set variables dependent on dimension of cell
+ const int cellDim = _quadrature->cellDim();
+ int tensorSize = 0;
+ Elasticity::totalStrain_fn_type calcTotalStrainFn;
+ if (1 == cellDim) {
+ tensorSize = 1;
+ calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain1D;
+ } else if (2 == cellDim) {
+ tensorSize = 3;
+ calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain2D;
+ } else if (3 == cellDim) {
+ tensorSize = 6;
+ calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain3D;
+ } else
+ assert(0);
+
+ // Get cell information
+ const ALE::Obj<ALE::Mesh::label_sequence>& cells =
+ mesh->getLabelStratum("material-id", _material->id());
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get sections
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+
+ const int cellVecSize = numBasis*spaceDim;
+ double_array dispCell(cellVecSize);
+
+ // Allocate vector for total strain
+ std::vector<double_array> totalStrain(numQuadPts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ totalStrain[iQuad].resize(tensorSize);
+ totalStrain[iQuad] = 0.0;
+ } // for
+
+#ifdef FASTER
+ const int dispTTag = _dispTTags[_material->id()];
+#endif
+
+ int c_index = 0;
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter, ++c_index) {
+ // Compute geometry information for current cell
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
+
+ // Restrict input fields to cell
+#ifdef FASTER
+ mesh->restrict(disp, dispTTag, c_index, &dispCell[0], cellVecSize);
+#else
+ mesh->restrict(disp, *c_iter, &dispCell[0], cellVecSize);
+#endif
+
+ // Get cell geometry information that depends on cell
+ const double_array& basisDeriv = _quadrature->basisDeriv();
+
+ // Compute strains
+ calcTotalStrainFn(&totalStrain, basisDeriv, dispCell, numBasis);
+
+ // Update material state
+ _material->updateState(totalStrain, *c_iter);
+ } // for
+
+ _material->useElasticBehavior(false);
+} // updateState
+
+// ----------------------------------------------------------------------
+// Verify configuration is acceptable.
+void
+pylith::feassemble::IntegratorElasticity::verifyConfiguration(
+ const ALE::Obj<Mesh>& mesh)
+{ // verifyConfiguration
+ assert(0 != _quadrature);
+ assert(0 != _material);
+
+ const int dimension = mesh->getDimension();
+
+ // check compatibility of mesh and material
+ if (_material->dimension() != dimension) {
+ std::ostringstream msg;
+ msg << "Material '" << _material->label()
+ << "' is incompatible with mesh.\n"
+ << "Dimension of mesh: " << dimension
+ << ", dimension of material: " << _material->dimension()
+ << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ // check compatibility of mesh and quadrature scheme
+ if (_quadrature->cellDim() != dimension) {
+ std::ostringstream msg;
+ msg << "Quadrature is incompatible with cells for material '"
+ << _material->label() << "'.\n"
+ << "Dimension of mesh: " << dimension
+ << ", dimension of quadrature: " << _quadrature->cellDim()
+ << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ const int numCorners = _quadrature->refGeometry().numCorners();
+ const ALE::Obj<ALE::Mesh::label_sequence>& cells =
+ mesh->getLabelStratum("material-id", _material->id());
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+ const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
+ assert(!sieve.isNull());
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ const int cellNumCorners = sieve->nCone(*c_iter, mesh->depth())->size();
+ if (numCorners != cellNumCorners) {
+ std::ostringstream msg;
+ msg << "Quadrature is incompatible with cell in material '"
+ << _material->label() << "'.\n"
+ << "Cell " << *c_iter << " has " << cellNumCorners
+ << " vertices but quadrature reference cell has "
+ << numCorners << " vertices.";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // for
+} // verifyConfiguration
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in residual for 1-D cells.
+void
+pylith::feassemble::IntegratorElasticity::_elasticityResidual1D(
+ const std::vector<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][0];
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const double N1 = wt*basisDeriv[iQuad*numBasis+iBasis ];
+ _cellVector[iBasis*spaceDim ] -= N1*s11;
+ } // for
+ } // for
+ PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*5));
+} // _elasticityResidual1D
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in residual for 2-D cells.
+void
+pylith::feassemble::IntegratorElasticity::_elasticityResidual2D(
+ const std::vector<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);
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double s11 = stress[iQuad][0];
+ const double s22 = stress[iQuad][1];
+ const double s12 = stress[iQuad][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
+ PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(8+2+9)));
+} // _elasticityResidual2D
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in residual for 3-D cells.
+void
+pylith::feassemble::IntegratorElasticity::_elasticityResidual3D(
+ const std::vector<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);
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double s11 = stress[iQuad][0];
+ const double s22 = stress[iQuad][1];
+ const double s33 = stress[iQuad][2];
+ const double s12 = stress[iQuad][3];
+ const double s23 = stress[iQuad][4];
+ const double s13 = stress[iQuad][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
+ PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(3+12)));
+} // _elasticityResidual3D
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in Jacobian for 1-D cells.
+void
+pylith::feassemble::IntegratorElasticity::_elasticityJacobian1D(
+ const std::vector<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][0];
+ 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
+ PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(2+numBasis*3)));
+} // _elasticityJacobian1D
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in Jacobian for 2-D cells.
+void
+pylith::feassemble::IntegratorElasticity::_elasticityJacobian2D(
+ const std::vector<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);
+
+ 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][0];
+ const double C1122 = elasticConsts[iQuad][1];
+ const double C1112 = elasticConsts[iQuad][2]/2.0;
+ const double C2222 = elasticConsts[iQuad][3];
+ const double C2212 = elasticConsts[iQuad][4]/2.0;
+ const double C1212 = elasticConsts[iQuad][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
+ PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(2+numBasis*(3*11+4))));
+} // _elasticityJacobian2D
+
+// ----------------------------------------------------------------------
+// Integrate elasticity term in Jacobian for 3-D cells.
+void
+pylith::feassemble::IntegratorElasticity::_elasticityJacobian3D(
+ const std::vector<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);
+
+ // 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][ 0];
+ const double C1122 = elasticConsts[iQuad][ 1];
+ const double C1133 = elasticConsts[iQuad][ 2];
+ const double C1112 = elasticConsts[iQuad][ 3]/2.0;
+ const double C1123 = elasticConsts[iQuad][ 4]/2.0;
+ const double C1113 = elasticConsts[iQuad][ 5]/2.0;
+ const double C2222 = elasticConsts[iQuad][ 6];
+ const double C2233 = elasticConsts[iQuad][ 7];
+ const double C2212 = elasticConsts[iQuad][ 8]/2.0;
+ const double C2223 = elasticConsts[iQuad][ 9]/2.0;
+ const double C2213 = elasticConsts[iQuad][10]/2.0;
+ const double C3333 = elasticConsts[iQuad][11];
+ const double C3312 = elasticConsts[iQuad][12]/2.0;
+ const double C3323 = elasticConsts[iQuad][13]/2.0;
+ const double C3313 = elasticConsts[iQuad][14]/2.0;
+ const double C1212 = elasticConsts[iQuad][15]/2.0;
+ const double C1223 = elasticConsts[iQuad][16]/2.0;
+ const double C1213 = elasticConsts[iQuad][17]/2.0;
+ const double C2323 = elasticConsts[iQuad][18]/2.0;
+ const double C2313 = elasticConsts[iQuad][19]/2.0;
+ const double C1313 = elasticConsts[iQuad][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
+ PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(3+numBasis*(6*26+9))));
+} // _elasticityJacobian3D
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2007-11-16 23:40:38 UTC (rev 8299)
@@ -0,0 +1,137 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/feassemble/IntegratorElasticity.hh
+ *
+ * @brief Object containing general elasticity operations for implicit
+ * and explicit time integration of the elasticity equation.
+ */
+
+#if !defined(pylith_feassemble_integratorelasticity_hh)
+#define pylith_feassemble_integratorelasticity_hh
+
+#include "Integrator.hh" // ISA Integrator
+#include "pylith/utils/array.hh" // USES std::vector, double_array
+
+namespace pylith {
+ namespace feassemble {
+ class IntegratorElasticity;
+ class TestIntegratorElasticity;
+ } // feassemble
+
+ namespace materials {
+ class ElasticMaterial;
+ } // feassemble
+} // pylith
+
+class pylith::feassemble::IntegratorElasticity : public Integrator
+{ // IntegratorElasticity
+ friend class TestIntegratorElasticity; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ IntegratorElasticity(void);
+
+ /// Destructor
+ ~IntegratorElasticity(void);
+
+ /** Set material.
+ *
+ * @param m Elastic material.
+ */
+ void material(materials::ElasticMaterial* m);
+
+ /** 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 field Current solution field.
+ * @param mesh Finite-element mesh
+ */
+ void updateState(const double t,
+ const ALE::Obj<real_section_type>& field,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Verify configuration is acceptable.
+ *
+ * @param mesh Finite-element mesh
+ */
+ void verifyConfiguration(const ALE::Obj<Mesh>& mesh);
+
+// PROTECTED METHODS ////////////////////////////////////////////////////
+protected :
+
+ /** Integrate elasticity term in residual for 1-D cells.
+ *
+ * @param stress Stress tensor for cell at quadrature points.
+ */
+ void _elasticityResidual1D(const std::vector<double_array>& stress);
+
+ /** Integrate elasticity term in residual for 2-D cells.
+ *
+ * @param stress Stress tensor for cell at quadrature points.
+ */
+ void _elasticityResidual2D(const std::vector<double_array>& stress);
+
+ /** Integrate elasticity term in residual for 3-D cells.
+ *
+ * @param stress Stress tensor for cell at quadrature points.
+ */
+ void _elasticityResidual3D(const std::vector<double_array>& stress);
+
+ /** Integrate elasticity term in Jacobian for 1-D cells.
+ *
+ * @param elasticConsts Matrix of elasticity constants at quadrature points.
+ */
+ void _elasticityJacobian1D(const std::vector<double_array>& elasticConsts);
+
+ /** Integrate elasticity term in Jacobian for 2-D cells.
+ *
+ * @param elasticConsts Matrix of elasticity constants at quadrature points.
+ */
+ void _elasticityJacobian2D(const std::vector<double_array>& elasticConsts);
+
+ /** Integrate elasticity term in Jacobian for 3-D cells.
+ *
+ * @param elasticConsts Matrix of elasticity constants at quadrature points.
+ */
+ void _elasticityJacobian3D(const std::vector<double_array>& elasticConsts);
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /// Not implemented.
+ IntegratorElasticity(const IntegratorElasticity& i);
+
+ /// Not implemented
+ const IntegratorElasticity& operator=(const IntegratorElasticity&);
+
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
+
+ /// Elastic material associated with integrator
+ materials::ElasticMaterial* _material;
+
+}; // IntegratorElasticity
+
+#endif // pylith_feassemble_integratorelasticity_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2007-11-16 23:40:38 UTC (rev 8299)
@@ -22,6 +22,7 @@
ElasticityImplicit.hh \
Integrator.hh \
Integrator.icc \
+ IntegratorElasticity.hh \
Elasticity.hh \
GeometryPoint1D.hh \
GeometryPoint2D.hh \
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2007-11-16 23:40:38 UTC (rev 8299)
@@ -131,8 +131,16 @@
materials::ElasticIsotropic3D material;
integrator.material(&material);
CPPUNIT_ASSERT_EQUAL(false, integrator._useSolnIncr);
- integrator.useSolnIncr(true);
- CPPUNIT_ASSERT_EQUAL(true, integrator._useSolnIncr);
+ try {
+ integrator.useSolnIncr(true);
+
+ // Should have thrown exception, so don't make it here.
+ CPPUNIT_ASSERT(false);
+ } catch (const std::logic_error& err) {
+ // Expect logic error so don't do anything.
+ } catch (...) {
+ CPPUNIT_ASSERT(false);
+ } // try/catch
} // testUseSolnIncr
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py 2007-11-16 23:40:38 UTC (rev 8299)
@@ -120,7 +120,11 @@
Test useSolnIncr().
"""
(mesh, integrator, fields) = self._initialize()
- integrator.useSolnIncr(True)
+ try:
+ integrator.useSolnIncr(True)
+ self.failIf(True)
+ except:
+ self.failIf(False)
return
More information about the cig-commits
mailing list