[cig-commits] r8300 - in short/3D/PyLith/trunk: libsrc
libsrc/feassemble unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Fri Nov 16 17:31:00 PST 2007
Author: brad
Date: 2007-11-16 17:30:59 -0800 (Fri, 16 Nov 2007)
New Revision: 8300
Added:
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.hh
Removed:
short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.hh
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.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/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
Log:
Moved Elasticity member functions to IntegratorElasticity. Removed Elasticity object. Changed name of TestElasticity to TestIntegratorElasticity (files and source code).
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2007-11-17 01:30:59 UTC (rev 8300)
@@ -36,7 +36,6 @@
faults/SlipTimeFn.cc \
feassemble/CellGeometry.cc \
feassemble/Constraint.cc \
- feassemble/Elasticity.cc \
feassemble/ElasticityExplicit.cc \
feassemble/ElasticityImplicit.cc \
feassemble/GeometryPoint1D.cc \
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.cc 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.cc 2007-11-17 01:30:59 UTC (rev 8300)
@@ -1,110 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-#include <portinfo>
-
-#include "Elasticity.hh" // implementation of class methods
-
-#include <assert.h> // USES assert()
-
-// ----------------------------------------------------------------------
-void
-pylith::feassemble::Elasticity::calcTotalStrain1D(
- std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis)
-{ // calcTotalStrain1D
- assert(0 != strain);
-
- const int dim = 1;
- const int numQuadPts = strain->size();
-
- assert(basisDeriv.size() == numQuadPts*numBasis*dim);
- assert(disp.size() == numBasis*dim);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- assert(1 == (*strain)[iQuad].size());
- (*strain)[iQuad] *= 0.0;
- for (int iBasis=0; iBasis < numBasis; ++iBasis)
- (*strain)[iQuad][0] +=
- basisDeriv[iQuad*numBasis+iBasis] * disp[iBasis];
- } // for
-} // calcTotalStrain1D
-
-// ----------------------------------------------------------------------
-void
-pylith::feassemble::Elasticity::calcTotalStrain2D(
- std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis)
-{ // calcTotalStrain2D
- assert(0 != strain);
-
- const int dim = 2;
- const int numQuadPts = strain->size();
-
- assert(basisDeriv.size() == numQuadPts*numBasis*dim);
- assert(disp.size() == numBasis*dim);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- assert(3 == (*strain)[iQuad].size());
- (*strain)[iQuad] *= 0.0;
- for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
- (*strain)[iQuad][0] += basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
- (*strain)[iQuad][1] += basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
- (*strain)[iQuad][2] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
- } // for
- } // for
-} // calcTotalStrain2D
-
-// ----------------------------------------------------------------------
-void
-pylith::feassemble::Elasticity::calcTotalStrain3D(
- std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis)
-{ // calcTotalStrain3D
- assert(0 != strain);
-
- const int dim = 3;
- const int numQuadPts = strain->size();
-
- assert(basisDeriv.size() == numQuadPts*numBasis*dim);
- assert(disp.size() == numBasis*dim);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- assert(6 == (*strain)[iQuad].size());
- (*strain)[iQuad] *= 0.0;
- for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
- (*strain)[iQuad][0] += basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
- (*strain)[iQuad][1] += basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
- (*strain)[iQuad][2] += basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+2];
- (*strain)[iQuad][3] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
- (*strain)[iQuad][4] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+1] +
- basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+2]);
- (*strain)[iQuad][5] +=
- 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim ] +
- basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+2]);
- } // for
- } // for
-} // calcTotalStrain3D
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.hh 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Elasticity.hh 2007-11-17 01:30:59 UTC (rev 8300)
@@ -1,95 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-/**
- * @file pylith/feassemble/Elasticity.hh
- *
- * @brief C++ Utility class for general operations associated with
- * integrating finite-element actions for solid elasticity.
- */
-
-#if !defined(pylith_feassemble_elasticity_hh)
-#define pylith_feassemble_elasticity_hh
-
-#include "pylith/utils/array.hh" // USES double_array, std::vector
-
-namespace pylith {
- namespace feassemble {
- class Elasticity;
- class TestElasticity;
- } // feassemble
-} // pylith
-
-class pylith::feassemble::Elasticity
-{ // Elasticity
- friend class TestElasticity; // unit testing
-
-// PUBLIC TYPEDEFS //////////////////////////////////////////////////////
-public :
-
- typedef void (*totalStrain_fn_type)(std::vector<double_array>*,
- const double_array&,
- const double_array&,
- const int);
-
-
-// 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 calcTotalStrain1D(std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis);
-
- /** Compute total strain in at quadrature points of a cell.
- *
- * @param strain Strain tensor at quadrature points.
- * @param basisDeriv Derivatives of basis functions at quadrature points.
- * @param disp Displacement at vertices of cell.
- * @param numBasis Number of basis functions for cell.
- */
-
- static
- void calcTotalStrain2D(std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis);
-
- /** Compute total strain in at quadrature points of a cell.
- *
- * @param strain Strain tensor at quadrature points.
- * @param basisDeriv Derivatives of basis functions at quadrature points.
- * @param disp Displacement at vertices of cell.
- * @param numBasis Number of basis functions for cell.
- */
-
- static
- void calcTotalStrain3D(std::vector<double_array>* strain,
- const double_array& basisDeriv,
- const double_array& disp,
- const int numBasis);
-
-}; // Elasticity
-
-#endif // pylith_feassemble_elasticity_hh
-
-// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2007-11-17 01:30:59 UTC (rev 8300)
@@ -15,12 +15,11 @@
#include "ElasticityExplicit.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/topology/FieldsManager.hh" // USES FieldsManager
-#include "pylith/utils/array.hh" // USES double_array
+#include "pylith/utils/array.hh" // USES double_array
#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
#include "petscmat.h" // USES PetscMat
@@ -94,23 +93,26 @@
// Set variables dependent on dimension of cell
const int cellDim = _quadrature->cellDim();
int tensorSize = 0;
- Elasticity::totalStrain_fn_type calcTotalStrainFn;
+ totalStrain_fn_type calcTotalStrainFn;
elasticityResidual_fn_type elasticityResidualFn;
if (1 == cellDim) {
tensorSize = 1;
elasticityResidualFn =
&pylith::feassemble::ElasticityExplicit::_elasticityResidual1D;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain1D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
} else if (2 == cellDim) {
tensorSize = 3;
elasticityResidualFn =
&pylith::feassemble::ElasticityExplicit::_elasticityResidual2D;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain2D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
} else if (3 == cellDim) {
tensorSize = 6;
elasticityResidualFn =
&pylith::feassemble::ElasticityExplicit::_elasticityResidual3D;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain3D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D;
} else
assert(0);
@@ -168,21 +170,21 @@
} // for
#ifdef FASTER
- if (_dispTTags.find(_material->id()) == _dispTTags.end()) {
- _dispTTags[_material->id()] =
+ if (_dispTags.find(_material->id()) == _dispTags.end()) {
+ _dispTags[_material->id()] =
mesh->calculateCustomAtlas(dispT, cells);
} // if
- const int dispTTag = _dispTTags[_material->id()];
+ const int dispTTag = _dispTags[_material->id()];
if (_dispTmdtTags.find(_material->id()) == _dispTmdtTags.end()) {
_dispTmdtTags[_material->id()] =
- dispTmdt->copyCustomAtlas(dispT, _dispTTags[_material->id()]);
+ dispTmdt->copyCustomAtlas(dispT, _dispTags[_material->id()]);
} // if
const int dispTmdtTag = _dispTmdtTags[_material->id()];
if (_residualTags.find(_material->id()) == _residualTags.end()) {
_residualTags[_material->id()] =
- residual->copyCustomAtlas(dispT, _dispTTags[_material->id()]);
+ residual->copyCustomAtlas(dispT, _dispTags[_material->id()]);
} // if
const int residualTag = _residualTags[_material->id()];
#endif
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.hh 2007-11-17 01:30:59 UTC (rev 8300)
@@ -125,7 +125,7 @@
topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh);
-// PRIVATE METHODS //////////////////////////////////////////////////////
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
/// Not implemented.
@@ -140,9 +140,7 @@
double _dtm1; ///< Time step for t-dt1 -> t
// Optimization
- std::map<int, int> _dispTTags; ///< Tags indexing dispT field
std::map<int, int> _dispTmdtTags; ///< Tags indexing dispTmdt field
- std::map<int, int> _residualTags; ///< tags indexing residual field
}; // ElasticityExplicit
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-11-17 01:30:59 UTC (rev 8300)
@@ -15,7 +15,6 @@
#include "ElasticityImplicit.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
@@ -112,23 +111,26 @@
// Set variables dependent on dimension of cell
const int cellDim = _quadrature->cellDim();
int tensorSize = 0;
- Elasticity::totalStrain_fn_type calcTotalStrainFn;
+ totalStrain_fn_type calcTotalStrainFn;
elasticityResidual_fn_type elasticityResidualFn;
if (1 == cellDim) {
tensorSize = 1;
elasticityResidualFn =
&pylith::feassemble::ElasticityImplicit::_elasticityResidual1D;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain1D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
} else if (2 == cellDim) {
tensorSize = 3;
elasticityResidualFn =
&pylith::feassemble::ElasticityImplicit::_elasticityResidual2D;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain2D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
} else if (3 == cellDim) {
tensorSize = 6;
elasticityResidualFn =
&pylith::feassemble::ElasticityImplicit::_elasticityResidual3D;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain3D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D;
} else
assert(0);
@@ -154,16 +156,17 @@
const int spaceDim = _quadrature->spaceDim();
#ifdef FASTER
- int c_index = 0;
-
- if (_dTags.find(_material->id()) == _dTags.end()) {
- _dTags[_material->id()] = mesh->calculateCustomAtlas(dispTBctpdt, cells);
- }
- if (_rTags.find(_material->id()) == _rTags.end()) {
- _rTags[_material->id()] = residual->copyCustomAtlas(dispTBctpdt, _dTags[_material->id()]);
- }
- const int dTag = _dTags[_material->id()];
- const int rTag = _rTags[_material->id()];
+ if (_dispTags.find(_material->id()) == _dispTags.end()) {
+ _dispTags[_material->id()] =
+ mesh->calculateCustomAtlas(dispTBctpdt, cells);
+ } // if
+ const int dispTBctpdtTag = _dispTags[_material->id()];
+
+ if (_residualTags.find(_material->id()) == _residualTags.end()) {
+ _residualTags[_material->id()] =
+ residual->copyCustomAtlas(dispTBctpdt, _dispTags[_material->id()]);
+ } // if
+ const int residualTag = _residualTags[_material->id()];
#endif
// Precompute the geometric and function space information
_quadrature->precomputeGeometry(mesh, coordinates, cells);
@@ -183,6 +186,7 @@
PetscLogEventEnd(setupEvent,0,0,0,0);
// Loop over cells
+ int c_index = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
++c_iter, ++c_index) {
@@ -202,7 +206,7 @@
// Restrict input fields to cell
PetscLogEventBegin(restrictEvent,0,0,0,0);
#ifdef FASTER
- mesh->restrict(dispTBctpdt, dTag, c_index, &dispTBctpdtCell[0],
+ mesh->restrict(dispTBctpdt, dispTBctpdtTag, c_index, &dispTBctpdtCell[0],
cellVecSize);
#else
mesh->restrict(dispTBctpdt, *c_iter, &dispTBctpdtCell[0], cellVecSize);
@@ -261,7 +265,7 @@
// Assemble cell contribution into field
PetscLogEventBegin(updateEvent,0,0,0,0);
#ifdef FASTER
- mesh->updateAdd(residual, rTag, c_index, _cellVector);
+ mesh->updateAdd(residual, residualTag, c_index, _cellVector);
#else
mesh->updateAdd(residual, *c_iter, _cellVector);
#endif
@@ -292,23 +296,26 @@
// Set variables dependent on dimension of cell
const int cellDim = _quadrature->cellDim();
int tensorSize = 0;
- Elasticity::totalStrain_fn_type calcTotalStrainFn;
+ totalStrain_fn_type calcTotalStrainFn;
elasticityJacobian_fn_type elasticityJacobianFn;
if (1 == cellDim) {
tensorSize = 1;
elasticityJacobianFn =
&pylith::feassemble::ElasticityImplicit::_elasticityJacobian1D;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain1D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
} else if (2 == cellDim) {
tensorSize = 3;
elasticityJacobianFn =
&pylith::feassemble::ElasticityImplicit::_elasticityJacobian2D;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain2D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
} else if (3 == cellDim) {
tensorSize = 6;
elasticityJacobianFn =
&pylith::feassemble::ElasticityImplicit::_elasticityJacobian3D;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain3D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D;
} else
assert(0);
@@ -357,10 +364,11 @@
} // for
#ifdef FASTER
- if (_dTags.find(_material->id()) == _dTags.end()) {
- _dTags[_material->id()] = mesh->calculateCustomAtlas(dispTBctpdt, cells);
- }
- const int dTag = _dTags[_material->id()];
+ if (_dispTags.find(_material->id()) == _dispTags.end()) {
+ _dispTags[_material->id()] =
+ mesh->calculateCustomAtlas(dispTBctpdt, cells);
+ } // if
+ const int dispTBctpdtTag = _dispTags[_material->id()];
#endif
// Loop over cells
@@ -379,7 +387,7 @@
// Restrict input fields to cell
#ifdef FASTER
- mesh->restrict(dispTBctpdt, dTag, c_index, &dispTBctpdtCell[0],
+ mesh->restrict(dispTBctpdt, dispTBctpdtTag, c_index, &dispTBctpdtCell[0],
cellVecSize);
#else
mesh->restrict(dispTBctpdt, *c_iter, &dispTBctpdtCell[0], cellVecSize);
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.hh 2007-11-17 01:30:59 UTC (rev 8300)
@@ -139,8 +139,7 @@
double _dtm1; ///< Time step for t-dt1 -> t
// Optimization
- std::map<int, int> _dTags; ///< Tags indexing dispTBctpdt field
- std::map<int, int> _rTags; ///< tags indexing residual field
+ std::map<int, int> _dispTBctpdtTags; ///< Tags indexing dispTBctpdt field.
}; // ElasticityImplicit
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2007-11-17 01:30:59 UTC (rev 8300)
@@ -15,7 +15,6 @@
#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
@@ -24,7 +23,7 @@
#include <assert.h> // USES assert()
#include <stdexcept> // USES std::runtime_error
-//#define FASTER
+#define FASTER
// ----------------------------------------------------------------------
// Constructor
@@ -80,21 +79,24 @@
// Set variables dependent on dimension of cell
const int cellDim = _quadrature->cellDim();
int tensorSize = 0;
- Elasticity::totalStrain_fn_type calcTotalStrainFn;
+ totalStrain_fn_type calcTotalStrainFn;
if (1 == cellDim) {
tensorSize = 1;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain1D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
} else if (2 == cellDim) {
tensorSize = 3;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain2D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
} else if (3 == cellDim) {
tensorSize = 6;
- calcTotalStrainFn = &pylith::feassemble::Elasticity::calcTotalStrain3D;
+ calcTotalStrainFn =
+ &pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D;
} else
assert(0);
// Get cell information
- const ALE::Obj<ALE::Mesh::label_sequence>& cells =
+ const ALE::Obj<Mesh::label_sequence>& cells =
mesh->getLabelStratum("material-id", _material->id());
assert(!cells.isNull());
const Mesh::label_sequence::iterator cellsEnd = cells->end();
@@ -120,9 +122,11 @@
} // for
#ifdef FASTER
- const int dispTTag = _dispTTags[_material->id()];
+ assert(_dispTags.find(_material->id()) != _dispTags.end());
+ const int dispTag = _dispTags[_material->id()];
#endif
-
+
+ // Loop over cells
int c_index = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
@@ -132,7 +136,7 @@
// Restrict input fields to cell
#ifdef FASTER
- mesh->restrict(disp, dispTTag, c_index, &dispCell[0], cellVecSize);
+ mesh->restrict(disp, dispTag, c_index, &dispCell[0], cellVecSize);
#else
mesh->restrict(disp, *c_iter, &dispCell[0], cellVecSize);
#endif
@@ -517,5 +521,95 @@
PetscLogFlopsNoCheck(numQuadPts*(1+numBasis*(3+numBasis*(6*26+9))));
} // _elasticityJacobian3D
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D(
+ std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis)
+{ // calcTotalStrain1D
+ assert(0 != strain);
+
+ const int dim = 1;
+ const int numQuadPts = strain->size();
+ assert(basisDeriv.size() == numQuadPts*numBasis*dim);
+ assert(disp.size() == numBasis*dim);
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ assert(1 == (*strain)[iQuad].size());
+ (*strain)[iQuad] *= 0.0;
+ for (int iBasis=0; iBasis < numBasis; ++iBasis)
+ (*strain)[iQuad][0] +=
+ basisDeriv[iQuad*numBasis+iBasis] * disp[iBasis];
+ } // for
+} // calcTotalStrain1D
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D(
+ std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis)
+{ // calcTotalStrain2D
+ assert(0 != strain);
+
+ const int dim = 2;
+ const int numQuadPts = strain->size();
+
+ assert(basisDeriv.size() == numQuadPts*numBasis*dim);
+ assert(disp.size() == numBasis*dim);
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ assert(3 == (*strain)[iQuad].size());
+ (*strain)[iQuad] *= 0.0;
+ for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
+ (*strain)[iQuad][0] += basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
+ (*strain)[iQuad][1] += basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
+ (*strain)[iQuad][2] +=
+ 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
+ } // for
+ } // for
+} // calcTotalStrain2D
+
+// ----------------------------------------------------------------------
+void
+pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D(
+ std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis)
+{ // calcTotalStrain3D
+ assert(0 != strain);
+
+ const int dim = 3;
+ const int numQuadPts = strain->size();
+
+ assert(basisDeriv.size() == numQuadPts*numBasis*dim);
+ assert(disp.size() == numBasis*dim);
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ assert(6 == (*strain)[iQuad].size());
+ (*strain)[iQuad] *= 0.0;
+ for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
+ (*strain)[iQuad][0] += basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
+ (*strain)[iQuad][1] += basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
+ (*strain)[iQuad][2] += basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+2];
+ (*strain)[iQuad][3] +=
+ 0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
+ (*strain)[iQuad][4] +=
+ 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+1] +
+ basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+2]);
+ (*strain)[iQuad][5] +=
+ 0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim ] +
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+2]);
+ } // for
+ } // for
+} // calcTotalStrain3D
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2007-11-17 01:30:59 UTC (rev 8300)
@@ -38,6 +38,15 @@
{ // IntegratorElasticity
friend class TestIntegratorElasticity; // unit testing
+// PUBLIC TYPEDEFS //////////////////////////////////////////////////////
+public :
+
+ typedef void (*totalStrain_fn_type)(std::vector<double_array>*,
+ const double_array&,
+ const double_array&,
+ const int);
+
+
// PUBLIC MEMBERS ///////////////////////////////////////////////////////
public :
@@ -114,7 +123,50 @@
*/
void _elasticityJacobian3D(const std::vector<double_array>& elasticConsts);
-// PRIVATE METHODS //////////////////////////////////////////////////////
+ /** 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 _calcTotalStrain1D(std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis);
+
+ /** Compute total strain in at quadrature points of a cell.
+ *
+ * @param strain Strain tensor at quadrature points.
+ * @param basisDeriv Derivatives of basis functions at quadrature points.
+ * @param disp Displacement at vertices of cell.
+ * @param numBasis Number of basis functions for cell.
+ */
+
+ static
+ void _calcTotalStrain2D(std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis);
+
+ /** Compute total strain in at quadrature points of a cell.
+ *
+ * @param strain Strain tensor at quadrature points.
+ * @param basisDeriv Derivatives of basis functions at quadrature points.
+ * @param disp Displacement at vertices of cell.
+ * @param numBasis Number of basis functions for cell.
+ */
+
+ static
+ void _calcTotalStrain3D(std::vector<double_array>* strain,
+ const double_array& basisDeriv,
+ const double_array& disp,
+ const int numBasis);
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
/// Not implemented.
@@ -129,6 +181,10 @@
/// Elastic material associated with integrator
materials::ElasticMaterial* _material;
+ // Optimization
+ std::map<int, int> _dispTags; ///< Tags indexing displacement field.
+ std::map<int, int> _residualTags; ///< Tags indexing residual field.
+
}; // IntegratorElasticity
#endif // pylith_feassemble_integratorelasticity_hh
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2007-11-17 01:30:59 UTC (rev 8300)
@@ -23,7 +23,6 @@
Integrator.hh \
Integrator.icc \
IntegratorElasticity.hh \
- Elasticity.hh \
GeometryPoint1D.hh \
GeometryPoint2D.hh \
GeometryPoint3D.hh \
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am 2007-11-17 01:30:59 UTC (rev 8300)
@@ -22,7 +22,6 @@
# Primary source files
testfeassemble_SOURCES = \
TestCellGeometry.cc \
- TestElasticity.cc \
TestElasticityExplicit.cc \
TestElasticityExplicit1DLinear.cc \
TestElasticityExplicit1DQuadratic.cc \
@@ -50,6 +49,7 @@
TestGeometryTet3D.cc \
TestGeometryHex3D.cc \
TestIntegrator.cc \
+ TestIntegratorElasticity.cc \
TestQuadrature.cc \
TestQuadrature0D.cc \
TestQuadrature1D.cc \
@@ -61,7 +61,6 @@
noinst_HEADERS = \
TestCellGeometry.hh \
- TestElasticity.hh \
TestElasticityExplicit.hh \
TestElasticityExplicit1DLinear.hh \
TestElasticityExplicit1DQuadratic.hh \
@@ -89,6 +88,7 @@
TestGeometryTet3D.hh \
TestGeometryHex3D.hh \
TestIntegrator.hh \
+ TestIntegratorElasticity.hh \
TestQuadrature.hh \
TestQuadrature0D.hh \
TestQuadrature1D.hh \
Deleted: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.cc 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.cc 2007-11-17 01:30:59 UTC (rev 8300)
@@ -1,181 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-#include <portinfo>
-
-#include "TestElasticity.hh" // Implementation of class methods
-
-#include "pylith/feassemble/Elasticity.hh" // USES Elasticity
-
-#include <math.h> // USES fabs()
-
-#include <stdexcept>
-// ----------------------------------------------------------------------
-CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestElasticity );
-
-// ----------------------------------------------------------------------
-// Test calcTotalStrain1D().
-void
-pylith::feassemble::TestElasticity::testCalcTotalStrain1D(void)
-{ // testCalcTotalStrain1D
- // N0 = 0.5 * (1 - x)
- // N1 = 0.5 * (1 + x)
- // dN0/dx = -0.5
- // dN1/dx = +0.5
- // Let quad pt 0 be dN/dx, let quad pt 1 be 0.5*dN/dx
- const int dim = 1;
- const int numBasis = 2;
- const int numQuadPts = 2;
- const double basisDerivVals[] = {
- -0.50, 0.50,
- -0.25, 0.25 };
- const int tensorSize = 1;
-
- // Let u(x) = 1 + 0.5 * x
- const double dispVals[] = { 0.5, 1.5 };
- const double strainE[] = { 0.5, 0.25 };
-
- std::vector<double_array> strain(numQuadPts);
- for (int i=0; i < numQuadPts; ++i)
- strain[i].resize(tensorSize);
-
- double_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
- double_array disp(dispVals, numBasis*dim);
-
- Elasticity::calcTotalStrain1D(&strain, basisDeriv, disp, numBasis);
-
- const double tolerance = 1.0e-06;
- CPPUNIT_ASSERT_EQUAL(numQuadPts, int(strain.size()));
- for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
- CPPUNIT_ASSERT_EQUAL(tensorSize, int(strain[iQuad].size()));
- for (int iStrain=0; iStrain < tensorSize; ++iStrain)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i++], strain[iQuad][iStrain],
- tolerance);
- } // for
-} // testCalcTotalStrain1D
-
-// ----------------------------------------------------------------------
-// Test calcTotalStrain2D().
-void
-pylith::feassemble::TestElasticity::testCalcTotalStrain2D(void)
-{ // testCalcTotalStrain2D
- // N0 = x
- // N1 = y
- // N2 = 1 - x - y
- // dN0/dx = +1.0, dN0/dy = 0.0
- // dN1/dx = 0.0, dN1/dy = +1.0
- // dN2/dx = -1.0, dN2/dy = -1.0
- // Let quad pt 0 be derivatives, let quad pt 1 be 2.0*derivatives
- const int dim = 2;
- const int numBasis = 3;
- const int numQuadPts = 2;
- const double basisDerivVals[] = {
- +1.0, 0.0, 0.0, +1.0, -1.0, -1.0,
- +2.0, 0.0, 0.0, +2.0, -2.0, -2.0
- };
- const int tensorSize = 3;
-
- // Let ux(x,y) = +0.4 + 0.3*x + 0.8*y
- // Ley uy(x,y) = -2.0 + 0.5*x - 0.2*y
- const double dispVals[] = {
- 0.7, -1.5,
- 1.2, -2.2,
- 0.4, -2.0
- };
- const double strainE[] = {
- 0.3, -0.2, 0.65,
- 0.6, -0.4, 1.3
- };
-
- std::vector<double_array> strain(numQuadPts);
- for (int i=0; i < numQuadPts; ++i)
- strain[i].resize(tensorSize);
-
- double_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
- double_array disp(dispVals, numBasis*dim);
-
- Elasticity::calcTotalStrain2D(&strain, basisDeriv, disp, numBasis);
-
- const double tolerance = 1.0e-06;
- CPPUNIT_ASSERT_EQUAL(numQuadPts, int(strain.size()));
- for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
- CPPUNIT_ASSERT_EQUAL(tensorSize, int(strain[iQuad].size()));
- for (int iStrain=0; iStrain < tensorSize; ++iStrain)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i++], strain[iQuad][iStrain],
- tolerance);
- } // for
-} // testCalcTotalStrain2D
-
-// ----------------------------------------------------------------------
-// Test calcTotalStrain3D().
-void
-pylith::feassemble::TestElasticity::testCalcTotalStrain3D(void)
-{ // testCalcTotalStrain3D
- // N0 = x
- // N1 = y
- // N2 = z
- // N2 = 1 - x - y - z
- // dN0/dx = +1.0, dN0/dy = 0.0, dN0/dz = 0.0
- // dN1/dx = 0.0, dN1/dy = +1.0, dN0/dz = 0.0
- // dN2/dx = 0.0, dN2/dy = 0.0, dN2/dz = +1.0
- // dN3/dx = -1.0, dN3/dy = -1.0, dN3/dz = -1.0
- // Let quad pt 0 be derivatives, let quad pt 1 be 3.0*derivatives
- const int dim = 3;
- const int numBasis = 4;
- const int numQuadPts = 2;
- const double basisDerivVals[] = {
- +1.0, 0.0, 0.0, // Quad pt 0
- 0.0, +1.0, 0.0,
- 0.0, 0.0, +1.0,
- -1.0, -1.0, -1.0,
- +3.0, 0.0, 0.0, // Quad pt 1
- 0.0, +3.0, 0.0,
- 0.0, 0.0, +3.0,
- -3.0, -3.0, -3.0
- };
- const int tensorSize = 6;
-
- // Let ux(x,y,z) = +0.4 + 0.3*x + 0.8*y + 0.4*z
- // Ley uy(x,y,z) = -2.0 + 0.5*x - 0.2*y + 1.2*z
- // Ley uz(x,y,z) = -1.0 + 0.2*x - 0.7*y - 0.3*z
- const double dispVals[] = {
- 0.7, -1.5, -0.8,
- 1.2, -2.2, -1.7,
- 0.8, -0.8, -1.3,
- 0.4, -2.0, -1.0
- };
- const double strainE[] = {
- 0.3, -0.2, -0.3, 0.65, 0.25, 0.3,
- 0.9, -0.6, -0.9, 1.95, 0.75, 0.9
- };
-
- std::vector<double_array> strain(numQuadPts);
- for (int i=0; i < numQuadPts; ++i)
- strain[i].resize(tensorSize);
-
- double_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
- double_array disp(dispVals, numBasis*dim);
-
- Elasticity::calcTotalStrain3D(&strain, basisDeriv, disp, numBasis);
-
- const double tolerance = 1.0e-06;
- CPPUNIT_ASSERT_EQUAL(numQuadPts, int(strain.size()));
- for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
- CPPUNIT_ASSERT_EQUAL(tensorSize, int(strain[iQuad].size()));
- for (int iStrain=0; iStrain < tensorSize; ++iStrain)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i++], strain[iQuad][iStrain],
- tolerance);
- } // for
-} // testCalcTotalStrain3D
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.hh 2007-11-16 23:40:38 UTC (rev 8299)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.hh 2007-11-17 01:30:59 UTC (rev 8300)
@@ -1,63 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-/**
- * @file unittests/libtests/feassemble/TestElasticity.hh
- *
- * @brief C++ TestElasticity object
- *
- * C++ unit testing for Elasticity.
- */
-
-#if !defined(pylith_feassemble_testelasticity_hh)
-#define pylith_feassemble_testelasticity_hh
-
-#include <cppunit/extensions/HelperMacros.h>
-
-/// Namespace for pylith package
-namespace pylith {
- namespace feassemble {
- class TestElasticity;
- } // feassemble
-} // pylith
-
-/// C++ unit testing for Elasticity
-class pylith::feassemble::TestElasticity : public CppUnit::TestFixture
-{ // class TestElasticity
-
- // CPPUNIT TEST SUITE /////////////////////////////////////////////////
- CPPUNIT_TEST_SUITE( TestElasticity );
-
- CPPUNIT_TEST( testCalcTotalStrain1D );
- CPPUNIT_TEST( testCalcTotalStrain2D );
- CPPUNIT_TEST( testCalcTotalStrain3D );
-
- CPPUNIT_TEST_SUITE_END();
-
- // PUBLIC METHODS /////////////////////////////////////////////////////
-public :
-
- /// Test calcTotalStrain1D().
- void testCalcTotalStrain1D(void);
-
- /// Test calcTotalStrain2D().
- void testCalcTotalStrain2D(void);
-
- /// Test calcTotalStrain3D().
- void testCalcTotalStrain3D(void);
-
-}; // class TestElasticity
-
-#endif // pylith_feassemble_testelasticity_hh
-
-
-// End of file
Copied: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.cc (from rev 8298, short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.cc)
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.cc 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.cc 2007-11-17 01:30:59 UTC (rev 8300)
@@ -0,0 +1,184 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestIntegratorElasticity.hh" // Implementation of class methods
+
+#include "pylith/feassemble/IntegratorElasticity.hh" // USES IntegratorElasticity
+
+#include <math.h> // USES fabs()
+
+#include <stdexcept>
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestIntegratorElasticity );
+
+// ----------------------------------------------------------------------
+// Test calcTotalStrain1D().
+void
+pylith::feassemble::TestIntegratorElasticity::testCalcTotalStrain1D(void)
+{ // testCalcTotalStrain1D
+ // N0 = 0.5 * (1 - x)
+ // N1 = 0.5 * (1 + x)
+ // dN0/dx = -0.5
+ // dN1/dx = +0.5
+ // Let quad pt 0 be dN/dx, let quad pt 1 be 0.5*dN/dx
+ const int dim = 1;
+ const int numBasis = 2;
+ const int numQuadPts = 2;
+ const double basisDerivVals[] = {
+ -0.50, 0.50,
+ -0.25, 0.25 };
+ const int tensorSize = 1;
+
+ // Let u(x) = 1 + 0.5 * x
+ const double dispVals[] = { 0.5, 1.5 };
+ const double strainE[] = { 0.5, 0.25 };
+
+ std::vector<double_array> strain(numQuadPts);
+ for (int i=0; i < numQuadPts; ++i)
+ strain[i].resize(tensorSize);
+
+ double_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
+ double_array disp(dispVals, numBasis*dim);
+
+ IntegratorElasticity::_calcTotalStrain1D(&strain,
+ basisDeriv, disp, numBasis);
+
+ const double tolerance = 1.0e-06;
+ CPPUNIT_ASSERT_EQUAL(numQuadPts, int(strain.size()));
+ for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
+ CPPUNIT_ASSERT_EQUAL(tensorSize, int(strain[iQuad].size()));
+ for (int iStrain=0; iStrain < tensorSize; ++iStrain)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i++], strain[iQuad][iStrain],
+ tolerance);
+ } // for
+} // testCalcTotalStrain1D
+
+// ----------------------------------------------------------------------
+// Test calcTotalStrain2D().
+void
+pylith::feassemble::TestIntegratorElasticity::testCalcTotalStrain2D(void)
+{ // testCalcTotalStrain2D
+ // N0 = x
+ // N1 = y
+ // N2 = 1 - x - y
+ // dN0/dx = +1.0, dN0/dy = 0.0
+ // dN1/dx = 0.0, dN1/dy = +1.0
+ // dN2/dx = -1.0, dN2/dy = -1.0
+ // Let quad pt 0 be derivatives, let quad pt 1 be 2.0*derivatives
+ const int dim = 2;
+ const int numBasis = 3;
+ const int numQuadPts = 2;
+ const double basisDerivVals[] = {
+ +1.0, 0.0, 0.0, +1.0, -1.0, -1.0,
+ +2.0, 0.0, 0.0, +2.0, -2.0, -2.0
+ };
+ const int tensorSize = 3;
+
+ // Let ux(x,y) = +0.4 + 0.3*x + 0.8*y
+ // Ley uy(x,y) = -2.0 + 0.5*x - 0.2*y
+ const double dispVals[] = {
+ 0.7, -1.5,
+ 1.2, -2.2,
+ 0.4, -2.0
+ };
+ const double strainE[] = {
+ 0.3, -0.2, 0.65,
+ 0.6, -0.4, 1.3
+ };
+
+ std::vector<double_array> strain(numQuadPts);
+ for (int i=0; i < numQuadPts; ++i)
+ strain[i].resize(tensorSize);
+
+ double_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
+ double_array disp(dispVals, numBasis*dim);
+
+ IntegratorElasticity::_calcTotalStrain2D(&strain,
+ basisDeriv, disp, numBasis);
+
+ const double tolerance = 1.0e-06;
+ CPPUNIT_ASSERT_EQUAL(numQuadPts, int(strain.size()));
+ for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
+ CPPUNIT_ASSERT_EQUAL(tensorSize, int(strain[iQuad].size()));
+ for (int iStrain=0; iStrain < tensorSize; ++iStrain)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i++], strain[iQuad][iStrain],
+ tolerance);
+ } // for
+} // testCalcTotalStrain2D
+
+// ----------------------------------------------------------------------
+// Test calcTotalStrain3D().
+void
+pylith::feassemble::TestIntegratorElasticity::testCalcTotalStrain3D(void)
+{ // testCalcTotalStrain3D
+ // N0 = x
+ // N1 = y
+ // N2 = z
+ // N2 = 1 - x - y - z
+ // dN0/dx = +1.0, dN0/dy = 0.0, dN0/dz = 0.0
+ // dN1/dx = 0.0, dN1/dy = +1.0, dN0/dz = 0.0
+ // dN2/dx = 0.0, dN2/dy = 0.0, dN2/dz = +1.0
+ // dN3/dx = -1.0, dN3/dy = -1.0, dN3/dz = -1.0
+ // Let quad pt 0 be derivatives, let quad pt 1 be 3.0*derivatives
+ const int dim = 3;
+ const int numBasis = 4;
+ const int numQuadPts = 2;
+ const double basisDerivVals[] = {
+ +1.0, 0.0, 0.0, // Quad pt 0
+ 0.0, +1.0, 0.0,
+ 0.0, 0.0, +1.0,
+ -1.0, -1.0, -1.0,
+ +3.0, 0.0, 0.0, // Quad pt 1
+ 0.0, +3.0, 0.0,
+ 0.0, 0.0, +3.0,
+ -3.0, -3.0, -3.0
+ };
+ const int tensorSize = 6;
+
+ // Let ux(x,y,z) = +0.4 + 0.3*x + 0.8*y + 0.4*z
+ // Ley uy(x,y,z) = -2.0 + 0.5*x - 0.2*y + 1.2*z
+ // Ley uz(x,y,z) = -1.0 + 0.2*x - 0.7*y - 0.3*z
+ const double dispVals[] = {
+ 0.7, -1.5, -0.8,
+ 1.2, -2.2, -1.7,
+ 0.8, -0.8, -1.3,
+ 0.4, -2.0, -1.0
+ };
+ const double strainE[] = {
+ 0.3, -0.2, -0.3, 0.65, 0.25, 0.3,
+ 0.9, -0.6, -0.9, 1.95, 0.75, 0.9
+ };
+
+ std::vector<double_array> strain(numQuadPts);
+ for (int i=0; i < numQuadPts; ++i)
+ strain[i].resize(tensorSize);
+
+ double_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
+ double_array disp(dispVals, numBasis*dim);
+
+ IntegratorElasticity::_calcTotalStrain3D(&strain,
+ basisDeriv, disp, numBasis);
+
+ const double tolerance = 1.0e-06;
+ CPPUNIT_ASSERT_EQUAL(numQuadPts, int(strain.size()));
+ for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
+ CPPUNIT_ASSERT_EQUAL(tensorSize, int(strain[iQuad].size()));
+ for (int iStrain=0; iStrain < tensorSize; ++iStrain)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i++], strain[iQuad][iStrain],
+ tolerance);
+ } // for
+} // testCalcTotalStrain3D
+
+
+// End of file
Copied: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.hh (from rev 8298, short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.hh)
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticity.hh 2007-11-16 21:58:35 UTC (rev 8298)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.hh 2007-11-17 01:30:59 UTC (rev 8300)
@@ -0,0 +1,63 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/feassemble/TestIntegratorElasticity.hh
+ *
+ * @brief C++ TestIntegratorElasticity object
+ *
+ * C++ unit testing for IntegratorElasticity.
+ */
+
+#if !defined(pylith_feassemble_testintegratorelasticity_hh)
+#define pylith_feassemble_testintegratorelasticity_hh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace feassemble {
+ class TestIntegratorElasticity;
+ } // feassemble
+} // pylith
+
+/// C++ unit testing for Elasticity
+class pylith::feassemble::TestIntegratorElasticity : public CppUnit::TestFixture
+{ // class TestIntegratorElasticity
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestIntegratorElasticity );
+
+ CPPUNIT_TEST( testCalcTotalStrain1D );
+ CPPUNIT_TEST( testCalcTotalStrain2D );
+ CPPUNIT_TEST( testCalcTotalStrain3D );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Test calcTotalStrain1D().
+ void testCalcTotalStrain1D(void);
+
+ /// Test calcTotalStrain2D().
+ void testCalcTotalStrain2D(void);
+
+ /// Test calcTotalStrain3D().
+ void testCalcTotalStrain3D(void);
+
+}; // class TestIntegratorElasticity
+
+#endif // pylith_feassemble_testintegratorelasticity_hh
+
+
+// End of file
More information about the cig-commits
mailing list