[cig-commits] r15945 - short/3D/PyLith/branches/pylith-friction/libsrc/feassemble
brad at geodynamics.org
brad at geodynamics.org
Mon Nov 9 17:22:39 PST 2009
Author: brad
Date: 2009-11-09 17:22:38 -0800 (Mon, 09 Nov 2009)
New Revision: 15945
Modified:
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.hh
Log:
Fixed calcTotalStrain for large deformations.
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc 2009-11-10 00:21:51 UTC (rev 15944)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc 2009-11-10 01:22:38 UTC (rev 15945)
@@ -55,9 +55,7 @@
bool
pylith::feassemble::IntegratorElasticityLgDeform::needNewJacobian(void)
{ // needNewJacobian
- assert(0 != _material);
- if (!_needNewJacobian)
- _needNewJacobian = _material->needNewJacobian();
+ _needNewJacobian = IntegratorElasticity::needNewJacobian();
return _needNewJacobian;
} // needNewJacobian
@@ -101,6 +99,8 @@
// Allocate arrays for cell data.
double_array dispCell(numBasis*spaceDim);
double_array strainCell(numQuadPts*tensorSize);
+ double_array deformCell(numQuadPts*spaceDim*spaceDim);
+ deformCell = 0.0;
strainCell = 0.0;
// Get cell information
@@ -120,7 +120,6 @@
topology::Mesh::RestrictVisitor dispVisitor(*dispSection,
dispCell.size(), &dispCell[0]);
-#if !defined(PRECOMPUTE_GEOMETRY)
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
@@ -128,7 +127,6 @@
topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
coordinatesCell.size(),
&coordinatesCell[0]);
-#endif
// Loop over cells
for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
@@ -150,9 +148,12 @@
// Get cell geometry information that depends on cell
const double_array& basisDeriv = _quadrature->basisDeriv();
+ // Compute deformation tensor.
+ _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispCell,
+ numBasis, numQuadPts, spaceDim);
+
// Compute strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispCell,
- numBasis, numQuadPts);
+ calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
// Update material state
_material->updateStateVars(strainCell, *c_iter);
@@ -196,6 +197,7 @@
// Allocate arrays for cell data.
double_array dispCell(numBasis*spaceDim);
+ double_array deformCell(numQuadPts*spaceDim*spaceDim);
double_array strainCell(numQuadPts*tensorSize);
strainCell = 0.0;
double_array stressCell(numQuadPts*tensorSize);
@@ -218,7 +220,6 @@
topology::Mesh::RestrictVisitor dispVisitor(*dispSection,
dispCell.size(), &dispCell[0]);
-#if !defined(PRECOMPUTE_GEOMETRY)
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
@@ -226,7 +227,6 @@
topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
coordinatesCell.size(),
&coordinatesCell[0]);
-#endif
const ALE::Obj<RealSection>& fieldSection = field->section();
assert(!fieldSection.isNull());
@@ -251,10 +251,13 @@
// Get cell geometry information that depends on cell
const double_array& basisDeriv = _quadrature->basisDeriv();
+ // Compute deformation tensor.
+ _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispCell,
+ numBasis, numQuadPts, spaceDim);
+
// Compute strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispCell,
- numBasis, numQuadPts);
-
+ calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
+
if (!calcStress)
fieldSection->updatePoint(*c_iter, &strainCell[0]);
else {
@@ -628,96 +631,114 @@
} // _elasticityJacobian3D
// ----------------------------------------------------------------------
-void
+// Calculate Green-Lagrange strain tensor at quadrature points of a 1-D cell.
+void
pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain1D(
- double_array* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis,
- const int numQuadPts)
-{ // calcTotalStrain1D
+ double_array* strain,
+ const double_array& deform,
+ const int numQuadPts)
+{ // _calcTotalStrain1D
+ // Green-Lagrange strain tensor = 1/2 ( X^T X - I )
+ // X: deformation tensor
+ // I: identity matrix
+
assert(0 != strain);
const int dim = 1;
-
- assert(basisDeriv.size() == numQuadPts*numBasis*dim);
- assert(disp.size() == numBasis*dim);
+ const int strainSize = 1;
+ assert(deform.size() == numQuadPts*dim*dim);
+ assert(strain->size() == numQuadPts*strainSize);
- (*strain) = 0.0;
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- for (int iBasis=0; iBasis < numBasis; ++iBasis)
- (*strain)[iQuad] += basisDeriv[iQuad*numBasis+iBasis] * disp[iBasis];
- } // for
-} // calcTotalStrain1D
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ (*strain)[iQuad] = 0.5*(deform[iQuad]*deform[iQuad] - 1.0);
+} // _calcTotalStrain1D
+
// ----------------------------------------------------------------------
-void
+// Calculate Green-Lagrange strain tensor at quadrature points of a 2-D cell.
+void
pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain2D(
- double_array* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis,
- const int numQuadPts)
-{ // calcTotalStrain2D
+ double_array* strain,
+ const double_array& deform,
+ const int numQuadPts)
+{ // _calcTotalStrain2D
+ // Green-Lagrange strain tensor = 1/2 ( X^T X - I )
+ // X: deformation tensor
+ // I: identity matrix
+
assert(0 != strain);
-
+
const int dim = 2;
-
- assert(basisDeriv.size() == numQuadPts*numBasis*dim);
- assert(disp.size() == numBasis*dim);
-
- (*strain) = 0.0;
+ const int deformSize = dim*dim;
const int strainSize = 3;
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
- (*strain)[iQuad*strainSize+0] +=
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
- (*strain)[iQuad*strainSize+1] +=
- basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
- (*strain)[iQuad*strainSize+2] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
- } // for
-} // calcTotalStrain2D
+ assert(deform.size() == numQuadPts*deformSize);
+ assert(strain->size() == numQuadPts*strainSize);
+ for (int iQuad=0, iDeform=0, iStrain=0;
+ iQuad < numQuadPts;
+ ++iQuad, iDeform+=deformSize, iStrain+=strainSize) {
+ (*strain)[iStrain ] =
+ 0.5 * (deform[iDeform ]*deform[iDeform ] +
+ deform[iDeform+2]*deform[iDeform+2] - 1.0);
+ (*strain)[iStrain+1] =
+ 0.5 * (deform[iDeform+1]*deform[iDeform+1] +
+ deform[iDeform+3]*deform[iDeform+3] - 1.0);
+ (*strain)[iStrain+2] =
+ 0.5 * (deform[iDeform ]*deform[iDeform+2] +
+ deform[iDeform+1]*deform[iDeform+3]);
+ } // for
+} // _calcTotalStrain2D
+
// ----------------------------------------------------------------------
-void
+// Calculate Green-Lagrange strain tensor at quadrature points of a 3-D cell.
+void
pylith::feassemble::IntegratorElasticityLgDeform::_calcTotalStrain3D(
- double_array* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis,
- const int numQuadPts)
-{ // calcTotalStrain3D
+ double_array* strain,
+ const double_array& deform,
+ const int numQuadPts)
+{ // _calcTotalStrain3D
+ // Green-Lagrange strain tensor = 1/2 ( X^T X - I )
+ // X: deformation tensor
+ // I: identity matrix
+
assert(0 != strain);
const int dim = 3;
-
- assert(basisDeriv.size() == numQuadPts*numBasis*dim);
- assert(disp.size() == numBasis*dim);
-
- (*strain) = 0.0;
+ const int deformSize = dim*dim;
const int strainSize = 6;
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
- (*strain)[iQuad*strainSize+0] +=
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
- (*strain)[iQuad*strainSize+1] +=
- basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
- (*strain)[iQuad*strainSize+2] +=
- basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+2];
- (*strain)[iQuad*strainSize+3] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
- (*strain)[iQuad*strainSize+4] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+1] +
- basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+2]);
- (*strain)[iQuad*strainSize+5] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim ] +
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+2]);
- } // for
-} // calcTotalStrain3D
+ assert(deform.size() == numQuadPts*dim*dim);
+ assert(strain->size() == numQuadPts*strainSize);
+ for (int iQuad=0, iDeform=0, iStrain=0;
+ iQuad < numQuadPts;
+ ++iQuad, iDeform+=deformSize, iStrain+=strainSize) {
+ (*strain)[iStrain ] =
+ 0.5 * (deform[iDeform ]*deform[iDeform ] +
+ deform[iDeform+3]*deform[iDeform+3] +
+ deform[iDeform+6]*deform[iDeform+6] - 1.0);
+ (*strain)[iStrain+1] =
+ 0.5 * (deform[iDeform+1]*deform[iDeform+1] +
+ deform[iDeform+4]*deform[iDeform+4] +
+ deform[iDeform+7]*deform[iDeform+7] - 1.0);
+ (*strain)[iStrain+2] =
+ 0.5 * (deform[iDeform+2]*deform[iDeform+2] +
+ deform[iDeform+5]*deform[iDeform+5] +
+ deform[iDeform+8]*deform[iDeform+8] - 1.0);
+ (*strain)[iStrain+3] =
+ 0.5 * (deform[iDeform ]*deform[iDeform+1] +
+ deform[iDeform+3]*deform[iDeform+4] +
+ deform[iDeform+6]*deform[iDeform+7]);
+ (*strain)[iStrain+4] =
+ 0.5 * (deform[iDeform+1]*deform[iDeform+2] +
+ deform[iDeform+4]*deform[iDeform+5] +
+ deform[iDeform+7]*deform[iDeform+8]);
+ (*strain)[iStrain+5] =
+ 0.5 * (deform[iDeform+0]*deform[iDeform+2] +
+ deform[iDeform+3]*deform[iDeform+5] +
+ deform[iDeform+6]*deform[iDeform+8]);
+ } // for
+} // _calcTotalStrain3D
+
// ----------------------------------------------------------------------
// Calculate deformation tensor.
void
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.hh 2009-11-10 00:21:51 UTC (rev 15944)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.hh 2009-11-10 01:22:38 UTC (rev 15945)
@@ -39,8 +39,6 @@
typedef void (*totalStrain_fn_type)(double_array*,
const double_array&,
- const double_array&,
- const int,
const int);
@@ -126,49 +124,40 @@
*/
void _elasticityJacobian3D(const double_array& elasticConsts);
- /** Compute total strain in at quadrature points of a cell.
+ /** Calculate Green-Lagrange strain tensor at quadrature points of a
+ * 1-D cell.
*
- * @param strain Strain tensor at quadrature points.
- * @param basisDeriv Derivatives of basis functions at quadrature points.
- * @param disp Displacement at vertices of cell.
- * @param numBasis Number of basis functions for cell.
+ * @param strain Green-Lagrange strain tensor at quadrature points (output).
+ * @param deform Deformation tensor at quadrature points.
* @param numQuadPts Number of quadrature points.
*/
static
void _calcTotalStrain1D(double_array* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis,
+ const double_array& deform,
const int numQuadPts);
- /** Compute total strain in at quadrature points of a cell.
+ /** Calculate Green-Lagrange strain tensor at quadrature points of a
+ * 2-D cell.
*
- * @param strain Strain tensor at quadrature points.
- * @param basisDeriv Derivatives of basis functions at quadrature points.
- * @param disp Displacement at vertices of cell.
- * @param numBasis Number of basis functions for cell.
+ * @param strain Green-Lagrange strain tensor at quadrature points (output).
+ * @param deform Deformation tensor at quadrature points.
* @param numQuadPts Number of quadrature points.
*/
static
void _calcTotalStrain2D(double_array* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis,
+ const double_array& deform,
const int numQuadPts);
- /** Compute total strain in at quadrature points of a cell.
+ /** Calculate Green-Lagrange strain tensor at quadrature points of a
+ * 3-D cell.
*
- * @param strain Strain tensor at quadrature points.
- * @param basisDeriv Derivatives of basis functions at quadrature points.
- * @param disp Displacement at vertices of cell.
- * @param numBasis Number of basis functions for cell.
+ * @param strain Green-Lagrange strain tensor at quadrature points (output).
+ * @param deform Deformation tensor at quadrature points.
* @param numQuadPts Number of quadrature points.
*/
static
void _calcTotalStrain3D(double_array* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis,
+ const double_array& deform,
const int numQuadPts);
/** Calculate deformation tensor.
More information about the CIG-COMMITS
mailing list