[cig-commits] r6533 - in short/3D/PyLith/trunk: libsrc/feassemble
libsrc/materials libsrc/utils unittests/libtests/materials
brad at geodynamics.org
brad at geodynamics.org
Mon Apr 9 18:28:52 PDT 2007
Author: brad
Date: 2007-04-09 18:28:52 -0700 (Mon, 09 Apr 2007)
New Revision: 6533
Added:
short/3D/PyLith/trunk/libsrc/utils/stlfwd.hh
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh
short/3D/PyLith/trunk/libsrc/materials/Material.cc
short/3D/PyLith/trunk/libsrc/materials/Material.hh
short/3D/PyLith/trunk/libsrc/utils/Makefile.am
short/3D/PyLith/trunk/libsrc/utils/petscfwd.h
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
Log:
Switch to using C++ STL containers in materials. This reduces clutter. Also improved memory management.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc 2007-04-10 01:28:52 UTC (rev 6533)
@@ -20,6 +20,7 @@
#include "petscmat.h" // USES PetscMat
#include "spatialdata/spatialdb/SpatialDB.hh"
+#include <valarray> // USES std::valarray (double_array)
#include <assert.h> // USES assert()
#include <stdexcept> // USES std::runtime_error
@@ -119,10 +120,11 @@
// Compute action for cell
// Compute action for inertial terms
- const double* density = _material->calcDensity(*cellIter, numQuadPts);
+ const std::vector<double_array>& density =
+ _material->calcDensity(*cellIter);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt =
- quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
+ quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad][0] / dt2;
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
const double valI = wt*basis[iQ+iBasis];
@@ -154,62 +156,57 @@
// Compute action for elastic terms
if (1 == cellDim) {
// Compute total strains
- const int stressSize = _material->stressSize();
- assert(1 == stressSize);
- const int strainSize = stressSize * numQuadPts;
- double* totalStrain = (strainSize > 0) ? new double[strainSize] : 0;
- memset(totalStrain, 0, strainSize*sizeof(double));
+ const int tensorSize = 1;
+ std::vector<double_array> totalStrain(numQuadPts);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ totalStrain[iQuad].resize(tensorSize);
+ totalStrain[iQuad] *= 0.0;
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis)
- totalStrain[iQuad] += basisDeriv[iQ+iBasis] * dispTCell[iBasis];
+ totalStrain[iQuad][0] += basisDeriv[iQ+iBasis] * dispTCell[iBasis];
} // for
- const double* stress =
- _material->calcStress(*cellIter, totalStrain, numQuadPts,
- spaceDim);
+ const std::vector<double_array>& stress =
+ _material->calcStress(*cellIter, totalStrain);
// Compute elastic action
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const int iStress = iQuad*stressSize;
- const double s11 = stress[iStress ];
+ const double s11 = stress[iQuad][0];
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
const double N1 = wt*basisDeriv[iQ+iBasis*cellDim ];
_cellVector[iBlock ] -= N1*s11;
} // for
} // for
- delete[] totalStrain; totalStrain = 0;
PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*5));
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
+
} else if (2 == cellDim) {
// Compute total strains
- const int stressSize = _material->stressSize();
- assert(3 == stressSize);
- const int strainSize = stressSize * numQuadPts;
- double* totalStrain = (strainSize > 0) ? new double[strainSize] : 0;
- memset(totalStrain, 0, strainSize*sizeof(double));
+ const int tensorSize = 3;
+ std::vector<double_array> totalStrain(numQuadPts);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ totalStrain[iQuad].resize(tensorSize);
+ totalStrain[iQuad] *= 0.0;
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- totalStrain[iQuad ] +=
+ totalStrain[iQuad][0] +=
basisDeriv[iQ+iBasis ] * dispTCell[iBasis ];
- totalStrain[iQuad+1] +=
+ totalStrain[iQuad][1] +=
basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+1];
- totalStrain[iQuad+2] +=
+ totalStrain[iQuad][2] +=
0.5 * (basisDeriv[iQ+iBasis+1] * dispTCell[iBasis ] +
basisDeriv[iQ+iBasis ] * dispTCell[iBasis+1]);
} // for
} // for
- const double* stress =
- _material->calcStress(*cellIter, totalStrain, numQuadPts,
- spaceDim);
+ const std::vector<double_array>& stress =
+ _material->calcStress(*cellIter, totalStrain);
+
// Compute elastic action
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const int iStress = iQuad*stressSize;
- const double s11 = stress[iStress ];
- const double s22 = stress[iStress+1];
- const double s12 = stress[iStress+2];
+ 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; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
const double N1 = wt*basisDeriv[iQ+iBasis*cellDim ];
@@ -221,44 +218,44 @@
err = PetscLogFlops(numQuadPts*(1+numBasis*(8+2+9)));
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
+
} else if (3 == cellDim) {
// Compute total strains
- const int stressSize = _material->stressSize();
- assert(6 == stressSize);
- const int strainSize = stressSize*numQuadPts;
- double* totalStrain = (strainSize > 0) ? new double[strainSize] : 0;
- memset(totalStrain, 0, strainSize*sizeof(double));
+ const int tensorSize = 6;
+ std::vector<double_array> totalStrain(numQuadPts);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ totalStrain[iQuad].resize(tensorSize);
+ totalStrain[iQuad] *= 0.0;
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
- totalStrain[iQuad ] +=
+ totalStrain[iQuad][0] +=
basisDeriv[iQ+iBasis ] * dispTCell[iBasis ];
- totalStrain[iQuad+1] +=
+ totalStrain[iQuad][1] +=
basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+1];
- totalStrain[iQuad+2] +=
+ totalStrain[iQuad][2] +=
basisDeriv[iQ+iBasis+2] * dispTCell[iBasis+2];
- totalStrain[iQuad+3] +=
+ totalStrain[iQuad][3] +=
0.5 * (basisDeriv[iQ+iBasis+1] * dispTCell[iBasis ] +
basisDeriv[iQ+iBasis ] * dispTCell[iBasis+1]);
- totalStrain[iQuad+4] +=
+ totalStrain[iQuad][4] +=
0.5 * (basisDeriv[iQ+iBasis+2] * dispTCell[iBasis+1] +
basisDeriv[iQ+iBasis+1] * dispTCell[iBasis+2]);
- totalStrain[iQuad+5] +=
+ totalStrain[iQuad][5] +=
0.5 * (basisDeriv[iQ+iBasis+2] * dispTCell[iBasis ] +
basisDeriv[iQ+iBasis ] * dispTCell[iBasis+2]);
} // for
} // for
- const double* stress =
- _material->calcStress(*cellIter, totalStrain, numQuadPts, spaceDim);
+ const std::vector<double_array>& stress =
+ _material->calcStress(*cellIter, totalStrain);
+
// Compute elastic action
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const int iStress = iQuad*stressSize;
- const double s11 = stress[iStress ];
- const double s22 = stress[iStress+1];
- const double s33 = stress[iStress+2];
- const double s12 = stress[iStress+3];
- const double s23 = stress[iStress+4];
- const double s13 = stress[iStress+5];
+ 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; iBasis < numBasis; ++iBasis) {
const int iBlock = iBasis * spaceDim;
@@ -330,14 +327,15 @@
const double* jacobianDet = _quadrature->jacobianDet();
// Get material physical properties at quadrature points for this cell
- const double* density = _material->calcDensity(*cellIter, numQuadPts);
+ const std::vector<double_array>& density =
+ _material->calcDensity(*cellIter);
// Compute Jacobian for cell
// Compute Jacobian for inertial terms
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt =
- quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
+ quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad][0] / dt2;
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const double valI = wt*basis[iQ+iBasis];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2007-04-10 01:28:52 UTC (rev 6533)
@@ -14,6 +14,7 @@
#include "ElasticIsotropic3D.hh" // implementation of object methods
+#include <valarray> // USES std::valarray (double_array)
#include <assert.h> // USES assert()
// ----------------------------------------------------------------------
@@ -29,8 +30,8 @@
*/
class pylith::materials::_ElasticIsotropic3D {
public:
- /// Number of entries in stress tensor.
- static const int stressSize;
+ /// Number of entries in stress/strain tensors.
+ static const int tensorSize;
/// Number of entries in derivative of elasticity matrix.
static const int numElasticConsts;
@@ -54,7 +55,7 @@
static const int pidLambda;
}; // _ElasticIsotropic3D
-const int pylith::materials::_ElasticIsotropic3D::stressSize = 6;
+const int pylith::materials::_ElasticIsotropic3D::tensorSize = 6;
const int pylith::materials::_ElasticIsotropic3D::numElasticConsts = 21;
const int pylith::materials::_ElasticIsotropic3D::numDBValues = 3;
const char* pylith::materials::_ElasticIsotropic3D::namesDBValues[] =
@@ -91,22 +92,6 @@
} // copy constructor
// ----------------------------------------------------------------------
-// Get number of entries in stress tensor.
-int
-pylith::materials::ElasticIsotropic3D::stressSize(void) const
-{ // stressSize
- return _ElasticIsotropic3D::stressSize;
-} // stressSize
-
-// ----------------------------------------------------------------------
-// Get number of entries in elasticity matrix for material.
-int
-pylith::materials::ElasticIsotropic3D::numElasticConsts(void) const
-{ // numElasticConsts
- return _ElasticIsotropic3D::numElasticConsts;
-} // numElasticConsts
-
-// ----------------------------------------------------------------------
// Get names of values expected to be in database of parameters for
const char**
pylith::materials::ElasticIsotropic3D::_dbValues(void) const
@@ -141,15 +126,14 @@
// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
-pylith::materials::ElasticIsotropic3D::_dbToParameters(double* paramVals,
- const int numParams,
- const double* dbValues,
- const int numValues) const
-{ // computeParameters
+pylith::materials::ElasticIsotropic3D::_dbToParameters(double_array* paramVals,
+ const double_array& dbValues) const
+{ // _dbToParameters
assert(0 != paramVals);
+ const int numParams = paramVals->size();
assert(_ElasticIsotropic3D::numParameters == numParams);
- assert(0 != dbValues);
- assert(_ElasticIsotropic3D::numDBValues == numValues);
+ const int numDBValues = dbValues.size();
+ assert(_ElasticIsotropic3D::numDBValues == numDBValues);
const double density = dbValues[_ElasticIsotropic3D::didDensity];
const double vs = dbValues[_ElasticIsotropic3D::didVs];
@@ -158,41 +142,51 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- paramVals[_ElasticIsotropic3D::pidDensity] = density;
- paramVals[_ElasticIsotropic3D::pidMu] = mu;
- paramVals[_ElasticIsotropic3D::pidLambda] = lambda;
-} // computeParameters
+ (*paramVals)[_ElasticIsotropic3D::pidDensity] = density;
+ (*paramVals)[_ElasticIsotropic3D::pidMu] = mu;
+ (*paramVals)[_ElasticIsotropic3D::pidLambda] = lambda;
+} // _dbToParameters
// ----------------------------------------------------------------------
+// Get number of entries in stress tensor.
+int
+pylith::materials::ElasticIsotropic3D::_tensorSize(void) const
+{ // _tensorSize
+ return _ElasticIsotropic3D::tensorSize;
+} // _tensorSize
+
+// ----------------------------------------------------------------------
+// Get number of entries in elasticity matrix for material.
+int
+pylith::materials::ElasticIsotropic3D::_numElasticConsts(void) const
+{ // numElasticConsts
+ return _ElasticIsotropic3D::numElasticConsts;
+} // numElasticConsts
+
+// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::ElasticIsotropic3D::_calcDensity(double* const density,
- const int size,
- const double* parameters,
- const int numParameters)
+pylith::materials::ElasticIsotropic3D::_calcDensity(double_array* const density,
+ const double_array& parameters)
{ // _calcDensity
assert(0 != density);
- assert(1 == size);
- assert(0 != parameters);
+ assert(1 == density->size());
+ assert(_ElasticIsotropic3D::numParameters == parameters.size());
- *density = parameters[_ElasticIsotropic3D::pidDensity];
+ (*density)[0] = parameters[_ElasticIsotropic3D::pidDensity];
} // _calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor at location from parameters.
void
-pylith::materials::ElasticIsotropic3D::_calcStress(double* const stress,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim)
+pylith::materials::ElasticIsotropic3D::_calcStress(double_array* const stress,
+ const double_array& parameters,
+ const double_array& totalStrain)
{ // _calcStress
assert(0 != stress);
- assert(_ElasticIsotropic3D::stressSize == size);
- assert(0 != parameters);
- assert(_ElasticIsotropic3D::numParameters == numParameters);
- assert(spaceDim == _dimension);
+ assert(_ElasticIsotropic3D::tensorSize == stress->size());
+ assert(_ElasticIsotropic3D::numParameters == parameters.size());
+ assert(_ElasticIsotropic3D::tensorSize == totalStrain.size());
const double density = parameters[_ElasticIsotropic3D::pidDensity];
const double mu = parameters[_ElasticIsotropic3D::pidMu];
@@ -209,29 +203,26 @@
const double s123 = lambda * (e11 + e22 + e33);
- stress[0] = s123 + 2.0*mu*e11;
- stress[1] = s123 + 2.0*mu*e22;
- stress[2] = s123 + 2.0*mu*e33;
- stress[3] = 2.0 * mu * e12;
- stress[4] = 2.0 * mu * e23;
- stress[5] = 2.0 * mu * e13;
+ (*stress)[0] = s123 + 2.0*mu*e11;
+ (*stress)[1] = s123 + 2.0*mu*e22;
+ (*stress)[2] = s123 + 2.0*mu*e33;
+ (*stress)[3] = 2.0 * mu * e12;
+ (*stress)[4] = 2.0 * mu * e23;
+ (*stress)[5] = 2.0 * mu * e13;
} // _calcStress
// ----------------------------------------------------------------------
// Compute derivative of elasticity matrix at location from parameters.
void
-pylith::materials::ElasticIsotropic3D::_calcElasticConsts(double* const elasticConsts,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim)
+pylith::materials::ElasticIsotropic3D::_calcElasticConsts(
+ double_array* const elasticConsts,
+ const double_array& parameters,
+ const double_array& totalStrain)
{ // _calcElasticConsts
assert(0 != elasticConsts);
- assert(_ElasticIsotropic3D::numElasticConsts == size);
- assert(0 != parameters);
- assert(_ElasticIsotropic3D::numParameters == numParameters);
- assert(spaceDim == _dimension);
+ assert(_ElasticIsotropic3D::numElasticConsts == elasticConsts->size());
+ assert(_ElasticIsotropic3D::numParameters == parameters.size());
+ assert(_ElasticIsotropic3D::tensorSize == totalStrain.size());
const double density = parameters[_ElasticIsotropic3D::pidDensity];
const double mu = parameters[_ElasticIsotropic3D::pidMu];
@@ -239,27 +230,27 @@
const double lambda2mu = lambda + 2.0 * mu;
- elasticConsts[ 0] = lambda2mu; // C1111
- elasticConsts[ 1] = lambda; // C1122
- elasticConsts[ 2] = lambda; // C1133
- elasticConsts[ 3] = 0; // C1112
- elasticConsts[ 4] = 0; // C1123
- elasticConsts[ 5] = 0; // C1113
- elasticConsts[ 6] = lambda2mu; // C2222
- elasticConsts[ 7] = lambda; // C2233
- elasticConsts[ 8] = 0; // C2212
- elasticConsts[ 9] = 0; // C2223
- elasticConsts[10] = 0; // C2213
- elasticConsts[11] = lambda2mu; // C3333
- elasticConsts[12] = 0; // C3312
- elasticConsts[13] = 0; // C3323
- elasticConsts[14] = 0; // C3313
- elasticConsts[15] = 2.0 * mu; // C1212
- elasticConsts[16] = 0; // C1223
- elasticConsts[17] = 0; // C1213
- elasticConsts[18] = 2.0 * mu; // C2323
- elasticConsts[19] = 0; // C2313
- elasticConsts[20] = 2.0 * mu; // C1313
+ (*elasticConsts)[ 0] = lambda2mu; // C1111
+ (*elasticConsts)[ 1] = lambda; // C1122
+ (*elasticConsts)[ 2] = lambda; // C1133
+ (*elasticConsts)[ 3] = 0; // C1112
+ (*elasticConsts)[ 4] = 0; // C1123
+ (*elasticConsts)[ 5] = 0; // C1113
+ (*elasticConsts)[ 6] = lambda2mu; // C2222
+ (*elasticConsts)[ 7] = lambda; // C2233
+ (*elasticConsts)[ 8] = 0; // C2212
+ (*elasticConsts)[ 9] = 0; // C2223
+ (*elasticConsts)[10] = 0; // C2213
+ (*elasticConsts)[11] = lambda2mu; // C3333
+ (*elasticConsts)[12] = 0; // C3312
+ (*elasticConsts)[13] = 0; // C3323
+ (*elasticConsts)[14] = 0; // C3313
+ (*elasticConsts)[15] = 2.0 * mu; // C1212
+ (*elasticConsts)[16] = 0; // C1223
+ (*elasticConsts)[17] = 0; // C1213
+ (*elasticConsts)[18] = 2.0 * mu; // C2323
+ (*elasticConsts)[19] = 0; // C2313
+ (*elasticConsts)[20] = 2.0 * mu; // C1313
} // _calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh 2007-04-10 01:28:52 UTC (rev 6533)
@@ -55,26 +55,6 @@
*/
ElasticMaterial* clone(void) const;
- /** Get number of entries in stress tensor.
- *
- * 1-D = 1
- * 2-D = 3
- * 3-D = 6
- *
- * @returns Number of entries in stress tensor.
- */
- int stressSize(void) const;
-
- /** Get number of entries in derivative of elasticity matrix.
- *
- * 1-D = 1
- * 2-D = 6
- * 3-D = 21
- *
- * @returns Number of entries in derivative of elasticity matrix.
- */
- int numElasticConsts(void) const;
-
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -116,60 +96,59 @@
* parameterNames().
*
* @param paramVals Array of parameters
- * @param numParams Number of parameters
* @param dbValues Array of database values
- * @param numValues Number of database values
*/
- void _dbToParameters(double* paramVals,
- const int numParams,
- const double* dbValues,
- const int numValues) const;
+ void _dbToParameters(double_array* paramVals,
+ const double_array& dbValues) const;
+ /** Get number of entries in stress/strain tensors.
+ *
+ * 1-D = 1
+ * 2-D = 3
+ * 3-D = 6
+ *
+ * @returns Number of entries in stress/strain tensors.
+ */
+ int _tensorSize(void) const;
+
+ /** Get number of entries in derivative of elasticity matrix.
+ *
+ * 1-D = 1
+ * 2-D = 6
+ * 3-D = 21
+ *
+ * @returns Number of entries in derivative of elasticity matrix.
+ */
+ int _numElasticConsts(void) const;
+
/** Compute density from parameters.
*
* @param density Array for density
- * @param size Size of array for density
* @param parameters Parameters at location
- * @param numParameters Number of parameters
*/
- void _calcDensity(double* const density,
- const int size,
- const double* parameters,
- const int numParameters);
+ void _calcDensity(double_array* const density,
+ const double_array& parameters);
/** Compute stress tensor from parameters.
*
* @param stress Array for stress tensor
- * @param size Size of array for stress tensor
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
- void _calcStress(double* const stress,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim);
+ void _calcStress(double_array* const stress,
+ const double_array& parameters,
+ const double_array& totalStrain);
/** Compute derivatives of elasticity matrix from parameters.
*
* @param elasticConsts Array for elastic constants
- * @param size Size of array
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
- void _calcElasticConsts(double* const elasticConsts,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim);
+ void _calcElasticConsts(double_array* const elasticConsts,
+ const double_array& parameters,
+ const double_array& totalStrain);
-
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2007-04-10 01:28:52 UTC (rev 6533)
@@ -18,14 +18,13 @@
#include <petscmesh.h> // USES Mesh
+#include <valarray> // USES std::valarray (double_array)
#include <assert.h> // USES assert()
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::ElasticMaterial::ElasticMaterial(void) :
- _density(0),
- _stress(0),
- _elasticConsts(0)
+ _numQuadPts(0)
{ // constructor
} // constructor
@@ -33,98 +32,69 @@
// Destructor.
pylith::materials::ElasticMaterial::~ElasticMaterial(void)
{ // destructor
- delete[] _density; _density = 0;
- delete[] _stress; _stress = 0;
- delete[] _elasticConsts; _elasticConsts = 0;
} // destructor
// ----------------------------------------------------------------------
// Copy constructor.
pylith::materials::ElasticMaterial::ElasticMaterial(const ElasticMaterial& m) :
Material(m),
- _density(0),
- _stress(0),
- _elasticConsts(0)
+ _numQuadPts(0)
{ // copy constructor
} // copy constructor
// ----------------------------------------------------------------------
// Compute density for cell at quadrature points.
-const double*
-pylith::materials::ElasticMaterial::calcDensity(
- const Mesh::point_type& cell,
- const int numQuadPts)
+const std::vector<pylith::double_array>&
+pylith::materials::ElasticMaterial::calcDensity(const Mesh::point_type& cell)
{ // calcDensity
- assert(0 != _density);
- int size = numQuadPts;
- memset(_density, 0, size*sizeof(double));
+ _getParameters(cell);
- double* paramsCell = 0;
- _getParameters(¶msCell, cell, numQuadPts);
- const int numParams = _numParameters();
- for (int iQuad=0, indexP=0; iQuad < numQuadPts; ++iQuad, indexP+=numParams)
- _calcDensity(&_density[iQuad], 1, ¶msCell[indexP], numParams);
- delete[] paramsCell; paramsCell = 0;
+ const int numQuadPts = _numQuadPts;
+ assert(_paramsCell.size() == numQuadPts);
+ assert(_density.size() == numQuadPts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ _calcDensity(&_density[iQuad], _paramsCell[iQuad]);
+
return _density;
} // calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor for cell at quadrature points.
-const double*
+const std::vector<pylith::double_array>&
pylith::materials::ElasticMaterial::calcStress(const Mesh::point_type& cell,
- const double* totalStrain,
- const int numQuadPts,
- const int spaceDim)
+ const std::vector<double_array>& totalStrain)
{ // calcStress
- assert(0 != totalStrain);
- assert(0 != _stress);
+ _getParameters(cell);
+
+ const int numQuadPts = _numQuadPts;
+ assert(_paramsCell.size() == numQuadPts);
+ assert(totalStrain.size() == numQuadPts);
+ assert(_stress.size() == numQuadPts);
- int size = numQuadPts * stressSize();
- memset(_stress, 0, size*sizeof(double));
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ _calcStress(&_stress[iQuad], _paramsCell[iQuad], totalStrain[iQuad]);
- double* paramsCell = 0;
- _getParameters(¶msCell, cell, numQuadPts);
- const int numParams = _numParameters();
- const int stressOffset = stressSize();
- for (int iQuad=0, iStress=0, indexP=0;
- iQuad < numQuadPts;
- ++iQuad, iStress+=stressOffset, indexP+=numParams)
- _calcStress(&_stress[iStress], stressOffset,
- ¶msCell[indexP], numParams,
- &totalStrain[iStress], spaceDim);
- delete[] paramsCell; paramsCell = 0;
-
return _stress;
} // calcStress
// ----------------------------------------------------------------------
// Compute derivative of elasticity matrix for cell at quadrature points.
-const double*
+const std::vector<pylith::double_array>&
pylith::materials::ElasticMaterial::calcDerivElastic(
- const Mesh::point_type& cell,
- const double* totalStrain,
- const int numQuadPts,
- const int spaceDim)
+ const Mesh::point_type& cell,
+ const std::vector<double_array>& totalStrain)
{ // calcDerivElastic
- assert(0 != totalStrain);
- assert(0 != _elasticConsts);
+ _getParameters(cell);
- int size = numQuadPts * numElasticConsts();
- memset(_elasticConsts, 0, size*sizeof(double));
+ const int numQuadPts = _numQuadPts;
+ assert(_paramsCell.size() == numQuadPts);
+ assert(totalStrain.size() == numQuadPts);
+ assert(_elasticConsts.size() == numQuadPts);
- double* paramsCell = 0;
- _getParameters(¶msCell, cell, numQuadPts);
- const int numParams = _numParameters();
- const int stressOffset = stressSize();
- const int constOffset = numElasticConsts();
- for (int iQuad=0, iStress=0, iConst=0, indexP=0;
- iQuad < numQuadPts;
- ++iQuad, iStress+=stressOffset, iConst += constOffset, indexP+=numParams)
- _calcElasticConsts(&_elasticConsts[iConst], constOffset,
- ¶msCell[indexP], numParams,
- &totalStrain[iStress], spaceDim);
- delete[] paramsCell; paramsCell = 0;
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+ _calcElasticConsts(&_elasticConsts[iQuad], _paramsCell[iQuad],
+ totalStrain[iQuad]);
return _elasticConsts;
} // calcDerivElastic
@@ -134,38 +104,43 @@
void
pylith::materials::ElasticMaterial::initCellData(const int numQuadPts)
{ // initCellData
- int size = numQuadPts;
- delete[] _density; _density = (size > 0) ? new double[size] : 0;
+ if (_numQuadPts != numQuadPts) {
+ _numQuadPts = numQuadPts;
- size = numQuadPts * stressSize();
- delete[] _stress; _stress = (size > 0) ? new double[size] : 0;
-
- size = numQuadPts * numElasticConsts();
- delete[] _elasticConsts; _elasticConsts = (size > 0) ? new double[size] : 0;
+ _paramsCell.resize(numQuadPts);
+ _density.resize(numQuadPts);
+ _stress.resize(numQuadPts);
+ _elasticConsts.resize(numQuadPts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ _paramsCell[iQuad].resize(_numParameters());
+ _density[iQuad].resize(1);
+ _stress[iQuad].resize(_tensorSize());
+ _elasticConsts[iQuad].resize(_numElasticConsts());
+ } // for
+ } // if
} // initCellData
// ----------------------------------------------------------------------
// Get parameters for cell.
void
-pylith::materials::ElasticMaterial::_getParameters(double** paramsCell,
- const Mesh::point_type& cell,
- const int numQuadPts)
+pylith::materials::ElasticMaterial::_getParameters(const Mesh::point_type& cell)
{ // _getParameters
+ const int numQuadPts = _numQuadPts;
+ assert(_paramsCell.size() == numQuadPts);
+
const int numParams = _numParameters();
const char** paramNames = _parameterNames();
-
- const int size = numParams * numQuadPts;
- delete[] *paramsCell; *paramsCell = (size > 0) ? new double[size] : 0;
+
for (int iParam=0; iParam < numParams; ++iParam) {
const ALE::Obj<real_section_type> parameter =
_parameters->getReal(paramNames[iParam]);
- assert(numQuadPts == parameter->getFiberDimension(cell));
+ assert(parameter->getFiberDimension(cell) == numQuadPts);
const real_section_type::value_type* parameterCell =
parameter->restrictPoint(cell);
for (int iQuadPt=0; iQuadPt < numQuadPts; ++iQuadPt)
- (*paramsCell)[iQuadPt*numParams+iParam] = parameterCell[iQuadPt];
+ _paramsCell[iQuadPt][iParam] = parameterCell[iQuadPt];
} // for
} // _getParameters
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh 2007-04-10 01:28:52 UTC (rev 6533)
@@ -23,6 +23,8 @@
#include "Material.hh" // ISA Material
#include <petscmesh.h> // USES Mesh
+#include "pylith/utils/stlfwd.hh" // USES double_array
+#include <vector> // USES std::vector
/// Namespace for pylith package
namespace pylith {
@@ -64,18 +66,15 @@
/** Compute density for cell at quadrature points.
*
* @param cell Finite-element cell
- * @param patch Finite-element patch
- * @param numQuadPts Number of quadrature points (consistency check)
*
* @returns Array of density values at cell's quadrature points.
*/
- const double* calcDensity(const Mesh::point_type& cell,
- const int numQuadPts);
+ const std::vector<double_array>&
+ calcDensity(const Mesh::point_type& cell);
/** Get stress tensor at quadrature points.
*
- * Index into array of elasticity constants:
- * index = iQuadPt*STRESSSIZE + iStress
+ * Size of array of stress tensors = [numQuadPts][tensorSize].
*
* Order of stresses for 3-D:
* 0: S11, 1: S22, 2: S33, 3: S12, 4: S23, 5: S13,
@@ -87,33 +86,18 @@
* 0: S11
*
* @param cell Finite-element cell
- * @param patch Finite-element patch
* @param totalStrain Total strain tensor at quadrature points
- * @param numQuadPts Number of quadrature points (consistency check)
- * @param spaceDim Spatial dimension (consistency check)
+ * [numQuadPts][tensorSize]
*
* @returns Array of stresses at cell's quadrature points.
*/
- const double* calcStress(const Mesh::point_type& cell,
- const double* totalStrain,
- const int numQuadPts,
- const int spaceDim);
+ const std::vector<double_array>&
+ calcStress(const Mesh::point_type& cell,
+ const std::vector<double_array>& totalStrain);
- /** Get number of entries in stress tensor for material.
- *
- * 1-D = 1
- * 2-D = 3
- * 3-D = 6
- *
- * @returns Number of entries in stress tensor
- */
- virtual
- int stressSize(void) const = 0;
-
/** Compute derivative of elasticity matrix for cell at quadrature points.
*
- * Index into array of elasticity constants:
- * index = iQuadPt*NUMELASTCONSTS + iConstant
+ * Size of array of elasticity constants = [numQuadPts][numElasticConsts]
*
* Order of elasticity constants for 3-D:
* 0: C1111, 1: C1122, 2: C1133, 3: C1112, 4: C1123, 5: C1113,
@@ -132,27 +116,13 @@
* 0: C1111
*
* @param cell Finite-element cell
- * @param patch Finite-element patch
* @param totalStrain Total strain tensor at quadrature points
- * @param numQuadPts Number of quadrature points (consistency check)
- * @param spaceDim Spatial dimension (consistency check)
+ * [numQuadPts][tensorSize]
*/
- const double* calcDerivElastic(const Mesh::point_type& cell,
- const double* totalStrain,
- const int numQuadPts,
- const int spaceDim);
+ const std::vector<double_array>&
+ calcDerivElastic(const Mesh::point_type& cell,
+ const std::vector<double_array>& totalStrain);
- /** Get number of entries in derivatives of elasticity matrix for material.
- *
- * 1-D = 1
- * 2-D = 6
- * 3-D = 21
- *
- * @returns Number of entries in derivative of elasticity matrix
- */
- virtual
- int numElasticConsts(void) const = 0;
-
/** Initialize arrays holding cell data.
*
* @param numQuadPts Number of quadrature points
@@ -168,52 +138,50 @@
*/
ElasticMaterial(const ElasticMaterial& m);
+ /** Get number of entries in stress/strain tensors.
+ *
+ * @returns Size of stress/strain tensors.
+ */
+ virtual
+ int _tensorSize(void) const = 0;
+
+ /** Get number of entries in elastic constants array.
+ *
+ * @returns Number of entries in array of elastic constants.
+ */
+ virtual
+ int _numElasticConsts(void) const = 0;
+
/** Compute density from parameters.
*
* @param density Array for density
- * @param size Size of array for density
* @param parameters Parameters at location
- * @param numParameters Number of parameters
*/
virtual
- void _calcDensity(double* const density,
- const int size,
- const double* parameters,
- const int numParameters) = 0;
+ void _calcDensity(double_array* const density,
+ const double_array& parameters) = 0;
/** Compute stress tensor from parameters.
*
* @param stress Array for stress tensor
- * @param size Size of array for stress tensor
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
virtual
- void _calcStress(double* const stress,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim) = 0;
+ void _calcStress(double_array* const stress,
+ const double_array& parameters,
+ const double_array& totalStrain) = 0;
/** Compute derivatives of elasticity matrix from parameters.
*
* @param elasticConsts Array for elastic constants
- * @param size Size of array
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
virtual
- void _calcElasticConsts(double* const elasticConsts,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim) = 0;
+ void _calcElasticConsts(double_array* const elasticConsts,
+ const double_array& parameters,
+ const double_array& totalStrain) = 0;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
@@ -226,42 +194,41 @@
/** Get parameters for cell.
*
- * Parameters are returned in paramsCells.
- *
- * size = numQuadPts * numParams
- * index = iQuad*numParams + iParam
- *
- * @param paramsCell Array of parameters for cell
* @param cell Finite-element cell
- * @param numQuadPts Number of quadrature points (consistency check)
*/
- void _getParameters(double** paramsCells,
- const Mesh::point_type& cell,
- const int numQuadPts);
+ void _getParameters(const Mesh::point_type& cell);
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
+ int _numQuadPts; ///< Number of quadrature points
+
+ /** Parameters at quadrature points for current cell.
+ *
+ * size = [numQuadPts][numParams]
+ */
+ std::vector<double_array> _paramsCell;
+
/** Density value at quadrature points for current cell.
*
- * size = numQuadPts
- * index = iQuadPt
+ * size = [numQuadPts][1]
+ * index = [iQuadPt][0]
*/
- double* _density;
+ std::vector<double_array> _density;
/** Stress tensor at quadrature points for current cell.
*
- * size = stressSize*numQuadPts
- * index = iQuadPt*stressSize+iStress
+ * size = [numQuadPts][tensorSize]
+ * index = [iQuadPt][iStress]
*/
- double* _stress;
+ std::vector<double_array> _stress;
/** Elasticity matrix at quadrature points for current cell.
*
- * size = numElasticConsts*numQuadPts
- * index = iQuadPt*numElasticConsts+iConstant
+ * size = [numQuadPts][numElasticConsts]
+ * index = [iQuadPt][iConstant]
*/
- double* _elasticConsts;
+ std::vector<double_array> _elasticConsts;
}; // class ElasticMaterial
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc 2007-04-10 01:28:52 UTC (rev 6533)
@@ -14,6 +14,7 @@
#include "ElasticPlaneStrain.hh" // implementation of object methods
+#include <valarray> // USES std::valarray (double_array)
#include <assert.h> // USES assert()
// ----------------------------------------------------------------------
@@ -30,7 +31,7 @@
class pylith::materials::_ElasticPlaneStrain {
public:
/// Number of entries in stress tensor.
- static const int stressSize;
+ static const int tensorSize;
/// Number of elastic constants (for general 3-D elastic material)
static const int numElasticConsts;
@@ -54,7 +55,7 @@
static const int pidLambda;
}; // _ElasticPlaneStrain
-const int pylith::materials::_ElasticPlaneStrain::stressSize = 3;
+const int pylith::materials::_ElasticPlaneStrain::tensorSize = 3;
const int pylith::materials::_ElasticPlaneStrain::numElasticConsts = 6;
const int pylith::materials::_ElasticPlaneStrain::numDBValues = 3;
const char* pylith::materials::_ElasticPlaneStrain::namesDBValues[] =
@@ -91,22 +92,6 @@
} // copy constructor
// ----------------------------------------------------------------------
-// Get number of entries in stress tensor.
-int
-pylith::materials::ElasticPlaneStrain::stressSize(void) const
-{ // stressSize
- return _ElasticPlaneStrain::stressSize;
-} // stressSize
-
-// ----------------------------------------------------------------------
-// Get number of elastic constants for material.
-int
-pylith::materials::ElasticPlaneStrain::numElasticConsts(void) const
-{ // numElasticConsts
- return _ElasticPlaneStrain::numElasticConsts;
-} // numElasticConsts
-
-// ----------------------------------------------------------------------
// Get names of values expected to be in database of parameters for
const char**
pylith::materials::ElasticPlaneStrain::_dbValues(void) const
@@ -141,15 +126,14 @@
// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
-pylith::materials::ElasticPlaneStrain::_dbToParameters(double* paramVals,
- const int numParams,
- const double* dbValues,
- const int numValues) const
+pylith::materials::ElasticPlaneStrain::_dbToParameters(double_array* const paramVals,
+ const double_array& dbValues) const
{ // computeParameters
assert(0 != paramVals);
+ const int numParams = paramVals->size();
assert(_ElasticPlaneStrain::numParameters == numParams);
- assert(0 != dbValues);
- assert(_ElasticPlaneStrain::numDBValues == numValues);
+ const int numDBValues = dbValues.size();
+ assert(_ElasticPlaneStrain::numDBValues == numDBValues);
const double density = dbValues[_ElasticPlaneStrain::didDensity];
const double vs = dbValues[_ElasticPlaneStrain::didVs];
@@ -158,41 +142,51 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- paramVals[_ElasticPlaneStrain::pidDensity] = density;
- paramVals[_ElasticPlaneStrain::pidMu] = mu;
- paramVals[_ElasticPlaneStrain::pidLambda] = lambda;
+ (*paramVals)[_ElasticPlaneStrain::pidDensity] = density;
+ (*paramVals)[_ElasticPlaneStrain::pidMu] = mu;
+ (*paramVals)[_ElasticPlaneStrain::pidLambda] = lambda;
} // computeParameters
// ----------------------------------------------------------------------
+// Get number of entries in stress tensor.
+int
+pylith::materials::ElasticPlaneStrain::_tensorSize(void) const
+{ // _tensorSize
+ return _ElasticPlaneStrain::tensorSize;
+} // _tensorSize
+
+// ----------------------------------------------------------------------
+// Get number of elastic constants for material.
+int
+pylith::materials::ElasticPlaneStrain::_numElasticConsts(void) const
+{ // _numElasticConsts
+ return _ElasticPlaneStrain::numElasticConsts;
+} // _numElasticConsts
+
+// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::ElasticPlaneStrain::_calcDensity(double* const density,
- const int size,
- const double* parameters,
- const int numParameters)
+pylith::materials::ElasticPlaneStrain::_calcDensity(double_array* const density,
+ const double_array& parameters)
{ // calcDensity
assert(0 != density);
- assert(1 == size);
- assert(0 != parameters);
+ assert(1 == density->size());
+ assert(_ElasticPlaneStrain::numParameters == parameters.size());
- *density = parameters[_ElasticPlaneStrain::pidDensity];
+ (*density)[0] = parameters[_ElasticPlaneStrain::pidDensity];
} // calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor at location from parameters.
void
-pylith::materials::ElasticPlaneStrain::_calcStress(double* const stress,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim)
+pylith::materials::ElasticPlaneStrain::_calcStress(double_array* const stress,
+ const double_array& parameters,
+ const double_array& totalStrain)
{ // _calcStress
assert(0 != stress);
- assert(_ElasticPlaneStrain::stressSize == size);
- assert(0 != parameters);
- assert(_ElasticPlaneStrain::numParameters == numParameters);
- assert(spaceDim == _dimension);
+ assert(_ElasticPlaneStrain::tensorSize == stress->size());
+ assert(_ElasticPlaneStrain::numParameters == parameters.size());
+ assert(_ElasticPlaneStrain::tensorSize == totalStrain.size());
const double density = parameters[_ElasticPlaneStrain::pidDensity];
const double mu = parameters[_ElasticPlaneStrain::pidMu];
@@ -206,25 +200,22 @@
const double s12 = lambda * (e11 + e22);
- stress[0] = s12 + 2.0*mu*e11;
- stress[1] = s12 + 2.0*mu*e22;
- stress[2] = 2.0 * mu * e12;
+ (*stress)[0] = s12 + 2.0*mu*e11;
+ (*stress)[1] = s12 + 2.0*mu*e22;
+ (*stress)[2] = 2.0 * mu * e12;
} // _calcStress
// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::ElasticPlaneStrain::_calcElasticConsts(double* const elasticConsts,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim)
+pylith::materials::ElasticPlaneStrain::_calcElasticConsts(double_array* const elasticConsts,
+ const double_array& parameters,
+ const double_array& totalStrain)
{ // calcElasticConsts
assert(0 != elasticConsts);
- assert(_ElasticPlaneStrain::numElasticConsts == size);
- assert(0 != parameters);
- assert(_ElasticPlaneStrain::numParameters == numParameters);
+ assert(_ElasticPlaneStrain::numElasticConsts == elasticConsts->size());
+ assert(_ElasticPlaneStrain::numParameters == parameters.size());
+ assert(_ElasticPlaneStrain::tensorSize == totalStrain.size());
const double density = parameters[_ElasticPlaneStrain::pidDensity];
const double mu = parameters[_ElasticPlaneStrain::pidMu];
@@ -232,12 +223,12 @@
const double lambda2mu = lambda + 2.0*mu;
- elasticConsts[0] = lambda2mu; // C1111
- elasticConsts[1] = lambda; // C1122
- elasticConsts[2] = 0; // C1112
- elasticConsts[3] = lambda2mu; // C2222
- elasticConsts[4] = 0; // C2212
- elasticConsts[5] = 2.0*mu; // C1212
+ (*elasticConsts)[0] = lambda2mu; // C1111
+ (*elasticConsts)[1] = lambda; // C1122
+ (*elasticConsts)[2] = 0; // C1112
+ (*elasticConsts)[3] = lambda2mu; // C2222
+ (*elasticConsts)[4] = 0; // C2212
+ (*elasticConsts)[5] = 2.0*mu; // C1212
} // calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh 2007-04-10 01:28:52 UTC (rev 6533)
@@ -55,26 +55,6 @@
*/
ElasticMaterial* clone(void) const;
- /** Get number of entries in stress tensor.
- *
- * 1-D = 1
- * 2-D = 3
- * 3-D = 6
- *
- * @returns Number of entries in stress tensor.
- */
- int stressSize(void) const;
-
- /** Get number of elastic constants for material.
- *
- * 1-D = 1
- * 2-D = 6
- * 3-D = 21
- *
- * @returns Number of elastic constants
- */
- int numElasticConsts(void) const;
-
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -116,58 +96,58 @@
* parameterNames().
*
* @param paramVals Array of parameters
- * @param numParams Number of parameters
* @param dbValues Array of database values
- * @param numValues Number of database values
*/
- void _dbToParameters(double* paramVals,
- const int numParams,
- const double* dbValues,
- const int numValues) const;
+ void _dbToParameters(double_array* const paramVals,
+ const double_array& dbValues) const;
+ /** Get number of entries in stress tensor.
+ *
+ * 1-D = 1
+ * 2-D = 3
+ * 3-D = 6
+ *
+ * @returns Number of entries in stress tensor.
+ */
+ int _tensorSize(void) const;
+
+ /** Get number of elastic constants for material.
+ *
+ * 1-D = 1
+ * 2-D = 6
+ * 3-D = 21
+ *
+ * @returns Number of elastic constants
+ */
+ int _numElasticConsts(void) const;
+
/** Compute density from parameters.
*
* @param density Array for density
- * @param size Size of array for density
* @param parameters Parameters at location
- * @param numParameters Number of parameters
*/
- void _calcDensity(double* const density,
- const int size,
- const double* parameters,
- const int numParameters);
+ void _calcDensity(double_array* const density,
+ const double_array& parameters);
/** Compute stress tensor from parameters.
*
* @param stress Array for stress tensor
- * @param size Size of array for stress tensor
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
- void _calcStress(double* const stress,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim);
+ void _calcStress(double_array* const stress,
+ const double_array& parameters,
+ const double_array& totalStrain);
/** Compute derivatives of elasticity matrix from parameters.
*
* @param elasticConsts Array for elastic constants
- * @param size Size of array
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
- void _calcElasticConsts(double* const elasticConsts,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim);
+ void _calcElasticConsts(double_array* const elasticConsts,
+ const double_array& parameters,
+ const double_array& totalStrain);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc 2007-04-10 01:28:52 UTC (rev 6533)
@@ -14,6 +14,7 @@
#include "ElasticPlaneStress.hh" // implementation of object methods
+#include <valarray> // USES std::valarray (double_array)
#include <assert.h> // USES assert()
// ----------------------------------------------------------------------
@@ -30,7 +31,7 @@
class pylith::materials::_ElasticPlaneStress {
public:
/// Number of entries in stress tensor.
- static const int stressSize;
+ static const int tensorSize;
/// Number of elastic constants (for general 3-D elastic material)
static const int numElasticConsts;
@@ -54,7 +55,7 @@
static const int pidLambda;
}; // _ElasticPlaneStress
-const int pylith::materials::_ElasticPlaneStress::stressSize = 3;
+const int pylith::materials::_ElasticPlaneStress::tensorSize = 3;
const int pylith::materials::_ElasticPlaneStress::numElasticConsts = 6;
const int pylith::materials::_ElasticPlaneStress::numDBValues = 3;
const char* pylith::materials::_ElasticPlaneStress::namesDBValues[] =
@@ -91,22 +92,6 @@
} // copy constructor
// ----------------------------------------------------------------------
-// Get number of entries in stress tensor.
-int
-pylith::materials::ElasticPlaneStress::stressSize(void) const
-{ // stressSize
- return _ElasticPlaneStress::stressSize;
-} // stressSize
-
-// ----------------------------------------------------------------------
-// Get number of elastic constants for material.
-int
-pylith::materials::ElasticPlaneStress::numElasticConsts(void) const
-{ // numElasticConsts
- return _ElasticPlaneStress::numElasticConsts;
-} // numElasticConsts
-
-// ----------------------------------------------------------------------
// Get names of values expected to be in database of parameters for
const char**
pylith::materials::ElasticPlaneStress::_dbValues(void) const
@@ -141,15 +126,14 @@
// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
-pylith::materials::ElasticPlaneStress::_dbToParameters(double* paramVals,
- const int numParams,
- const double* dbValues,
- const int numValues) const
-{ // computeParameters
+pylith::materials::ElasticPlaneStress::_dbToParameters(double_array* paramVals,
+ const double_array& dbValues) const
+{ // _dbToParameters
assert(0 != paramVals);
+ const int numParams = paramVals->size();
assert(_ElasticPlaneStress::numParameters == numParams);
- assert(0 != dbValues);
- assert(_ElasticPlaneStress::numDBValues == numValues);
+ const int numDBValues = dbValues.size();
+ assert(_ElasticPlaneStress::numDBValues == numDBValues);
const double density = dbValues[_ElasticPlaneStress::didDensity];
const double vs = dbValues[_ElasticPlaneStress::didVs];
@@ -158,41 +142,51 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- paramVals[_ElasticPlaneStress::pidDensity] = density;
- paramVals[_ElasticPlaneStress::pidMu] = mu;
- paramVals[_ElasticPlaneStress::pidLambda] = lambda;
-} // computeParameters
+ (*paramVals)[_ElasticPlaneStress::pidDensity] = density;
+ (*paramVals)[_ElasticPlaneStress::pidMu] = mu;
+ (*paramVals)[_ElasticPlaneStress::pidLambda] = lambda;
+} // _dbToParameters
// ----------------------------------------------------------------------
+// Get number of entries in stress tensor.
+int
+pylith::materials::ElasticPlaneStress::_tensorSize(void) const
+{ // _tensorSize
+ return _ElasticPlaneStress::tensorSize;
+} // _tensorSize
+
+// ----------------------------------------------------------------------
+// Get number of elastic constants for material.
+int
+pylith::materials::ElasticPlaneStress::_numElasticConsts(void) const
+{ // _numElasticConsts
+ return _ElasticPlaneStress::numElasticConsts;
+} // _numElasticConsts
+
+// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::ElasticPlaneStress::_calcDensity(double* const density,
- const int size,
- const double* parameters,
- const int numParameters)
+pylith::materials::ElasticPlaneStress::_calcDensity(double_array* const density,
+ const double_array& parameters)
{ // calcDensity
assert(0 != density);
- assert(1 == size);
- assert(0 != parameters);
+ assert(1 == density->size());
+ assert(_ElasticPlaneStress::numParameters == parameters.size());
- *density = parameters[_ElasticPlaneStress::pidDensity];
+ (*density)[0] = parameters[_ElasticPlaneStress::pidDensity];
} // calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor at location from parameters.
void
-pylith::materials::ElasticPlaneStress::_calcStress(double* const stress,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim)
+pylith::materials::ElasticPlaneStress::_calcStress(double_array* const stress,
+ const double_array& parameters,
+ const double_array& totalStrain)
{ // _calcStress
assert(0 != stress);
- assert(_ElasticPlaneStress::stressSize == size);
- assert(0 != parameters);
- assert(_ElasticPlaneStress::numParameters == numParameters);
- assert(spaceDim == _dimension);
+ assert(_ElasticPlaneStress::tensorSize == stress->size());
+ assert(_ElasticPlaneStress::numParameters == parameters.size());
+ assert(_ElasticPlaneStress::tensorSize == totalStrain.size());
const double density = parameters[_ElasticPlaneStress::pidDensity];
const double mu = parameters[_ElasticPlaneStress::pidMu];
@@ -205,25 +199,22 @@
const double e22 = totalStrain[1];
const double e12 = totalStrain[2];
- stress[0] = (4.0*mu*lambdamu * e11 + 2.0*mu*lambda * e22)/lambda2mu;
- stress[1] = (2.0*mu*lambda * e11 + 4.0*mu*lambdamu * e22)/lambda2mu;
- stress[2] = 2.0 * mu * e12;
+ (*stress)[0] = (4.0*mu*lambdamu * e11 + 2.0*mu*lambda * e22)/lambda2mu;
+ (*stress)[1] = (2.0*mu*lambda * e11 + 4.0*mu*lambdamu * e22)/lambda2mu;
+ (*stress)[2] = 2.0 * mu * e12;
} // _calcStress
// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::ElasticPlaneStress::_calcElasticConsts(double* const elasticConsts,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim)
+pylith::materials::ElasticPlaneStress::_calcElasticConsts(double_array* const elasticConsts,
+ const double_array& parameters,
+ const double_array& totalStrain)
{ // calcElasticConsts
assert(0 != elasticConsts);
- assert(_ElasticPlaneStress::numElasticConsts == size);
- assert(0 != parameters);
- assert(_ElasticPlaneStress::numParameters == numParameters);
+ assert(_ElasticPlaneStress::numElasticConsts == elasticConsts->size());
+ assert(_ElasticPlaneStress::numParameters == parameters.size());
+ assert(_ElasticPlaneStress::tensorSize == totalStrain.size());
const double density = parameters[_ElasticPlaneStress::pidDensity];
const double mu = parameters[_ElasticPlaneStress::pidMu];
@@ -232,12 +223,12 @@
const double lambda2mu = lambda + 2.0 * mu;
const double c11 = 4.0 * mu * (lambda + mu) / lambda2mu;
- elasticConsts[0] = c11; // C1111
- elasticConsts[1] = 2.0 * mu * lambda / lambda2mu; // C1122
- elasticConsts[2] = 0; // C1112
- elasticConsts[3] = c11; // C2222
- elasticConsts[4] = 0; // C2212
- elasticConsts[5] = 2.0 * mu; // C1212
+ (*elasticConsts)[0] = c11; // C1111
+ (*elasticConsts)[1] = 2.0 * mu * lambda / lambda2mu; // C1122
+ (*elasticConsts)[2] = 0; // C1112
+ (*elasticConsts)[3] = c11; // C2222
+ (*elasticConsts)[4] = 0; // C2212
+ (*elasticConsts)[5] = 2.0 * mu; // C1212
} // calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh 2007-04-10 01:28:52 UTC (rev 6533)
@@ -55,26 +55,6 @@
*/
ElasticMaterial* clone(void) const;
- /** Get number of entries in stress tensor.
- *
- * 1-D = 1
- * 2-D = 3
- * 3-D = 6
- *
- * @returns Number of entries in stress tensor.
- */
- int stressSize(void) const;
-
- /** Get number of elastic constants for material.
- *
- * 1-D = 1
- * 2-D = 6
- * 3-D = 21
- *
- * @returns Number of elastic constants
- */
- int numElasticConsts(void) const;
-
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -116,58 +96,58 @@
* parameterNames().
*
* @param paramVals Array of parameters
- * @param numParams Number of parameters
* @param dbValues Array of database values
- * @param numValues Number of database values
*/
- void _dbToParameters(double* paramVals,
- const int numParams,
- const double* dbValues,
- const int numValues) const;
+ void _dbToParameters(double_array* const paramVals,
+ const double_array& dbValues) const;
+ /** Get number of entries in stress tensor.
+ *
+ * 1-D = 1
+ * 2-D = 3
+ * 3-D = 6
+ *
+ * @returns Number of entries in stress tensor.
+ */
+ int _tensorSize(void) const;
+
+ /** Get number of elastic constants for material.
+ *
+ * 1-D = 1
+ * 2-D = 6
+ * 3-D = 21
+ *
+ * @returns Number of elastic constants
+ */
+ int _numElasticConsts(void) const;
+
/** Compute density from parameters.
*
* @param density Array for density
- * @param size Size of array for density
* @param parameters Parameters at location
- * @param numParameters Number of parameters
*/
- void _calcDensity(double* const density,
- const int size,
- const double* parameters,
- const int numParameters);
+ void _calcDensity(double_array* const density,
+ const double_array& parameters);
/** Compute stress tensor from parameters.
*
* @param stress Array for stress tensor
- * @param size Size of array for stress tensor
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
- void _calcStress(double* const stress,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim);
+ void _calcStress(double_array* const stress,
+ const double_array& parameters,
+ const double_array& totalStrain);
/** Compute derivatives of elasticity matrix from parameters.
*
* @param elasticConsts Array for elastic constants
- * @param size Size of array
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
- void _calcElasticConsts(double* const elasticConsts,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim);
+ void _calcElasticConsts(double_array* const elasticConsts,
+ const double_array& parameters,
+ const double_array& totalStrain);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc 2007-04-10 01:28:52 UTC (rev 6533)
@@ -14,6 +14,7 @@
#include "ElasticStrain1D.hh" // implementation of object methods
+#include <valarray> // USES std::valarray (double_array)
#include <assert.h> // USES assert()
// ----------------------------------------------------------------------
@@ -30,7 +31,7 @@
class pylith::materials::_ElasticStrain1D {
public:
/// Number of entries in stress tensor.
- static const int stressSize;
+ static const int tensorSize;
/// Number of entries in derivative of elasticity matrix.
static const int numElasticConsts;
@@ -52,7 +53,7 @@
static const int pidLambda2mu;
}; // _ElasticStrain1D
-const int pylith::materials::_ElasticStrain1D::stressSize = 1;
+const int pylith::materials::_ElasticStrain1D::tensorSize = 1;
const int pylith::materials::_ElasticStrain1D::numElasticConsts = 1;
const int pylith::materials::_ElasticStrain1D::numDBValues = 2;
const char* pylith::materials::_ElasticStrain1D::namesDBValues[] =
@@ -87,22 +88,6 @@
} // copy constructor
// ----------------------------------------------------------------------
-// Get number of entries in stress tensor.
-int
-pylith::materials::ElasticStrain1D::stressSize(void) const
-{ // stressSize
- return _ElasticStrain1D::stressSize;
-} // stressSize
-
-// ----------------------------------------------------------------------
-// Get number of elastic constants for material.
-int
-pylith::materials::ElasticStrain1D::numElasticConsts(void) const
-{ // numElasticConsts
- return _ElasticStrain1D::numElasticConsts;
-} // numElasticConsts
-
-// ----------------------------------------------------------------------
// Get names of values expected to be in database of parameters for
const char**
pylith::materials::ElasticStrain1D::_dbValues(void) const
@@ -137,82 +122,88 @@
// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
-pylith::materials::ElasticStrain1D::_dbToParameters(double* paramVals,
- const int numParams,
- const double* dbValues,
- const int numValues) const
-{ // computeParameters
+pylith::materials::ElasticStrain1D::_dbToParameters(double_array* const paramVals,
+ const double_array& dbValues) const
+{ // _dbToParameters
assert(0 != paramVals);
+ const int numParams = paramVals->size();
assert(_ElasticStrain1D::numParameters == numParams);
- assert(0 != dbValues);
- assert(_ElasticStrain1D::numDBValues == numValues);
+ const int numDBValues = dbValues.size();
+ assert(_ElasticStrain1D::numDBValues == numDBValues);
const double density = dbValues[_ElasticStrain1D::didDensity];
const double vp = dbValues[_ElasticStrain1D::didVp];
const double lambda2mu = density * vp*vp;
- paramVals[_ElasticStrain1D::pidDensity] = density;
- paramVals[_ElasticStrain1D::pidLambda2mu] = lambda2mu;
-} // computeParameters
+ (*paramVals)[_ElasticStrain1D::pidDensity] = density;
+ (*paramVals)[_ElasticStrain1D::pidLambda2mu] = lambda2mu;
+} // _dbToParameters
// ----------------------------------------------------------------------
+// Get number of entries in stress tensor.
+int
+pylith::materials::ElasticStrain1D::_tensorSize(void) const
+{ // _tensorSize
+ return _ElasticStrain1D::tensorSize;
+} // _tensorSize
+
+// ----------------------------------------------------------------------
+// Get number of elastic constants for material.
+int
+pylith::materials::ElasticStrain1D::_numElasticConsts(void) const
+{ // _numElasticConsts
+ return _ElasticStrain1D::numElasticConsts;
+} // _numElasticConsts
+
+// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::ElasticStrain1D::_calcDensity(double* const density,
- const int size,
- const double* parameters,
- const int numParameters)
+pylith::materials::ElasticStrain1D::_calcDensity(double_array* const density,
+ const double_array& parameters)
{ // _calcDensity
assert(0 != density);
- assert(1 == size);
- assert(0 != parameters);
+ assert(1 == density->size());
+ assert(_ElasticStrain1D::numParameters == parameters.size());
- *density = parameters[_ElasticStrain1D::pidDensity];
+ (*density)[0] = parameters[_ElasticStrain1D::pidDensity];
} // _calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor at location from parameters.
void
-pylith::materials::ElasticStrain1D::_calcStress(double* const stress,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim)
+pylith::materials::ElasticStrain1D::_calcStress(double_array* const stress,
+ const double_array& parameters,
+ const double_array& totalStrain)
{ // _calcStress
assert(0 != stress);
- assert(_ElasticStrain1D::stressSize == size);
- assert(0 != parameters);
- assert(_ElasticStrain1D::numParameters == numParameters);
- assert(spaceDim == _dimension);
+ assert(_ElasticStrain1D::tensorSize == stress->size());
+ assert(_ElasticStrain1D::numParameters == parameters.size());
+ assert(_ElasticStrain1D::tensorSize == totalStrain.size());
const double density = parameters[_ElasticStrain1D::pidDensity];
const double lambda2mu = parameters[_ElasticStrain1D::pidLambda2mu];
const double e11 = totalStrain[0];
- stress[0] = lambda2mu * e11;
+ (*stress)[0] = lambda2mu * e11;
} // _calcStress
// ----------------------------------------------------------------------
// Compute derivative of elasticity matrix at location from parameters.
void
-pylith::materials::ElasticStrain1D::_calcElasticConsts(double* const elasticConsts,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim)
+pylith::materials::ElasticStrain1D::_calcElasticConsts(
+ double_array* const elasticConsts,
+ const double_array& parameters,
+ const double_array& totalStrain)
{ // _calcElasticConsts
assert(0 != elasticConsts);
- assert(_ElasticStrain1D::numElasticConsts == size);
- assert(0 != parameters);
- assert(_ElasticStrain1D::numParameters == numParameters);
- assert(spaceDim == _dimension);
+ assert(_ElasticStrain1D::numElasticConsts == elasticConsts->size());
+ assert(_ElasticStrain1D::numParameters == parameters.size());
+ assert(_ElasticStrain1D::tensorSize == totalStrain.size());
const double density = parameters[_ElasticStrain1D::pidDensity];
const double lambda2mu = parameters[_ElasticStrain1D::pidLambda2mu];
- elasticConsts[0] = lambda2mu;
+ (*elasticConsts)[0] = lambda2mu;
} // _calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh 2007-04-10 01:28:52 UTC (rev 6533)
@@ -54,26 +54,6 @@
*/
ElasticMaterial* clone(void) const;
- /** Get number of entries in stress tensor.
- *
- * 1-D = 1
- * 2-D = 3
- * 3-D = 6
- *
- * @returns Number of entries in stress tensor.
- */
- int stressSize(void) const;
-
- /** Get number of elastic constants for material.
- *
- * 1-D = 1
- * 2-D = 6
- * 3-D = 21
- *
- * @returns Number of elastic constants
- */
- int numElasticConsts(void) const;
-
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -115,58 +95,58 @@
* parameterNames().
*
* @param paramVals Array of parameters
- * @param numParams Number of parameters
* @param dbValues Array of database values
- * @param numValues Number of database values
*/
- void _dbToParameters(double* paramVals,
- const int numParams,
- const double* dbValues,
- const int numValues) const;
+ void _dbToParameters(double_array* const paramVals,
+ const double_array& dbValues) const;
+ /** Get number of entries in stress tensor.
+ *
+ * 1-D = 1
+ * 2-D = 3
+ * 3-D = 6
+ *
+ * @returns Number of entries in stress tensor.
+ */
+ int _tensorSize(void) const;
+
+ /** Get number of elastic constants for material.
+ *
+ * 1-D = 1
+ * 2-D = 6
+ * 3-D = 21
+ *
+ * @returns Number of elastic constants
+ */
+ int _numElasticConsts(void) const;
+
/** Compute density from parameters.
*
* @param density Array for density
- * @param size Size of array for density
* @param parameters Parameters at location
- * @param numParameters Number of parameters
*/
- void _calcDensity(double* const density,
- const int size,
- const double* parameters,
- const int numParameters);
+ void _calcDensity(double_array* const density,
+ const double_array& parameters);
/** Compute stress tensor from parameters.
*
* @param stress Array for stress tensor
- * @param size Size of array for stress tensor
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
- void _calcStress(double* const stress,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim);
+ void _calcStress(double_array* const stress,
+ const double_array& parameters,
+ const double_array& totalStrain);
/** Compute derivatives of elasticity matrix from parameters.
*
* @param elasticConsts Array for elastic constants
- * @param size Size of array
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
- void _calcElasticConsts(double* const elasticConsts,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim);
+ void _calcElasticConsts(double_array* const elasticConsts,
+ const double_array& parameters,
+ const double_array& totalStrain);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc 2007-04-10 01:28:52 UTC (rev 6533)
@@ -14,6 +14,7 @@
#include "ElasticStress1D.hh" // implementation of object methods
+#include <valarray> // USES std::valarray (double_array)
#include <assert.h> // USES assert()
// ----------------------------------------------------------------------
@@ -30,7 +31,7 @@
class pylith::materials::_ElasticStress1D {
public:
/// Number of entries in stress tensor.
- static const int stressSize;
+ static const int tensorSize;
/// Number of entries in derivative of elasticity matrix.
static const int numElasticConsts;
@@ -54,7 +55,7 @@
static const int pidLambda;
}; // _ElasticStress1D
-const int pylith::materials::_ElasticStress1D::stressSize = 1;
+const int pylith::materials::_ElasticStress1D::tensorSize = 1;
const int pylith::materials::_ElasticStress1D::numElasticConsts = 1;
const int pylith::materials::_ElasticStress1D::numDBValues = 3;
const char* pylith::materials::_ElasticStress1D::namesDBValues[] =
@@ -91,22 +92,6 @@
} // copy constructor
// ----------------------------------------------------------------------
-// Get number of entries in stress tensor.
-int
-pylith::materials::ElasticStress1D::stressSize(void) const
-{ // stressSize
- return _ElasticStress1D::stressSize;
-} // stressSize
-
-// ----------------------------------------------------------------------
-// Get number of elastic constants for material.
-int
-pylith::materials::ElasticStress1D::numElasticConsts(void) const
-{ // numElasticConsts
- return _ElasticStress1D::numElasticConsts;
-} // numElasticConsts
-
-// ----------------------------------------------------------------------
// Get names of values expected to be in database of parameters for
const char**
pylith::materials::ElasticStress1D::_dbValues(void) const
@@ -141,15 +126,14 @@
// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
-pylith::materials::ElasticStress1D::_dbToParameters(double* paramVals,
- const int numParams,
- const double* dbValues,
- const int numValues) const
+pylith::materials::ElasticStress1D::_dbToParameters(double_array* const paramVals,
+ const double_array& dbValues) const
{ // computeParameters
assert(0 != paramVals);
+ const int numParams = paramVals->size();
assert(_ElasticStress1D::numParameters == numParams);
- assert(0 != dbValues);
- assert(_ElasticStress1D::numDBValues == numValues);
+ const int numDBValues = dbValues.size();
+ assert(_ElasticStress1D::numDBValues == numDBValues);
const double density = dbValues[_ElasticStress1D::didDensity];
const double vs = dbValues[_ElasticStress1D::didVs];
@@ -158,71 +142,78 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- paramVals[_ElasticStress1D::pidDensity] = density;
- paramVals[_ElasticStress1D::pidMu] = mu;
- paramVals[_ElasticStress1D::pidLambda] = lambda;
+ (*paramVals)[_ElasticStress1D::pidDensity] = density;
+ (*paramVals)[_ElasticStress1D::pidMu] = mu;
+ (*paramVals)[_ElasticStress1D::pidLambda] = lambda;
} // computeParameters
// ----------------------------------------------------------------------
+// Get number of entries in stress tensor.
+int
+pylith::materials::ElasticStress1D::_tensorSize(void) const
+{ // _tensorSize
+ return _ElasticStress1D::tensorSize;
+} // _tensorSize
+
+// ----------------------------------------------------------------------
+// Get number of elastic constants for material.
+int
+pylith::materials::ElasticStress1D::_numElasticConsts(void) const
+{ // _numElasticConsts
+ return _ElasticStress1D::numElasticConsts;
+} // _numElasticConsts
+
+// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::ElasticStress1D::_calcDensity(double* const density,
- const int size,
- const double* parameters,
- const int numParameters)
+pylith::materials::ElasticStress1D::_calcDensity(double_array* const density,
+ const double_array& parameters)
{ // _calcDensity
assert(0 != density);
- assert(1 == size);
- assert(0 != parameters);
+ assert(1 == density->size());
+ assert(_ElasticStress1D::numParameters == parameters.size());
- *density = parameters[_ElasticStress1D::pidDensity];
+ (*density)[0] = parameters[_ElasticStress1D::pidDensity];
} // _calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor at location from parameters.
void
-pylith::materials::ElasticStress1D::_calcStress(double* const stress,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim)
+pylith::materials::ElasticStress1D::_calcStress(double_array* const stress,
+ const double_array& parameters,
+ const double_array& totalStrain)
{ // _calcStress
assert(0 != stress);
- assert(_ElasticStress1D::stressSize == size);
- assert(0 != parameters);
- assert(_ElasticStress1D::numParameters == numParameters);
- assert(spaceDim == _dimension);
+ assert(_ElasticStress1D::tensorSize == stress->size());
+ assert(_ElasticStress1D::numParameters == parameters.size());
+ assert(_ElasticStress1D::tensorSize == totalStrain.size());
const double density = parameters[_ElasticStress1D::pidDensity];
const double mu = parameters[_ElasticStress1D::pidMu];
const double lambda = parameters[_ElasticStress1D::pidLambda];
const double e11 = totalStrain[0];
- stress[0] = mu * (3.0*lambda+2.0*mu) / (lambda + mu) * e11;
+ (*stress)[0] = mu * (3.0*lambda+2.0*mu) / (lambda + mu) * e11;
} // _calcStress
// ----------------------------------------------------------------------
// Compute derivative of elasticity matrix at location from parameters.
void
-pylith::materials::ElasticStress1D::_calcElasticConsts(double* const elasticConsts,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim)
+pylith::materials::ElasticStress1D::_calcElasticConsts(
+ double_array* const elasticConsts,
+ const double_array& parameters,
+ const double_array& totalStrain)
{ // _calcElasticConsts
assert(0 != elasticConsts);
- assert(_ElasticStress1D::numElasticConsts == size);
- assert(0 != parameters);
- assert(_ElasticStress1D::numParameters == numParameters);
- assert(spaceDim == _dimension);
+ assert(_ElasticStress1D::numElasticConsts == elasticConsts->size());
+ assert(_ElasticStress1D::numParameters == parameters.size());
+ assert(_ElasticStress1D::tensorSize == totalStrain.size());
const double density = parameters[_ElasticStress1D::pidDensity];
const double mu = parameters[_ElasticStress1D::pidMu];
const double lambda = parameters[_ElasticStress1D::pidLambda];
- elasticConsts[0] = mu * (3.0*lambda+2.0*mu) / (lambda + mu);
+ (*elasticConsts)[0] = mu * (3.0*lambda+2.0*mu) / (lambda + mu);
} // _calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh 2007-04-10 01:28:52 UTC (rev 6533)
@@ -55,26 +55,6 @@
*/
ElasticMaterial* clone(void) const;
- /** Get number of entries in stress tensor.
- *
- * 1-D = 1
- * 2-D = 3
- * 3-D = 6
- *
- * @returns Number of entries in stress tensor.
- */
- int stressSize(void) const;
-
- /** Get number of elastic constants for material.
- *
- * 1-D = 1
- * 2-D = 6
- * 3-D = 21
- *
- * @returns Number of elastic constants
- */
- int numElasticConsts(void) const;
-
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -116,58 +96,58 @@
* parameterNames().
*
* @param paramVals Array of parameters
- * @param numParams Number of parameters
* @param dbValues Array of database values
- * @param numValues Number of database values
*/
- void _dbToParameters(double* paramVals,
- const int numParams,
- const double* dbValues,
- const int numValues) const;
+ void _dbToParameters(double_array* const paramVals,
+ const double_array& dbValues) const;
+ /** Get number of entries in stress tensor.
+ *
+ * 1-D = 1
+ * 2-D = 3
+ * 3-D = 6
+ *
+ * @returns Number of entries in stress tensor.
+ */
+ int _tensorSize(void) const;
+
+ /** Get number of elastic constants for material.
+ *
+ * 1-D = 1
+ * 2-D = 6
+ * 3-D = 21
+ *
+ * @returns Number of elastic constants
+ */
+ int _numElasticConsts(void) const;
+
/** Compute density from parameters.
*
* @param density Array for density
- * @param size Size of array for density
* @param parameters Parameters at location
- * @param numParameters Number of parameters
*/
- void _calcDensity(double* const density,
- const int size,
- const double* parameters,
- const int numParameters);
+ void _calcDensity(double_array* const density,
+ const double_array& parameters);
/** Compute stress tensor from parameters.
*
* @param stress Array for stress tensor
- * @param size Size of array for stress tensor
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
- void _calcStress(double* const stress,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim);
+ void _calcStress(double_array* const stress,
+ const double_array& parameters,
+ const double_array& totalStrain);
/** Compute derivatives of elasticity matrix from parameters.
*
* @param elasticConsts Array for elastic constants
- * @param size Size of array
* @param parameters Parameters at locations.
- * @param numParameters Number of parameters.
* @param totalStrain Total strain at locations.
- * @param spaceDim Spatial dimension for locations.
*/
- void _calcElasticConsts(double* const elasticConsts,
- const int size,
- const double* parameters,
- const int numParameters,
- const double* totalStrain,
- const int spaceDim);
+ void _calcElasticConsts(double_array* const elasticConsts,
+ const double_array& parameters,
+ const double_array& totalStrain);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.cc 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.cc 2007-04-10 01:28:52 UTC (rev 6533)
@@ -21,6 +21,9 @@
#include <petscmesh.h> // USES Mesh
+#include <vector> // USES std::vector
+#include <valarray> // USES std::valarray (double_array)
+
#include <assert.h> // USES assert()
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
@@ -92,8 +95,7 @@
const int numParams = _numParameters();
const char** paramNames = _parameterNames();
- ALE::Obj<real_section_type>* paramSections =
- (numParams > 0) ? new ALE::Obj<real_section_type>[numParams] : 0;
+ std::vector<ALE::Obj<real_section_type> > paramSections(numParams);
for (int iParam=0; iParam < numParams; ++iParam) {
_parameters->addReal(paramNames[iParam]);
@@ -108,11 +110,11 @@
_db->queryVals(_dbValues(), numValues);
// Loop over cells
- double* queryData = (numValues > 0) ? new double[numValues] : 0;
- double* paramData = (numParams > 0) ? new double[numParams] : 0;
- double** cellData = (numParams > 0) ? new double*[numParams] : 0;
+ double_array queryData(numValues);
+ double_array paramData(numParams);
+ std::vector<double_array> cellData(numParams);
for (int iParam = 0; iParam < numParams; ++iParam)
- cellData[iParam] = (numQuadPts > 0) ? new double[numQuadPts] : 0;
+ cellData[iParam] = double_array(numQuadPts);
for (ALE::Mesh::label_sequence::iterator cellIter=cells->begin();
cellIter != cellsEnd;
++cellIter) {
@@ -126,7 +128,7 @@
for (int iQuadPt=0, index=0;
iQuadPt < numQuadPts;
++iQuadPt, index+=spaceDim) {
- const int err = _db->query(queryData, numValues, &quadPts[index],
+ const int err = _db->query(&queryData[0], numValues, &quadPts[index],
spaceDim, cs);
if (err) {
std::ostringstream msg;
@@ -136,32 +138,17 @@
msg << " " << quadPts[index+spaceDim];
msg << ") in material " << _label << "\n"
<< "using spatial database " << _db->label() << ".";
-
- // Cleanup, then throw exception
- for (int iParam=0; iParam < numParams; ++iParam) {
- delete[] cellData[iParam]; cellData[iParam] = 0;
- } // for
- delete[] cellData; cellData = 0;
- delete[] queryData; queryData = 0;
- delete[] paramData; paramData = 0;
throw std::runtime_error(msg.str());
} // if
- _dbToParameters(paramData, numParams, queryData, numValues);
+ _dbToParameters(¶mData, queryData);
for (int iParam=0; iParam < numParams; ++iParam)
cellData[iParam][iQuadPt] = paramData[iParam];
} // for
// Assemble cell contribution into fields
for (int iParam=0; iParam < numParams; ++iParam)
- mesh->updateAdd(paramSections[iParam], *cellIter, cellData[iParam]);
+ mesh->updateAdd(paramSections[iParam], *cellIter, &cellData[iParam][0]);
} // for
- for (int iParam=0; iParam < numParams; ++iParam) {
- delete[] cellData[iParam]; cellData[iParam] = 0;
- } // for
- delete[] cellData; cellData = 0;
- delete[] queryData; queryData = 0;
- delete[] paramData; paramData = 0;
- delete[] paramSections; paramSections = 0;
// Close database
_db->close();
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.hh 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.hh 2007-04-10 01:28:52 UTC (rev 6533)
@@ -20,6 +20,7 @@
#if !defined(pylith_materials_material_hh)
#define pylith_materials_material_hh
+#include "pylith/utils/stlfwd.hh"
#include <string> // HASA std::string
/// Namespace for pylith package
@@ -160,10 +161,8 @@
* @param numValues Number of database values
*/
virtual
- void _dbToParameters(double* paramVals,
- const int numParams,
- const double* dbValues,
- const int numValues) const = 0;
+ void _dbToParameters(double_array* const paramVals,
+ const double_array& dbValues) const = 0;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/utils/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/utils/Makefile.am 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/utils/Makefile.am 2007-04-10 01:28:52 UTC (rev 6533)
@@ -14,7 +14,8 @@
include $(top_srcdir)/subpackage.am
subpkginclude_HEADERS = \
- petscfwd.h
+ petscfwd.h \
+ stlfwd.hh
noinst_HEADERS =
Modified: short/3D/PyLith/trunk/libsrc/utils/petscfwd.h
===================================================================
--- short/3D/PyLith/trunk/libsrc/utils/petscfwd.h 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/utils/petscfwd.h 2007-04-10 01:28:52 UTC (rev 6533)
@@ -11,13 +11,13 @@
//
/**
- * @file pylith/utils/petscfwd.hh
+ * @file pylith/utils/petscfwd.h
*
* @brief Forward declarations for Petsc objects.
*/
-#if !defined(pylith_feassemble_petscfwd_hh)
-#define pylith_feassemble_petscfwd_hh
+#if !defined(pylith_utils_petscfwd_h)
+#define pylith_utils_petscfwd_h
/// forward declaration for PETSc Mat
@@ -33,6 +33,6 @@
typedef int PetscErrorCode;
-#endif
+#endif // pylith_utils_petscfwd_h
// End of file
Added: short/3D/PyLith/trunk/libsrc/utils/stlfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/utils/stlfwd.hh 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/libsrc/utils/stlfwd.hh 2007-04-10 01:28:52 UTC (rev 6533)
@@ -0,0 +1,55 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/utils/stlfwd.hh
+ *
+ * @brief Forward declarations for C++ STL objects.
+ *
+ * For simple types (i.e., int and double) std::valarray provides some
+ * features that std::vector does not have, such as operating on the
+ * whole array at once.
+ */
+
+#if !defined(pylith_utils_stlfwd_hh)
+#define pylith_utils_stlfwd_hh
+
+/// Forward declaration of STL vector
+namespace std {
+ // std::vector
+ template<typename T> class allocator;
+ template<typename T, typename U> class vector;
+
+ // std::valarray
+ template<typename T> class valarray;
+} // std
+
+/// Aliases
+namespace pylith {
+ /// Alias for std::vector<int>
+ typedef std::vector<int, std::allocator<int> > int_vector;
+
+ /// Alias for std::vector<double>
+ typedef std::vector<double, std::allocator<double> > double_vector;
+
+ /// Alias for std::valarray<int>
+ typedef std::valarray<int> int_array;
+
+ /// Alias for std::valarray<double>
+ typedef std::valarray<double> double_array;
+
+} // pylith
+
+
+#endif // pylith_utils_stlfwd_hh
+
+// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2007-04-10 01:28:52 UTC (rev 6533)
@@ -20,6 +20,8 @@
#include "pylith/feassemble/ParameterManager.hh" // USES ParameterManager
+#include <valarray> // USES std::valarray (double_array)
+
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestElasticMaterial );
@@ -117,15 +119,14 @@
parameterLambda->updateAddPoint(*cellIter, cellData);
material.initCellData(numQuadPts);
- const double* density = material.calcDensity(*cellIter, numQuadPts);
+ const std::vector<double_array>& density = material.calcDensity(*cellIter);
const double tolerance = 1.0e-06;
const double* densityE = data.density;
- CPPUNIT_ASSERT(0 != density);
CPPUNIT_ASSERT(0 != densityE);
const double size = numQuadPts;
for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, density[i]/densityE[i], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, density[i][0]/densityE[i], tolerance);
} // testCalcProperties
// ----------------------------------------------------------------------
@@ -202,21 +203,45 @@
cellData[1] = data.parameterData[5];
parameterLambda->updateAddPoint(*cellIter, cellData);
+ int strainSize = 0;
+ switch(data.dimension)
+ { // switch
+ case 1 :
+ strainSize = 1;
+ break;
+ case 2 :
+ strainSize = 3;
+ break;
+ case 3 :
+ strainSize = 6;
+ break;
+ } // switch
+ std::vector<double_array> strain(numQuadPts);
+ for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
+ strain[iQuad].resize(strainSize);
+ for (int iStrain=0; iStrain < strainSize; ++iStrain, ++i)
+ strain[iQuad][iStrain] = data.strain[i];
+ } // for
+
material.initCellData(numQuadPts);
- const double* stress = material.calcStress(*cellIter, data.strain,
- numQuadPts, data.dimension);
+ const std::vector<double_array>& stress =
+ material.calcStress(*cellIter, strain);
const double tolerance = 1.0e-06;
const double* stressE = data.stress;
- CPPUNIT_ASSERT(0 != stress);
CPPUNIT_ASSERT(0 != stressE);
- const double size = numQuadPts * material.stressSize();
- for (int i=0; i < size; ++i)
- if (fabs(stressE[i]) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, stress[i]/stressE[i],
+
+ for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
+ const int stressSize = stress[iQuad].size();
+ for (int iStress=0; iStress < stressSize; ++iStress, ++i) {
+ if (fabs(stressE[i]) > tolerance)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, stress[iQuad][iStress]/stressE[i],
+ tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(stressE[i], stress[iQuad][iStress],
tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(stressE[i], stress[i], tolerance);
+ } // for
+ } // for
} // testCalcStress
// ----------------------------------------------------------------------
@@ -293,23 +318,47 @@
cellData[1] = data.parameterData[5];
parameterLambda->updateAddPoint(*cellIter, cellData);
+ int strainSize = 0;
+ switch(data.dimension)
+ { // switch
+ case 1 :
+ strainSize = 1;
+ break;
+ case 2 :
+ strainSize = 3;
+ break;
+ case 3 :
+ strainSize = 6;
+ break;
+ } // switch
+ std::vector<double_array> strain(numQuadPts);
+ for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
+ strain[iQuad].resize(strainSize);
+ for (int iStrain=0; iStrain < strainSize; ++iStrain, ++i)
+ strain[iQuad][iStrain] = data.strain[i];
+ } // for
+
material.initCellData(numQuadPts);
- const double* elasticConsts =
- material.calcDerivElastic(*cellIter, data.strain,
- numQuadPts, data.dimension);
+ const std::vector<double_array>& elasticConsts =
+ material.calcDerivElastic(*cellIter, strain);
const double tolerance = 1.0e-06;
const double* elasticConstsE = data.elasticConsts;
- CPPUNIT_ASSERT(0 != elasticConsts);
CPPUNIT_ASSERT(0 != elasticConstsE);
- const double size = numQuadPts * material.numElasticConsts();
- for (int i=0; i < size; ++i)
+
+ for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
+ const int numElasticConsts = elasticConsts[iQuad].size();
+ for (int iConst=0; iConst < numElasticConsts; ++iConst, ++i) {
if (fabs(elasticConstsE[i]) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, elasticConsts[i]/elasticConstsE[i],
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,
+ elasticConsts[iQuad][iConst]/elasticConstsE[i],
tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(elasticConstsE[i], elasticConsts[i],
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(elasticConstsE[i],
+ elasticConsts[iQuad][iConst],
tolerance);
+ } // for
+ } // for
} // testCalcDerivElastic
// ----------------------------------------------------------------------
@@ -319,13 +368,13 @@
{ // testInitCellData
ElasticIsotropic3D material;
- CPPUNIT_ASSERT(0 == material._density);
- CPPUNIT_ASSERT(0 == material._elasticConsts);
+ CPPUNIT_ASSERT(0 == material._density.size());
+ CPPUNIT_ASSERT(0 == material._elasticConsts.size());
const int numQuadPts = 4;
material.initCellData(numQuadPts);
- CPPUNIT_ASSERT(0 != material._density);
- CPPUNIT_ASSERT(0 != material._stress);
- CPPUNIT_ASSERT(0 != material._elasticConsts);
+ CPPUNIT_ASSERT(numQuadPts == material._density.size());
+ CPPUNIT_ASSERT(numQuadPts == material._stress.size());
+ CPPUNIT_ASSERT(numQuadPts == material._elasticConsts.size());
} // testInitCellData
// ----------------------------------------------------------------------
@@ -335,14 +384,19 @@
ElasticMaterial* material,
const ElasticMaterialData& data) const
{ // _testCalcDensity
- double density = 0;
- material->_calcDensity(&density, 1, data.parameterData, data.numParameters);
+ const int numParameters = data.numParameters;
+ double_array parameters(numParameters);
+ for (int i=0; i < numParameters; ++i)
+ parameters[i] = data.parameterData[i];
+
+ double_array density(1);
+ material->_calcDensity(&density, parameters);
const double* densityE = data.density;
CPPUNIT_ASSERT(0 != densityE);
const double tolerance = 1.0e-06;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, density/densityE[0], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, density[0]/densityE[0], tolerance);
} // _testCalcDensity
// ----------------------------------------------------------------------
@@ -352,25 +406,41 @@
ElasticMaterial* material,
const ElasticMaterialData& data) const
{ // _testCalcElasticConsts
- const int size = material->stressSize();
- double* stress = (size > 0) ? new double[size] : 0;
- material->_calcStress(stress, size,
- data.parameterData, data.numParameters,
- data.strain, data.dimension);
+ const int numParameters = data.numParameters;
+ double_array parameters(numParameters);
+ for (int i=0; i < numParameters; ++i)
+ parameters[i] = data.parameterData[i];
+
+ int stressSize = 0;
+ switch(data.dimension)
+ { // switch
+ case 1 :
+ stressSize = 1;
+ break;
+ case 2 :
+ stressSize = 3;
+ break;
+ case 3 :
+ stressSize = 6;
+ break;
+ } // switch
+ double_array stress(stressSize);
+ double_array strain(stressSize);
+ for (int i=0; i < stressSize; ++i)
+ strain[i] = data.strain[i];
+ material->_calcStress(&stress, parameters, strain);
const double* stressE = data.stress;
- CPPUNIT_ASSERT(0 != stress);
CPPUNIT_ASSERT(0 != stressE);
const double tolerance = 1.0e-06;
- for (int i=0; i < size; ++i)
+ for (int i=0; i < stressSize; ++i)
if (fabs(stressE[i]) > tolerance)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, stress[i]/stressE[i],
tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(stressE[i], stress[i],
tolerance);
- delete[] stress; stress = 0;
} // _testCalcStress
// ----------------------------------------------------------------------
@@ -380,25 +450,45 @@
ElasticMaterial* material,
const ElasticMaterialData& data) const
{ // _testCalcElasticConsts
- const int size = material->numElasticConsts();
- double* elasticConsts = (size > 0) ? new double[size] : 0;
- material->_calcElasticConsts(elasticConsts, size,
- data.parameterData, data.numParameters,
- data.strain, data.dimension);
+ const int numParameters = data.numParameters;
+ double_array parameters(numParameters);
+ for (int i=0; i < numParameters; ++i)
+ parameters[i] = data.parameterData[i];
+
+ int numConsts = 0;
+ int strainSize = 0;
+ switch(data.dimension)
+ { // switch
+ case 1 :
+ numConsts = 1;
+ strainSize = 1;
+ break;
+ case 2 :
+ numConsts = 6;
+ strainSize = 3;
+ break;
+ case 3 :
+ numConsts = 21;
+ strainSize = 6;
+ break;
+ } // switch
+ double_array elasticConsts(numConsts);
+ double_array strain(strainSize);
+ for (int i=0; i < strainSize; ++i)
+ strain[i] = data.strain[i];
+ material->_calcElasticConsts(&elasticConsts, parameters, strain);
const double* elasticConstsE = data.elasticConsts;
- CPPUNIT_ASSERT(0 != elasticConsts);
CPPUNIT_ASSERT(0 != elasticConstsE);
const double tolerance = 1.0e-06;
- for (int i=0; i < size; ++i)
+ for (int i=0; i < numConsts; ++i)
if (fabs(elasticConstsE[i]) > tolerance)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, elasticConsts[i]/elasticConstsE[i],
tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(elasticConstsE[i], elasticConsts[i],
tolerance);
- delete[] elasticConsts; elasticConsts = 0;
} // _testCalcElasticConsts
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc 2007-04-09 20:45:03 UTC (rev 6532)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc 2007-04-10 01:28:52 UTC (rev 6533)
@@ -25,6 +25,7 @@
#include <petscmesh.h> // USES PETSc Mesh
+#include <valarray> // USES std::valarray (double_array)
#include <math.h> // USES assert()
// ----------------------------------------------------------------------
@@ -197,22 +198,20 @@
{ // _testDBToParameters
CPPUNIT_ASSERT(0 != material);
- const int numDBValues = data.numDBValues;
- double* const dbData = data.dbData;
+ double_array dbData(data.numDBValues);
+ for (int i=0; i < data.numDBValues; ++i)
+ dbData[i] = data.dbData[i];
+
const int numParameters = data.numParameters;
double* const parameterDataE = data.parameterData;
- double* parameterData = (numParameters > 0) ? new double[numParameters] : 0;
- material->_dbToParameters(parameterData, numParameters,
- dbData, numDBValues);
+ double_array parameterData(numParameters);
+ material->_dbToParameters(¶meterData, dbData);
const double tolerance = 1.0e-06;
for (int i=0; i < numParameters; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,
- parameterData[i]/parameterDataE[i],
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, parameterData[i]/parameterDataE[i],
tolerance);
-
- delete[] parameterData; parameterData = 0;
} // _testDBToParameters
// ----------------------------------------------------------------------
More information about the cig-commits
mailing list