[cig-commits] r8262 - short/3D/PyLith/trunk/libsrc/feassemble
brad at geodynamics.org
brad at geodynamics.org
Fri Nov 9 15:48:53 PST 2007
Author: brad
Date: 2007-11-09 15:48:53 -0800 (Fri, 09 Nov 2007)
New Revision: 8262
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
Log:
Started debugging in earnest explicit time stepping. Worked on incorporating all of the optimizations from ElasticityImplicit. Need to refactor elasticity functions if possible.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2007-11-09 23:47:43 UTC (rev 8261)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2007-11-09 23:48:53 UTC (rev 8262)
@@ -20,6 +20,7 @@
#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
#include "pylith/utils/array.hh" // USES double_array
+#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
#include "petscmat.h" // USES PetscMat
#include "spatialdata/spatialdb/SpatialDB.hh"
@@ -98,6 +99,10 @@
topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh)
{ // integrateResidual
+ /// Member prototype for _elasticityResidualXD()
+ typedef void (pylith::feassemble::ElasticityExplicit::*elasticityResidual_fn_type)
+ (const std::vector<double_array>&);
+
assert(0 != _quadrature);
assert(0 != _material);
assert(!residual.isNull());
@@ -106,6 +111,29 @@
PetscErrorCode err = 0;
+ // Set variables dependent on dimension of cell
+ const int cellDim = _quadrature->cellDim();
+ int tensorSize = 0;
+ Elasticity::totalStrain_fn_type calcTotalStrainFn;
+ elasticityResidual_fn_type elasticityResidualFn;
+ if (1 == cellDim) {
+ tensorSize = 1;
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityExplicit::_elasticityResidual1D;
+ calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain1D;
+ } else if (2 == cellDim) {
+ tensorSize = 3;
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityExplicit::_elasticityResidual2D;
+ calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain2D;
+ } else if (3 == cellDim) {
+ tensorSize = 6;
+ elasticityResidualFn =
+ &pylith::feassemble::ElasticityExplicit::_elasticityResidual3D;
+ 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());
@@ -117,10 +145,10 @@
mesh->getRealSection("coordinates");
assert(!coordinates.isNull());
const ALE::Obj<real_section_type>& dispTpdt = fields->getHistoryItem(0);
+ assert(!dispTpdt.isNull());
const ALE::Obj<real_section_type>& dispT = fields->getHistoryItem(1);
+ assert(!dispT.isNull());
const ALE::Obj<real_section_type>& dispTmdt = fields->getHistoryItem(2);
- assert(!dispTpdt.isNull());
- assert(!dispT.isNull());
assert(!dispTmdt.isNull());
// Get parameters used in integration.
@@ -134,7 +162,6 @@
assert(quadWts.size() == numQuadPts);
const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
- const int cellDim = _quadrature->cellDim();
/** :TODO:
*
@@ -146,6 +173,9 @@
if (cellDim != spaceDim)
throw std::logic_error("Not implemented yet.");
+ // Precompute the geometric and function space information
+ _quadrature->precomputeGeometry(mesh, coordinates, cells);
+
// Allocate vectors for cell values.
_initCellVector();
const int cellVecSize = numBasis*spaceDim;
@@ -154,28 +184,18 @@
double_array dispTmdtCell(cellVecSize);
// Allocate vector for total strain
- int tensorSize = 0;
- if (1 == cellDim)
- tensorSize = 1;
- else if (2 == cellDim)
- tensorSize = 3;
- else if (3 == cellDim)
- tensorSize = 6;
- else {
- std::cerr << "Unknown case for cellDim '" << cellDim << "'." << std::endl;
- assert(0);
- } // else
std::vector<double_array> totalStrain(numQuadPts);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
totalStrain[iQuad].resize(tensorSize);
totalStrain[iQuad] = 0.0;
} // for
+ int c = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
++c_iter) {
// Compute geometry information for current cell
- _quadrature->computeGeometry(mesh, coordinates, *c_iter);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
// Get state variables for cell.
_material->getStateVarsCell(*c_iter, numQuadPts);
@@ -184,6 +204,8 @@
_resetCellVector();
// Restrict input fields to cell
+ // :TODO: Replace with optimized restricts().
+ // NEED to add optimization tags to ElasticityExplicit object.
mesh->restrict(dispTpdt, *c_iter, &dispTpdtCell[0], cellVecSize);
mesh->restrict(dispT, *c_iter, &dispTCell[0], cellVecSize);
mesh->restrict(dispTmdt, *c_iter, &dispTmdtCell[0], cellVecSize);
@@ -213,84 +235,12 @@
} // for
PetscLogFlopsNoCheck(numQuadPts*(3+numBasis*(1+numBasis*(6*spaceDim))));
- // Compute action for elastic terms
- if (1 == cellDim) {
- // Compute stresses
- Elasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
- const std::vector<double_array>& stress =
- _material->calcStress(totalStrain);
+ // Compute B(transpose) * sigma, first computing strains
+ calcTotalStrainFn(&totalStrain, basisDeriv, dispTCell, numBasis);
+ const std::vector<double_array>& stress =
+ _material->calcStress(totalStrain);
+ CALL_MEMBER_FN(*this, elasticityResidualFn)(stress);
- // Compute elastic action
- 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));
-
- } else if (2 == cellDim) {
- // Compute stresses
- Elasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
- const std::vector<double_array>& stress =
- _material->calcStress(totalStrain);
-
- // Compute elastic action
- 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)));
-
- } else if (3 == cellDim) {
- // Compute stresses
- Elasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
- const std::vector<double_array>& stress =
- _material->calcStress(totalStrain);
-
- // Compute elastic action
- 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)));
- } else {
- std::cerr << "Unknown case for cellDim '" << cellDim << "'."
- << std::endl;
- assert(0);
- } // if/else
// Assemble cell contribution into field
mesh->updateAdd(residual, *c_iter, _cellVector);
} // for
@@ -334,15 +284,20 @@
const double_array& quadWts = _quadrature->quadWts();
const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ // Precompute the geometric and function space information
+ _quadrature->precomputeGeometry(mesh, coordinates, cells);
+
// Allocate vector for cell values (if necessary)
_initCellMatrix();
+ int c = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
++c_iter) {
// Compute geometry information for current cell
- _quadrature->computeGeometry(mesh, coordinates, *c_iter);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
// Get state variables for cell.
_material->getStateVarsCell(*c_iter, numQuadPts);
@@ -398,6 +353,7 @@
const ALE::Obj<real_section_type>& disp,
const ALE::Obj<Mesh>& mesh)
{ // updateState
+ assert(0 != _quadrature);
assert(0 != _material);
assert(!disp.isNull());
@@ -405,6 +361,22 @@
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());
@@ -420,34 +392,23 @@
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
- const int cellDim = _quadrature->cellDim();
const int cellVecSize = numBasis*spaceDim;
double_array dispCell(cellVecSize);
// Allocate vector for total strain
- int tensorSize = 0;
- if (1 == cellDim)
- tensorSize = 1;
- else if (2 == cellDim)
- tensorSize = 3;
- else if (3 == cellDim)
- tensorSize = 6;
- else {
- std::cerr << "Unknown case for cellDim '" << cellDim << "'." << std::endl;
- assert(0);
- } // else
std::vector<double_array> totalStrain(numQuadPts);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
totalStrain[iQuad].resize(tensorSize);
totalStrain[iQuad] = 0.0;
} // for
+ int c = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
++c_iter) {
// Compute geometry information for current cell
- _quadrature->computeGeometry(mesh, coordinates, *c_iter);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
// Restrict input fields to cell
mesh->restrict(disp, *c_iter, &dispCell[0], cellVecSize);
@@ -455,21 +416,10 @@
// Get cell geometry information that depends on cell
const double_array& basisDeriv = _quadrature->basisDeriv();
- // Compute action for elastic terms
- if (1 == cellDim)
- Elasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
- dispCell, numBasis);
- else if (2 == cellDim)
- Elasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
- dispCell, numBasis);
- else if (3 == cellDim)
- Elasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
- dispCell, numBasis);
- else {
- std::cerr << "Unknown case for cellDim '" << cellDim << "'."
- << std::endl;
- assert(0);
- } // else
+ // Compute strains
+ calcTotalStrainFn(&totalStrain, basisDeriv, dispCell, numBasis);
+
+ // Update material state
_material->updateState(totalStrain, *c_iter);
} // for
} // updateState
@@ -529,5 +479,317 @@
} // 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-09 23:47:43 UTC (rev 8261)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2007-11-09 23:48:53 UTC (rev 8262)
@@ -56,6 +56,7 @@
#define pylith_feassemble_elasticityexplicit_hh
#include "Integrator.hh" // ISA Integrator
+#include "pylith/utils/array.hh" // USES std::vector, double_array
namespace pylith {
namespace feassemble {
@@ -159,6 +160,45 @@
// 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);
More information about the cig-commits
mailing list