[cig-commits] r6641 - in short/3D/PyLith/trunk: . libsrc
libsrc/feassemble
brad at geodynamics.org
brad at geodynamics.org
Mon Apr 23 21:12:43 PDT 2007
Author: brad
Date: 2007-04-23 21:12:43 -0700 (Mon, 23 Apr 2007)
New Revision: 6641
Added:
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
Log:
Added IntegratorElasticity to factor out computation of total strain from IntegratorImplicit and IntegratorExplicit.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-04-24 04:09:08 UTC (rev 6640)
+++ short/3D/PyLith/trunk/TODO 2007-04-24 04:12:43 UTC (rev 6641)
@@ -28,6 +28,8 @@
a. Double check loops for calcConstant() and calcJacobian()
b. Create unit test (construction of constant term and Jacobian)
+ c. Add unit test for IntegratorElasticity::calcTotalStrain
+
2. Implement HDF5 output
3. Implement PML boundary conditions
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2007-04-24 04:09:08 UTC (rev 6640)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2007-04-24 04:12:43 UTC (rev 6641)
@@ -30,6 +30,7 @@
faults/SlipTimeFn.cc \
feassemble/ExplicitElasticity.cc \
feassemble/Integrator.cc \
+ feassemble/IntegratorElasticity.cc \
feassemble/IntegratorExplicit.cc \
feassemble/Quadrature.cc \
feassemble/Quadrature1D.cc \
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc 2007-04-24 04:09:08 UTC (rev 6640)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc 2007-04-24 04:12:43 UTC (rev 6641)
@@ -15,6 +15,7 @@
#include "ExplicitElasticity.hh" // implementation of class methods
#include "Quadrature.hh" // USES Quadrature
+#include "IntegratorElasticity.hh" // USES IntegratorElasticity
#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
#include "pylith/utils/array.hh" // USES double_array
@@ -174,10 +175,9 @@
// Compute action for elastic terms
if (1 == cellDim) {
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- for (int iBasis=0; iBasis < numBasis; ++iBasis)
- totalStrain[iQuad][0] +=
- basisDeriv[iQuad*numBasis+iBasis] * dispTCell[iBasis];
+ // Compute stresses
+ IntegratorElasticity::calcTotalStrain(&totalStrain, basisDeriv,
+ dispTCell, cellDim, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -196,18 +196,9 @@
throw std::runtime_error("Logging PETSc flops failed.");
} else if (2 == cellDim) {
- // Compute total strains
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- totalStrain[iQuad][0] +=
- basisDeriv[iQ+iBasis ] * dispTCell[iBasis ];
- totalStrain[iQuad][1] +=
- basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+1];
- totalStrain[iQuad][2] +=
- 0.5 * (basisDeriv[iQ+iBasis+1] * dispTCell[iBasis ] +
- basisDeriv[iQ+iBasis ] * dispTCell[iBasis+1]);
- } // for
- } // for
+ // Compute stresses
+ IntegratorElasticity::calcTotalStrain(&totalStrain, basisDeriv,
+ dispTCell, cellDim, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -230,26 +221,9 @@
throw std::runtime_error("Logging PETSc flops failed.");
} else if (3 == cellDim) {
- // Compute total strains
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- totalStrain[iQuad][0] +=
- basisDeriv[iQ+iBasis ] * dispTCell[iBasis ];
- totalStrain[iQuad][1] +=
- basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+1];
- totalStrain[iQuad][2] +=
- basisDeriv[iQ+iBasis+2] * dispTCell[iBasis+2];
- totalStrain[iQuad][3] +=
- 0.5 * (basisDeriv[iQ+iBasis+1] * dispTCell[iBasis ] +
- basisDeriv[iQ+iBasis ] * dispTCell[iBasis+1]);
- totalStrain[iQuad][4] +=
- 0.5 * (basisDeriv[iQ+iBasis+2] * dispTCell[iBasis+1] +
- basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+2]);
- totalStrain[iQuad][5] +=
- 0.5 * (basisDeriv[iQ+iBasis+2] * dispTCell[iBasis ] +
- basisDeriv[iQ+iBasis ] * dispTCell[iBasis+2]);
- } // for
- } // for
+ // Compute stresses
+ IntegratorElasticity::calcTotalStrain(&totalStrain, basisDeriv,
+ dispTCell, cellDim, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
Added: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2007-04-24 04:09:08 UTC (rev 6640)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2007-04-24 04:12:43 UTC (rev 6641)
@@ -0,0 +1,82 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "IntegratorElasticity.hh" // implementation of class methods
+
+#include <assert.h> // USES assert()
+#include <stdexcept> // USES std::logic_error
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::IntegratorElasticity::calcTotalStrain(
+ std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double* disp,
+ const int dimension,
+ const int numBasis)
+{ // calcTotalStrain
+ assert(0 != strain);
+
+ const int numQuadPts = strain->size();
+
+ assert(basisDeriv.size() == numQuadPts*numBasis*dimension);
+
+ if (3 == dimension) {
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ assert(6 == (*strain)[iQuad].size());
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ strain[iQuad][0] +=
+ basisDeriv[iQ+iBasis ] * disp[iBasis ];
+ strain[iQuad][1] +=
+ basisDeriv[iQ+iBasis+1] * disp[iBasis+1];
+ strain[iQuad][2] +=
+ basisDeriv[iQ+iBasis+2] * disp[iBasis+2];
+ strain[iQuad][3] +=
+ 0.5 * (basisDeriv[iQ+iBasis+1] * disp[iBasis ] +
+ basisDeriv[iQ+iBasis ] * disp[iBasis+1]);
+ strain[iQuad][4] +=
+ 0.5 * (basisDeriv[iQ+iBasis+2] * disp[iBasis+1] +
+ basisDeriv[iQ+iBasis+1] * disp[iBasis+2]);
+ strain[iQuad][5] +=
+ 0.5 * (basisDeriv[iQ+iBasis+2] * disp[iBasis ] +
+ basisDeriv[iQ+iBasis ] * disp[iBasis+2]);
+ } // for
+ } // for
+ } else if (2 == dimension) {
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ assert(3 == (*strain)[iQuad].size());
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ strain[iQuad][0] +=
+ basisDeriv[iQ+iBasis ] * disp[iBasis ];
+ strain[iQuad][1] +=
+ basisDeriv[iQ+iBasis+1] * disp[iBasis+1];
+ strain[iQuad][2] +=
+ 0.5 * (basisDeriv[iQ+iBasis+1] * disp[iBasis ] +
+ basisDeriv[iQ+iBasis ] * disp[iBasis+1]);
+ } // for
+ } // for
+ } else if (1 == dimension) {
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ assert(1 == (*strain)[iQuad].size());
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ (*strain)[iQuad][0] +=
+ basisDeriv[iQuad*numBasis+iBasis] * disp[iBasis];
+ } // for
+ } else {
+ throw std::logic_error("Dimension out of range.");
+ } // else
+} // calcTotalStrain
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2007-04-24 04:09:08 UTC (rev 6640)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2007-04-24 04:12:43 UTC (rev 6641)
@@ -0,0 +1,59 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/feassemble/IntegratorElasticity.hh
+ *
+ * @brief C++ Utility class for general operations associated with
+ * integrating finite-element actions for solid elasticity.
+ */
+
+#if !defined(pylith_feassemble_integratorelasticity_hh)
+#define pylith_feassemble_integratorelasticity_hh
+
+#include "pylith/utils/array.hh" // USES double_array, std::vector
+
+namespace pylith {
+ namespace feassemble {
+ class IntegratorElasticity;
+ class TestIntegratorElasticity;
+ } // feassemble
+} // pylith
+
+class pylith::feassemble::IntegratorElasticity
+{ // IntegratorElasticity
+ friend class TestIntegratorElasticity; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Compute total strain in at quadrature points of a cell.
+ *
+ * @param strain Strain tensor at quadrature points.
+ * @param basisDeriv Derivatives of basis functions at quadrature points.
+ * @param disp Displacement at vertices of cell.
+ * @param dimension Dimension of cell.
+ * @param numBasis Number of basis functions for cell.
+ */
+
+ static
+ void calcTotalStrain(std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double* disp,
+ const int dimension,
+ const int numBasis);
+
+}; // 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-04-24 04:09:08 UTC (rev 6640)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2007-04-24 04:12:43 UTC (rev 6641)
@@ -17,6 +17,7 @@
ExplicitElasticity.hh \
ExplicitElasticity.icc \
Integrator.hh \
+ IntegratorElasticity.hh \
IntegratorExplicit.hh \
Quadrature.hh \
Quadrature.icc \
More information about the cig-commits
mailing list