[cig-commits] r8984 - in short/3D/PyLith/trunk:
examples/twocells/twotet4-geoproj libsrc/feassemble
libsrc/materials libsrc/meshio unittests/libtests/feassemble
unittests/libtests/materials unittests/libtests/materials/data
brad at geodynamics.org
brad at geodynamics.org
Thu Jan 3 15:00:11 PST 2008
Author: brad
Date: 2008-01-03 15:00:07 -0800 (Thu, 03 Jan 2008)
New Revision: 8984
Modified:
short/3D/PyLith/trunk/examples/twocells/twotet4-geoproj/dislocation_slip.spatialdb
short/3D/PyLith/trunk/examples/twocells/twotet4-geoproj/pylithapp.cfg
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
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/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc
short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticIsotropic3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStrain.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStress.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStrain1D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStress1D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh
Log:
Optimized managing physical properties for materials. Switched from using a section for each parameter to a single section for all parameters.
Modified: short/3D/PyLith/trunk/examples/twocells/twotet4-geoproj/dislocation_slip.spatialdb
===================================================================
--- short/3D/PyLith/trunk/examples/twocells/twotet4-geoproj/dislocation_slip.spatialdb 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/examples/twocells/twotet4-geoproj/dislocation_slip.spatialdb 2008-01-03 23:00:07 UTC (rev 8984)
@@ -14,11 +14,11 @@
value-units = m m m
// The values are specified at one spatial location.
- num-locs = 1
+ num-locs = 3
// The dimension of the spatial distribution is 2, since the same data
// is specified on a plane.
- data-dim = 2
+ data-dim = 1
// The spatial dimension of the database is 3.
space-dim = 3
Modified: short/3D/PyLith/trunk/examples/twocells/twotet4-geoproj/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/twocells/twotet4-geoproj/pylithapp.cfg 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/examples/twocells/twotet4-geoproj/pylithapp.cfg 2008-01-03 23:00:07 UTC (rev 8984)
@@ -86,7 +86,7 @@
# The properties for this material are obtained from the SCEC CVM-H.
db = spatialdata.spatialdb.SCECCVMH
-db.data_dir = /Users/willic3/geoframe/tools/vx52/bin
+db.data_dir = /Users/brad/data/sceccvm-h/vx52/bin
# We are doing 3D quadrature for a tetrahedron.
quadrature = pylith.feassemble.quadrature.Quadrature3D
@@ -108,5 +108,5 @@
log_summary = true
ksp_max_it = 2000
ksp_gmres_restart = 100
-start_in_debugger = true
+#start_in_debugger = true
debugger_timeout = 100
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -81,7 +81,7 @@
{ // integrateResidual
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityExplicit::*elasticityResidual_fn_type)
- (const std::vector<double_array>&);
+ (const double_array&);
assert(0 != _quadrature);
assert(0 != _material);
@@ -165,11 +165,8 @@
double_array dispTmdtCell(cellVecSize);
// Allocate vector for total strain
- std::vector<double_array> totalStrain(numQuadPts);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- totalStrain[iQuad].resize(tensorSize);
- totalStrain[iQuad] = 0.0;
- } // for
+ double_array totalStrain(numQuadPts*tensorSize);
+ totalStrain = 0.0;
#ifdef FASTER
fields->createCustomAtlas("material-id", materialId);
@@ -210,11 +207,10 @@
const double_array& jacobianDet = _quadrature->jacobianDet();
// Compute action for inertial terms
- const std::vector<double_array>& density =
- _material->calcDensity();
+ const double_array& density = _material->calcDensity();
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt =
- quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad][0] / dt2;
+ quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
for (int iBasis=0; iBasis < numBasis; ++iBasis) {
const double valI = wt*basis[iQuad*numBasis+iBasis];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
@@ -229,9 +225,9 @@
PetscLogFlopsNoCheck(numQuadPts*(3+numBasis*(1+numBasis*(6*spaceDim))));
// Compute B(transpose) * sigma, first computing strains
- calcTotalStrainFn(&totalStrain, basisDeriv, dispTCell, numBasis);
- const std::vector<double_array>& stress =
- _material->calcStress(totalStrain);
+ calcTotalStrainFn(&totalStrain, basisDeriv, dispTCell,
+ numBasis, numQuadPts);
+ const double_array& stress = _material->calcStress(totalStrain);
CALL_MEMBER_FN(*this, elasticityResidualFn)(stress);
// Assemble cell contribution into field
@@ -307,12 +303,12 @@
const double_array& jacobianDet = _quadrature->jacobianDet();
// Get material physical properties at quadrature points for this cell
- const std::vector<double_array>& density = _material->calcDensity();
+ const double_array& density = _material->calcDensity();
// Compute Jacobian for inertial terms
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt =
- quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad][0] / dt2;
+ quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / 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/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -80,7 +80,7 @@
{ // integrateResidual
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityImplicit::*elasticityResidual_fn_type)
- (const std::vector<double_array>&);
+ (const double_array&);
assert(0 != _quadrature);
assert(0 != _material);
@@ -174,11 +174,8 @@
//double_array gravCell(cellVecSize);
// Allocate vector for total strain
- std::vector<double_array> totalStrain(numQuadPts);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- totalStrain[iQuad].resize(tensorSize);
- totalStrain[iQuad] = 0.0;
- } // for
+ double_array totalStrain(numQuadPts*tensorSize);
+ totalStrain = 0.0;
PetscLogEventEnd(setupEvent,0,0,0,0);
// Loop over cells
@@ -222,7 +219,7 @@
// with gravity vector.
// Get density at quadrature points for this cell
- const std::vector<double_array>& density = _material->calcDensity();
+ const double_array& density = _material->calcDensity();
// Compute action for element body forces
if (!grav.isNull()) {
@@ -243,9 +240,9 @@
// Compute B(transpose) * sigma, first computing strains
PetscLogEventBegin(stressEvent,0,0,0,0);
- calcTotalStrainFn(&totalStrain, basisDeriv, dispTBctpdtCell, numBasis);
- const std::vector<double_array>& stress =
- _material->calcStress(totalStrain);
+ calcTotalStrainFn(&totalStrain, basisDeriv, dispTBctpdtCell,
+ numBasis, numQuadPts);
+ const double_array& stress = _material->calcStress(totalStrain);
PetscLogEventEnd(stressEvent,0,0,0,0);
PetscLogEventBegin(computeEvent,0,0,0,0);
@@ -281,7 +278,7 @@
{ // integrateJacobian
/// Member prototype for _elasticityJacobianXD()
typedef void (pylith::feassemble::ElasticityImplicit::*elasticityJacobian_fn_type)
- (const std::vector<double_array>&);
+ (const double_array&);
assert(0 != _quadrature);
assert(0 != _material);
@@ -354,11 +351,8 @@
double_array dispTBctpdtCell(cellVecSize);
// Allocate vector for total strain
- std::vector<double_array> totalStrain(numQuadPts);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- totalStrain[iQuad].resize(tensorSize);
- totalStrain[iQuad] = 0.0;
- } // for
+ double_array totalStrain(numQuadPts*tensorSize);
+ totalStrain = 0.0;
#ifdef FASTER
fields->createCustomAtlas("material-id", materialId);
@@ -393,10 +387,11 @@
const double_array& jacobianDet = _quadrature->jacobianDet();
// Compute strains
- calcTotalStrainFn(&totalStrain, basisDeriv, dispTBctpdtCell, numBasis);
+ calcTotalStrainFn(&totalStrain, basisDeriv, dispTBctpdtCell,
+ numBasis, numQuadPts);
// Get "elasticity" matrix at quadrature points for this cell
- const std::vector<double_array>& elasticConsts =
+ const double_array& elasticConsts =
_material->calcDerivElastic(totalStrain);
CALL_MEMBER_FN(*this, elasticityJacobianFn)(elasticConsts);
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -117,11 +117,8 @@
double_array dispCell(cellVecSize);
// Allocate vector for total strain
- std::vector<double_array> totalStrain(numQuadPts);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- totalStrain[iQuad].resize(tensorSize);
- totalStrain[iQuad] = 0.0;
- } // for
+ double_array totalStrain(numQuadPts*tensorSize);
+ totalStrain = 0.0;
const ALE::Obj<real_section_type>& disp = fields->getSolution();
#ifdef FASTER
@@ -147,7 +144,7 @@
const double_array& basisDeriv = _quadrature->basisDeriv();
// Compute strains
- calcTotalStrainFn(&totalStrain, basisDeriv, dispCell, numBasis);
+ calcTotalStrainFn(&totalStrain, basisDeriv, dispCell, numBasis, numQuadPts);
// Update material state
_material->updateState(totalStrain, *c_iter);
@@ -215,7 +212,7 @@
// Integrate elasticity term in residual for 1-D cells.
void
pylith::feassemble::IntegratorElasticity::_elasticityResidual1D(
- const std::vector<double_array>& stress)
+ const double_array& stress)
{ // _elasticityResidual1D
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
@@ -230,7 +227,7 @@
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
+ const double s11 = stress[iQuad];
for (int iBasis=0; iBasis < numBasis; ++iBasis) {
const double N1 = wt*basisDeriv[iQuad*numBasis+iBasis ];
_cellVector[iBasis*spaceDim ] -= N1*s11;
@@ -243,7 +240,7 @@
// Integrate elasticity term in residual for 2-D cells.
void
pylith::feassemble::IntegratorElasticity::_elasticityResidual2D(
- const std::vector<double_array>& stress)
+ const double_array& stress)
{ // _elasticityResidual2D
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
@@ -255,12 +252,13 @@
assert(2 == cellDim);
assert(quadWts.size() == numQuadPts);
-
+ const int stressSize = 3;
+
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- const double s22 = stress[iQuad][1];
- const double s12 = stress[iQuad][2];
+ const double s11 = stress[iQuad*stressSize+0];
+ const double s22 = stress[iQuad*stressSize+1];
+ const double s12 = stress[iQuad*stressSize+2];
for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
iBasis < numBasis;
++iBasis) {
@@ -277,7 +275,7 @@
// Integrate elasticity term in residual for 3-D cells.
void
pylith::feassemble::IntegratorElasticity::_elasticityResidual3D(
- const std::vector<double_array>& stress)
+ const double_array& stress)
{ // _elasticityResidual3D
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
@@ -289,15 +287,16 @@
assert(3 == cellDim);
assert(quadWts.size() == numQuadPts);
+ const int stressSize = 6;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double s11 = stress[iQuad][0];
- const double s22 = stress[iQuad][1];
- const double s33 = stress[iQuad][2];
- const double s12 = stress[iQuad][3];
- const double s23 = stress[iQuad][4];
- const double s13 = stress[iQuad][5];
+ const double s11 = stress[iQuad*stressSize+0];
+ const double s22 = stress[iQuad*stressSize+1];
+ const double s33 = stress[iQuad*stressSize+2];
+ const double s12 = stress[iQuad*stressSize+3];
+ const double s23 = stress[iQuad*stressSize+4];
+ const double s13 = stress[iQuad*stressSize+5];
for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
iBasis < numBasis;
@@ -319,7 +318,7 @@
// Integrate elasticity term in Jacobian for 1-D cells.
void
pylith::feassemble::IntegratorElasticity::_elasticityJacobian1D(
- const std::vector<double_array>& elasticConsts)
+ const double_array& elasticConsts)
{ // _elasticityJacobian1D
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
@@ -334,7 +333,7 @@
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- const double C1111 = elasticConsts[iQuad][0];
+ const double C1111 = elasticConsts[iQuad];
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const double valI = wt*basisDeriv[iQ+iBasis]*C1111;
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
@@ -352,7 +351,7 @@
// Integrate elasticity term in Jacobian for 2-D cells.
void
pylith::feassemble::IntegratorElasticity::_elasticityJacobian2D(
- const std::vector<double_array>& elasticConsts)
+ const double_array& elasticConsts)
{ // _elasticityJacobian2D
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
@@ -364,19 +363,20 @@
assert(2 == cellDim);
assert(quadWts.size() == numQuadPts);
-
+ const int numConsts = 6;
+
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
// tau_ij = C_ijkl * e_kl
// = C_ijlk * 0.5 (u_k,l + u_l,k)
// = 0.5 * C_ijkl * (u_k,l + u_l,k)
// divide C_ijkl by 2 if k != l
- const double C1111 = elasticConsts[iQuad][0];
- const double C1122 = elasticConsts[iQuad][1];
- const double C1112 = elasticConsts[iQuad][2]/2.0;
- const double C2222 = elasticConsts[iQuad][3];
- const double C2212 = elasticConsts[iQuad][4]/2.0;
- const double C1212 = elasticConsts[iQuad][5]/2.0;
+ const double C1111 = elasticConsts[iQuad*numConsts+0];
+ const double C1122 = elasticConsts[iQuad*numConsts+1];
+ const double C1112 = elasticConsts[iQuad*numConsts+2]/2.0;
+ const double C2222 = elasticConsts[iQuad*numConsts+3];
+ const double C2212 = elasticConsts[iQuad*numConsts+4]/2.0;
+ const double C1212 = elasticConsts[iQuad*numConsts+5]/2.0;
for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
iBasis < numBasis;
++iBasis) {
@@ -415,7 +415,7 @@
// Integrate elasticity term in Jacobian for 3-D cells.
void
pylith::feassemble::IntegratorElasticity::_elasticityJacobian3D(
- const std::vector<double_array>& elasticConsts)
+ const double_array& elasticConsts)
{ // _elasticityJacobian3D
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
@@ -427,7 +427,8 @@
assert(3 == cellDim);
assert(quadWts.size() == numQuadPts);
-
+ const int numConsts = 21;
+
// Compute Jacobian for consistent tangent matrix
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
@@ -435,27 +436,27 @@
// = C_ijlk * 0.5 (u_k,l + u_l,k)
// = 0.5 * C_ijkl * (u_k,l + u_l,k)
// divide C_ijkl by 2 if k != l
- const double C1111 = elasticConsts[iQuad][ 0];
- const double C1122 = elasticConsts[iQuad][ 1];
- const double C1133 = elasticConsts[iQuad][ 2];
- const double C1112 = elasticConsts[iQuad][ 3]/2.0;
- const double C1123 = elasticConsts[iQuad][ 4]/2.0;
- const double C1113 = elasticConsts[iQuad][ 5]/2.0;
- const double C2222 = elasticConsts[iQuad][ 6];
- const double C2233 = elasticConsts[iQuad][ 7];
- const double C2212 = elasticConsts[iQuad][ 8]/2.0;
- const double C2223 = elasticConsts[iQuad][ 9]/2.0;
- const double C2213 = elasticConsts[iQuad][10]/2.0;
- const double C3333 = elasticConsts[iQuad][11];
- const double C3312 = elasticConsts[iQuad][12]/2.0;
- const double C3323 = elasticConsts[iQuad][13]/2.0;
- const double C3313 = elasticConsts[iQuad][14]/2.0;
- const double C1212 = elasticConsts[iQuad][15]/2.0;
- const double C1223 = elasticConsts[iQuad][16]/2.0;
- const double C1213 = elasticConsts[iQuad][17]/2.0;
- const double C2323 = elasticConsts[iQuad][18]/2.0;
- const double C2313 = elasticConsts[iQuad][19]/2.0;
- const double C1313 = elasticConsts[iQuad][20]/2.0;
+ const double C1111 = elasticConsts[iQuad*numConsts+ 0];
+ const double C1122 = elasticConsts[iQuad*numConsts+ 1];
+ const double C1133 = elasticConsts[iQuad*numConsts+ 2];
+ const double C1112 = elasticConsts[iQuad*numConsts+ 3] / 2.0;
+ const double C1123 = elasticConsts[iQuad*numConsts+ 4] / 2.0;
+ const double C1113 = elasticConsts[iQuad*numConsts+ 5] / 2.0;
+ const double C2222 = elasticConsts[iQuad*numConsts+ 6];
+ const double C2233 = elasticConsts[iQuad*numConsts+ 7];
+ const double C2212 = elasticConsts[iQuad*numConsts+ 8] / 2.0;
+ const double C2223 = elasticConsts[iQuad*numConsts+ 9] / 2.0;
+ const double C2213 = elasticConsts[iQuad*numConsts+10] / 2.0;
+ const double C3333 = elasticConsts[iQuad*numConsts+11];
+ const double C3312 = elasticConsts[iQuad*numConsts+12] / 2.0;
+ const double C3323 = elasticConsts[iQuad*numConsts+13] / 2.0;
+ const double C3313 = elasticConsts[iQuad*numConsts+14] / 2.0;
+ const double C1212 = elasticConsts[iQuad*numConsts+15] / 2.0;
+ const double C1223 = elasticConsts[iQuad*numConsts+16] / 2.0;
+ const double C1213 = elasticConsts[iQuad*numConsts+17] / 2.0;
+ const double C2323 = elasticConsts[iQuad*numConsts+18] / 2.0;
+ const double C2313 = elasticConsts[iQuad*numConsts+19] / 2.0;
+ const double C1313 = elasticConsts[iQuad*numConsts+20] / 2.0;
for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
iBasis < numBasis;
++iBasis) {
@@ -526,91 +527,92 @@
// ----------------------------------------------------------------------
void
pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D(
- std::vector<double_array>* strain,
+ double_array* strain,
const double_array& basisDeriv,
const double_array& disp,
- const int numBasis)
+ const int numBasis,
+ const int numQuadPts)
{ // calcTotalStrain1D
assert(0 != strain);
+
+ const int dim = 1;
- const int dim = 1;
- const int numQuadPts = strain->size();
-
assert(basisDeriv.size() == numQuadPts*numBasis*dim);
assert(disp.size() == numBasis*dim);
+ (*strain) = 0.0;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- assert(1 == (*strain)[iQuad].size());
- (*strain)[iQuad] *= 0.0;
for (int iBasis=0; iBasis < numBasis; ++iBasis)
- (*strain)[iQuad][0] +=
- basisDeriv[iQuad*numBasis+iBasis] * disp[iBasis];
+ (*strain)[iQuad] += basisDeriv[iQuad*numBasis+iBasis] * disp[iBasis];
} // for
} // calcTotalStrain1D
// ----------------------------------------------------------------------
void
pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D(
- std::vector<double_array>* strain,
+ double_array* strain,
const double_array& basisDeriv,
const double_array& disp,
- const int numBasis)
+ const int numBasis,
+ const int numQuadPts)
{ // calcTotalStrain2D
assert(0 != strain);
const int dim = 2;
- const int numQuadPts = strain->size();
assert(basisDeriv.size() == numQuadPts*numBasis*dim);
assert(disp.size() == numBasis*dim);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- assert(3 == (*strain)[iQuad].size());
- (*strain)[iQuad] *= 0.0;
+ (*strain) = 0.0;
+ const int strainSize = 3;
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
- (*strain)[iQuad][0] += basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
- (*strain)[iQuad][1] += basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
- (*strain)[iQuad][2] +=
+ (*strain)[iQuad*strainSize+0] +=
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
+ (*strain)[iQuad*strainSize+1] +=
+ basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
+ (*strain)[iQuad*strainSize+2] +=
0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
} // for
- } // for
} // calcTotalStrain2D
// ----------------------------------------------------------------------
void
pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D(
- std::vector<double_array>* strain,
+ double_array* strain,
const double_array& basisDeriv,
const double_array& disp,
- const int numBasis)
+ const int numBasis,
+ const int numQuadPts)
{ // calcTotalStrain3D
assert(0 != strain);
const int dim = 3;
- const int numQuadPts = strain->size();
assert(basisDeriv.size() == numQuadPts*numBasis*dim);
assert(disp.size() == numBasis*dim);
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- assert(6 == (*strain)[iQuad].size());
- (*strain)[iQuad] *= 0.0;
+ (*strain) = 0.0;
+ const int strainSize = 6;
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
for (int iBasis=0, iQ=iQuad*numBasis*dim; iBasis < numBasis; ++iBasis) {
- (*strain)[iQuad][0] += basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
- (*strain)[iQuad][1] += basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
- (*strain)[iQuad][2] += basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+2];
- (*strain)[iQuad][3] +=
+ (*strain)[iQuad*strainSize+0] +=
+ basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim ];
+ (*strain)[iQuad*strainSize+1] +=
+ basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+1];
+ (*strain)[iQuad*strainSize+2] +=
+ basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+2];
+ (*strain)[iQuad*strainSize+3] +=
0.5 * (basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim ] +
basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+1]);
- (*strain)[iQuad][4] +=
+ (*strain)[iQuad*strainSize+4] +=
0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim+1] +
basisDeriv[iQ+iBasis*dim+1] * disp[iBasis*dim+2]);
- (*strain)[iQuad][5] +=
+ (*strain)[iQuad*strainSize+5] +=
0.5 * (basisDeriv[iQ+iBasis*dim+2] * disp[iBasis*dim ] +
basisDeriv[iQ+iBasis*dim ] * disp[iBasis*dim+2]);
} // for
- } // for
} // calcTotalStrain3D
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -41,9 +41,10 @@
// PUBLIC TYPEDEFS //////////////////////////////////////////////////////
public :
- typedef void (*totalStrain_fn_type)(std::vector<double_array>*,
+ typedef void (*totalStrain_fn_type)(double_array*,
const double_array&,
const double_array&,
+ const int,
const int);
@@ -91,37 +92,37 @@
*
* @param stress Stress tensor for cell at quadrature points.
*/
- void _elasticityResidual1D(const std::vector<double_array>& stress);
+ void _elasticityResidual1D(const double_array& stress);
/** Integrate elasticity term in residual for 2-D cells.
*
* @param stress Stress tensor for cell at quadrature points.
*/
- void _elasticityResidual2D(const std::vector<double_array>& stress);
+ void _elasticityResidual2D(const double_array& stress);
/** Integrate elasticity term in residual for 3-D cells.
*
* @param stress Stress tensor for cell at quadrature points.
*/
- void _elasticityResidual3D(const std::vector<double_array>& stress);
+ void _elasticityResidual3D(const double_array& stress);
/** Integrate elasticity term in Jacobian for 1-D cells.
*
* @param elasticConsts Matrix of elasticity constants at quadrature points.
*/
- void _elasticityJacobian1D(const std::vector<double_array>& elasticConsts);
+ void _elasticityJacobian1D(const double_array& elasticConsts);
/** Integrate elasticity term in Jacobian for 2-D cells.
*
* @param elasticConsts Matrix of elasticity constants at quadrature points.
*/
- void _elasticityJacobian2D(const std::vector<double_array>& elasticConsts);
+ void _elasticityJacobian2D(const double_array& elasticConsts);
/** Integrate elasticity term in Jacobian for 3-D cells.
*
* @param elasticConsts Matrix of elasticity constants at quadrature points.
*/
- void _elasticityJacobian3D(const std::vector<double_array>& elasticConsts);
+ void _elasticityJacobian3D(const double_array& elasticConsts);
/** Compute total strain in at quadrature points of a cell.
*
@@ -130,13 +131,14 @@
* @param disp Displacement at vertices of cell.
* @param dimension Dimension of cell.
* @param numBasis Number of basis functions for cell.
+ * @param numQuadPts Number of quadrature points.
*/
-
static
- void _calcTotalStrain1D(std::vector<double_array>* strain,
+ void _calcTotalStrain1D(double_array* strain,
const double_array& basisDeriv,
const double_array& disp,
- const int numBasis);
+ const int numBasis,
+ const int numQuadPts);
/** Compute total strain in at quadrature points of a cell.
*
@@ -144,13 +146,14 @@
* @param basisDeriv Derivatives of basis functions at quadrature points.
* @param disp Displacement at vertices of cell.
* @param numBasis Number of basis functions for cell.
+ * @param numQuadPts Number of quadrature points.
*/
-
static
- void _calcTotalStrain2D(std::vector<double_array>* strain,
+ void _calcTotalStrain2D(double_array* strain,
const double_array& basisDeriv,
const double_array& disp,
- const int numBasis);
+ const int numBasis,
+ const int numQuadPts);
/** Compute total strain in at quadrature points of a cell.
*
@@ -158,13 +161,14 @@
* @param basisDeriv Derivatives of basis functions at quadrature points.
* @param disp Displacement at vertices of cell.
* @param numBasis Number of basis functions for cell.
+ * @param numQuadPts Number of quadrature points.
*/
-
static
- void _calcTotalStrain3D(std::vector<double_array>* strain,
+ void _calcTotalStrain3D(double_array* strain,
const double_array& basisDeriv,
const double_array& disp,
- const int numBasis);
+ const int numBasis,
+ const int numQuadPts);
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -45,12 +45,11 @@
/// Parameters
const int numParameters = 3;
const int numParamValues[] = { 1, 1, 1 };
- const char* namesParameters[] = { "density", "mu", "lambda" };
/// Indices (order) of parameters
const int pidDensity = 0;
- const int pidMu = 1;
- const int pidLambda = 2;
+ const int pidMu = pidDensity + 1;
+ const int pidLambda = pidMu + 1;
} // _ElasticIsotropic3D
} // materials
@@ -89,26 +88,17 @@
} // _numDBValues
// ----------------------------------------------------------------------
-// Get names of parameters for physical properties.
-const char**
-pylith::materials::ElasticIsotropic3D::_parameterNames(void) const
-{ // _parameterNames
- return _ElasticIsotropic3D::namesParameters;
-} // _parameterNames
-
-// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
-pylith::materials::ElasticIsotropic3D::_dbToParameters(std::vector<double_array>* paramVals,
- const double_array& dbValues) const
+pylith::materials::ElasticIsotropic3D::_dbToParameters(
+ double* const paramVals,
+ const int numParams,
+ const double_array& dbValues) const
{ // _dbToParameters
assert(0 != paramVals);
- const int numParams = paramVals->size();
- assert(_ElasticIsotropic3D::numParameters == numParams);
+ assert(_numParamsQuadPt == numParams);
const int numDBValues = dbValues.size();
assert(_ElasticIsotropic3D::numDBValues == numDBValues);
- for (int i=0; i < numParams; ++i)
- assert(_ElasticIsotropic3D::numParamValues[i] == (*paramVals)[i].size());
const double density = dbValues[_ElasticIsotropic3D::didDensity];
const double vs = dbValues[_ElasticIsotropic3D::didVs];
@@ -125,9 +115,9 @@
throw std::runtime_error(msg.str());
} // if
- (*paramVals)[_ElasticIsotropic3D::pidDensity][0] = density;
- (*paramVals)[_ElasticIsotropic3D::pidMu][0] = mu;
- (*paramVals)[_ElasticIsotropic3D::pidLambda][0] = lambda;
+ paramVals[_ElasticIsotropic3D::pidDensity] = density;
+ paramVals[_ElasticIsotropic3D::pidMu] = mu;
+ paramVals[_ElasticIsotropic3D::pidLambda] = lambda;
PetscLogFlopsNoCheck(6);
} // _dbToParameters
@@ -152,36 +142,38 @@
// Compute density at location from parameters.
void
pylith::materials::ElasticIsotropic3D::_calcDensity(
- double_array* const density,
- const std::vector<double_array>& parameters)
+ double* const density,
+ const double* parameters,
+ const int numParams)
{ // _calcDensity
assert(0 != density);
- assert(1 == density->size());
- assert(_ElasticIsotropic3D::numParameters == parameters.size());
- assert(1 == parameters[_ElasticIsotropic3D::pidDensity].size());
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
- (*density)[0] = parameters[_ElasticIsotropic3D::pidDensity][0];
+ density[0] = parameters[_ElasticIsotropic3D::pidDensity];
} // _calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor at location from parameters.
void
pylith::materials::ElasticIsotropic3D::_calcStress(
- double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+ double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _calcStress
assert(0 != stress);
- assert(_ElasticIsotropic3D::tensorSize == stress->size());
- assert(_ElasticIsotropic3D::numParameters == parameters.size());
- assert(_ElasticIsotropic3D::tensorSize == totalStrain.size());
- assert(1 == parameters[_ElasticIsotropic3D::pidDensity].size());
- assert(1 == parameters[_ElasticIsotropic3D::pidMu].size());
- assert(1 == parameters[_ElasticIsotropic3D::pidLambda].size());
+ assert(_ElasticIsotropic3D::tensorSize == stressSize);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_ElasticIsotropic3D::tensorSize == strainSize);
- const double density = parameters[_ElasticIsotropic3D::pidDensity][0];
- const double mu = parameters[_ElasticIsotropic3D::pidMu][0];
- const double lambda = parameters[_ElasticIsotropic3D::pidLambda][0];
+ const double density = parameters[_ElasticIsotropic3D::pidDensity];
+ const double mu = parameters[_ElasticIsotropic3D::pidMu];
+ const double lambda = parameters[_ElasticIsotropic3D::pidLambda];
const double mu2 = 2.0*mu;
@@ -194,12 +186,13 @@
const double s123 = lambda * (e11 + e22 + e33);
- (*stress)[0] = s123 + mu2*e11;
- (*stress)[1] = s123 + mu2*e22;
- (*stress)[2] = s123 + mu2*e33;
- (*stress)[3] = mu2 * e12;
- (*stress)[4] = mu2 * e23;
- (*stress)[5] = mu2 * e13;
+ stress[0] = s123 + mu2*e11;
+ stress[1] = s123 + mu2*e22;
+ stress[2] = s123 + mu2*e33;
+ stress[3] = mu2 * e12;
+ stress[4] = mu2 * e23;
+ stress[5] = mu2 * e13;
+
PetscLogFlopsNoCheck(13);
} // _calcStress
@@ -207,46 +200,49 @@
// Compute derivative of elasticity matrix at location from parameters.
void
pylith::materials::ElasticIsotropic3D::_calcElasticConsts(
- double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _calcElasticConsts
assert(0 != elasticConsts);
- assert(_ElasticIsotropic3D::numElasticConsts == elasticConsts->size());
- assert(_ElasticIsotropic3D::numParameters == parameters.size());
- assert(_ElasticIsotropic3D::tensorSize == totalStrain.size());
- assert(1 == parameters[_ElasticIsotropic3D::pidDensity].size());
- assert(1 == parameters[_ElasticIsotropic3D::pidMu].size());
- assert(1 == parameters[_ElasticIsotropic3D::pidLambda].size());
+ assert(_ElasticIsotropic3D::numElasticConsts == numElasticConsts);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_ElasticIsotropic3D::tensorSize == strainSize);
- const double density = parameters[_ElasticIsotropic3D::pidDensity][0];
- const double mu = parameters[_ElasticIsotropic3D::pidMu][0];
- const double lambda = parameters[_ElasticIsotropic3D::pidLambda][0];
+ const double density = parameters[_ElasticIsotropic3D::pidDensity];
+ const double mu = parameters[_ElasticIsotropic3D::pidMu];
+ const double lambda = parameters[_ElasticIsotropic3D::pidLambda];
const double mu2 = 2.0 * mu;
const double lambda2mu = lambda + mu2;
- (*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] = mu2; // C1212
- (*elasticConsts)[16] = 0; // C1223
- (*elasticConsts)[17] = 0; // C1213
- (*elasticConsts)[18] = mu2; // C2323
- (*elasticConsts)[19] = 0; // C2313
- (*elasticConsts)[20] = mu2; // 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] = mu2; // C1212
+ elasticConsts[16] = 0; // C1223
+ elasticConsts[17] = 0; // C1213
+ elasticConsts[18] = mu2; // C2323
+ elasticConsts[19] = 0; // C2313
+ elasticConsts[20] = mu2; // C1313
+
PetscLogFlopsNoCheck(2);
} // _calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -66,21 +66,17 @@
*/
int _numDBValues(void) const;
- /** Get names of parameters for physical properties.
- *
- * @returns Names of parameters
- */
- const char** _parameterNames(void) const;
-
/** Compute parameters from values in spatial database.
*
* Order of values in arrays matches order used in dbValues() and
* parameterNames().
*
* @param paramVals Array of parameters
+ * @param numParams Number of parameters at quadrature point.
* @param dbValues Array of database values
*/
- void _dbToParameters(std::vector<double_array>* paramVals,
+ void _dbToParameters(double* const paramVals,
+ const int numParams,
const double_array& dbValues) const;
/** Get number of entries in stress/strain tensors.
@@ -105,31 +101,45 @@
/** Compute density from parameters.
*
- * @param density Array for density
- * @param parameters Parameters at location
+ * @param density Array for density.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
*/
- void _calcDensity(double_array* const density,
- const std::vector<double_array>& parameters);
+ void _calcDensity(double* const density,
+ const double* parameters,
+ const int numParams);
/** Compute stress tensor from parameters.
*
- * @param stress Array for stress tensor
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcStress(double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
/** Compute derivatives of elasticity matrix from parameters.
*
- * @param elasticConsts Array for elastic constants
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcElasticConsts(double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcElasticConsts(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -21,8 +21,6 @@
#include <assert.h> // USES assert()
-#define FASTER
-
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::ElasticMaterial::ElasticMaterial(const int* numParamValues,
@@ -40,50 +38,59 @@
// ----------------------------------------------------------------------
// Compute density for cell at quadrature points.
-const std::vector<pylith::double_array>&
+const pylith::double_array&
pylith::materials::ElasticMaterial::calcDensity(void)
{ // calcDensity
const int numQuadPts = _numQuadPts;
- assert(_paramsCell.size() == numQuadPts);
+ const int numParamsQuadPt = _numParamsQuadPt;
+ assert(_paramsCell.size() == numQuadPts*numParamsQuadPt);
assert(_density.size() == numQuadPts);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- _calcDensity(&_density[iQuad], _paramsCell[iQuad]);
+ _calcDensity(&_density[iQuad],
+ &_paramsCell[iQuad*numParamsQuadPt], numParamsQuadPt);
return _density;
} // calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor for cell at quadrature points.
-const std::vector<pylith::double_array>&
-pylith::materials::ElasticMaterial::calcStress(
- const std::vector<double_array>& totalStrain)
+const pylith::double_array&
+pylith::materials::ElasticMaterial::calcStress(const double_array& totalStrain)
{ // calcStress
const int numQuadPts = _numQuadPts;
- assert(_paramsCell.size() == numQuadPts);
- assert(totalStrain.size() == numQuadPts);
- assert(_stress.size() == numQuadPts);
+ const int numParamsQuadPt = _numParamsQuadPt;
+ const int tensorSize = _tensorSize();
+ assert(_paramsCell.size() == numQuadPts*numParamsQuadPt);
+ assert(totalStrain.size() == numQuadPts*tensorSize);
+ assert(_stress.size() == numQuadPts*tensorSize);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- _calcStress(&_stress[iQuad], _paramsCell[iQuad], totalStrain[iQuad]);
+ _calcStress(&_stress[iQuad*tensorSize], tensorSize,
+ &_paramsCell[iQuad*numParamsQuadPt], numParamsQuadPt,
+ &totalStrain[iQuad*tensorSize], tensorSize);
return _stress;
} // calcStress
// ----------------------------------------------------------------------
// Compute derivative of elasticity matrix for cell at quadrature points.
-const std::vector<pylith::double_array>&
+const pylith::double_array&
pylith::materials::ElasticMaterial::calcDerivElastic(
- const std::vector<double_array>& totalStrain)
+ const double_array& totalStrain)
{ // calcDerivElastic
const int numQuadPts = _numQuadPts;
- assert(_paramsCell.size() == numQuadPts);
- assert(totalStrain.size() == numQuadPts);
- assert(_elasticConsts.size() == numQuadPts);
+ const int numParamsQuadPt = _numParamsQuadPt;
+ const int tensorSize = _tensorSize();
+ const int numElasticConsts = _numElasticConsts();
+ assert(_paramsCell.size() == numQuadPts*numParamsQuadPt);
+ assert(totalStrain.size() == numQuadPts*tensorSize);
+ assert(_elasticConsts.size() == numQuadPts*numElasticConsts);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- _calcElasticConsts(&_elasticConsts[iQuad], _paramsCell[iQuad],
- totalStrain[iQuad]);
+ _calcElasticConsts(&_elasticConsts[iQuad*numElasticConsts], numElasticConsts,
+ &_paramsCell[iQuad*numParamsQuadPt], numParamsQuadPt,
+ &totalStrain[iQuad*tensorSize], tensorSize);
return _elasticConsts;
} // calcDerivElastic
@@ -91,33 +98,21 @@
// ----------------------------------------------------------------------
// Get cell's state variable information from material's sections.
void
-pylith::materials::ElasticMaterial::getStateVarsCell(const Mesh::point_type& cell,
- const int numQuadPts)
+pylith::materials::ElasticMaterial::getStateVarsCell(
+ const Mesh::point_type& cell,
+ const int numQuadPts)
{ // getStateVarsCell
if (_numQuadPts != numQuadPts) {
+ const int numParamsQuadPt = _numParamsQuadPt;
+ const int tensorSize = _tensorSize();
+ const int numElasticConsts = _numElasticConsts();
+
_numQuadPts = numQuadPts;
- _paramsCell.resize(numQuadPts);
- _density.resize(numQuadPts);
- _stress.resize(numQuadPts);
- _elasticConsts.resize(numQuadPts);
-
- const int_array& numParamValues = _getNumParamValues();
- const int numParams = numParamValues.size();
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- _density[iQuad].resize(1);
- _stress[iQuad].resize(_tensorSize());
- _elasticConsts[iQuad].resize(_numElasticConsts());
- _paramsCell[iQuad].resize(numParams);
- for (int iParam=0; iParam < numParams; ++iParam)
- _paramsCell[iQuad][iParam].resize(numParamValues[iParam]);
- } // for
-
- int maxNumValues = 0;
- for (int iParam=0; iParam < numParams; ++iParam)
- if (numParamValues[iParam] > maxNumValues)
- maxNumValues = numParamValues[iParam];
- _parameterCell.resize(numQuadPts*maxNumValues);
+ _paramsCell.resize(numQuadPts*numParamsQuadPt);
+ _density.resize(numQuadPts*1);
+ _stress.resize(numQuadPts*tensorSize);
+ _elasticConsts.resize(numQuadPts*numElasticConsts);
} // if
_getParameters(cell);
@@ -126,42 +121,19 @@
// ----------------------------------------------------------------------
// Update state variables (for next time step).
void
-pylith::materials::ElasticMaterial::updateState(
- const std::vector<double_array>& totalStrain,
- const Mesh::point_type& cell)
+pylith::materials::ElasticMaterial::updateState(const double_array& totalStrain,
+ const Mesh::point_type& cell)
{ // updateState
const int numQuadPts = _numQuadPts;
+ const int numParamsQuadPt = _numParamsQuadPt;
+ const int tensorSize = _tensorSize();
getStateVarsCell(cell, numQuadPts);
- assert(_paramsCell.size() == numQuadPts);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- _updateState(&_paramsCell[iQuad], totalStrain[iQuad]);
-
- const int_array& numParamValues = _getNumParamValues();
- const int numParams = numParamValues.size();
- const char** paramNames = _parameterNames();
+ _updateState(&_paramsCell[iQuad*numParamsQuadPt], numParamsQuadPt,
+ &totalStrain[iQuad*tensorSize], tensorSize);
- for (int iParam=0; iParam < numParams; ++iParam) {
- const ALE::Obj<real_section_type>& parameter =
- _parameters->getReal(paramNames[iParam]);
- assert(!parameter.isNull());
-
- const int numValues = numParamValues[iParam];
- assert(_parameterCell.size() >= numQuadPts*numValues);
-
- assert(parameter->getFiberDimension(cell) == numQuadPts*numValues);
- for (int iQuadPt=0; iQuadPt < numQuadPts; ++iQuadPt)
-#ifdef FASTER
- memcpy(&_parameterCell[iQuadPt*numValues],
- &_paramsCell[iQuadPt][iParam][0],
- numValues*sizeof(double));
-#else
- for (int iValue=0; iValue < numValues; ++iValue)
- _parameterCell[iQuadPt*numValues+iValue] =
- _paramsCell[iQuadPt][iParam][iValue];
-#endif
- parameter->updatePoint(cell, &_parameterCell[0]);
- } // for
+ _parameters->updatePoint(cell, &_paramsCell[0]);
} // updateState
// ----------------------------------------------------------------------
@@ -169,43 +141,25 @@
void
pylith::materials::ElasticMaterial::_getParameters(const Mesh::point_type& cell)
{ // _getParameters
- assert(0 != _parameters);
+ assert(!_parameters.isNull());
const int numQuadPts = _numQuadPts;
- assert(_paramsCell.size() == numQuadPts);
-
- const int_array& numParamValues = _getNumParamValues();
- const int numParams = numParamValues.size();
- const char** paramNames = _parameterNames();
-
- for (int iParam=0; iParam < numParams; ++iParam) {
- const ALE::Obj<real_section_type>& parameter =
- _parameters->getReal(paramNames[iParam]);
- assert(!parameter.isNull());
-
- const int numValues = numParamValues[iParam];
- assert(parameter->getFiberDimension(cell) == numQuadPts*numValues);
- const real_section_type::value_type* parameterCell =
- parameter->restrictPoint(cell);
- for (int iQuadPt=0; iQuadPt < numQuadPts; ++iQuadPt)
-#ifdef FASTER
- memcpy(&_paramsCell[iQuadPt][iParam][0],
- ¶meterCell[iQuadPt*numValues],
- numValues*sizeof(double));
-#else
- for (int iValue=0; iValue < numValues; ++iValue)
- _paramsCell[iQuadPt][iParam][iValue] =
- parameterCell[iQuadPt*numValues+iValue];
-#endif
- } // for
+ const int numParamsQuadPt = _numParamsQuadPt;
+ assert(_paramsCell.size() == numQuadPts*numParamsQuadPt);
+ assert(_parameters->getFiberDimension(cell) == numQuadPts*numParamsQuadPt);
+ const real_section_type::value_type* parameterCell =
+ _parameters->restrictPoint(cell);
+ memcpy(&_paramsCell[0], parameterCell,
+ numQuadPts*numParamsQuadPt*sizeof(double));
} // _getParameters
// ----------------------------------------------------------------------
// Update parameters (for next time step).
void
-pylith::materials::ElasticMaterial::_updateState(
- std::vector<double_array>* const parameters,
- const double_array& totalStrain)
+pylith::materials::ElasticMaterial::_updateState(double* const parameters,
+ const int numParamsQuadPt,
+ const double* totalStrain,
+ const int strainSize)
{ // _updateState
} // _updateState
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -65,7 +65,7 @@
*
* @returns Array of density values at cell's quadrature points.
*/
- const std::vector<double_array>& calcDensity(void);
+ const double_array& calcDensity(void);
/** Get stress tensor at quadrature points.
*
@@ -85,8 +85,8 @@
*
* @returns Array of stresses at cell's quadrature points.
*/
- const std::vector<double_array>&
- calcStress(const std::vector<double_array>& totalStrain);
+ const double_array&
+ calcStress(const double_array& totalStrain);
/** Compute derivative of elasticity matrix for cell at quadrature points.
*
@@ -111,8 +111,8 @@
* @param totalStrain Total strain tensor at quadrature points
* [numQuadPts][tensorSize]
*/
- const std::vector<double_array>&
- calcDerivElastic(const std::vector<double_array>& totalStrain);
+ const double_array&
+ calcDerivElastic(const double_array& totalStrain);
/** Update state variables (for next time step).
*
@@ -120,7 +120,7 @@
* [numQuadPts][tensorSize]
* @param cell Finite element cell
*/
- void updateState(const std::vector<double_array>& totalStrain,
+ void updateState(const double_array& totalStrain,
const Mesh::point_type& cell);
/** Set whether elastic or inelastic constitutive relations are used.
@@ -157,43 +157,61 @@
/** Compute density from parameters.
*
- * @param density Array for density
- * @param parameters Parameters at location
+ * @param density Array for density.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
*/
virtual
- void _calcDensity(double_array* const density,
- const std::vector<double_array>& parameters) = 0;
+ void _calcDensity(double* const density,
+ const double* parameters,
+ const int numParams) = 0;
/** Compute stress tensor from parameters.
*
- * @param stress Array for stress tensor
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
* @param parameters Parameters at location.
+ * @param numParams Number of parameters.
* @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
virtual
- void _calcStress(double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain) = 0;
+ void _calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize) = 0;
/** Compute derivatives of elasticity matrix from parameters.
*
- * @param elasticConsts Array for elastic constants
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
* @param parameters Parameters at location.
+ * @param numParams Number of parameters.
* @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
virtual
- void _calcElasticConsts(double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain) = 0;
+ void _calcElasticConsts(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize) = 0;
/** Update parameters (for next time step).
*
* @param parameters Parameters at location.
+ * @param numParams Number of parameters.
* @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
virtual
- void _updateState(std::vector<double_array>* const parameters,
- const double_array& totalStrain);
+ void _updateState(double* const parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
@@ -220,38 +238,31 @@
/** Parameters at quadrature points for current cell.
*
- * size = [numQuadPts][numParams][numValues]
- * index = [iQuadPt][iParam][iValue]
+ * size = numQuadPts*numParamsQuadPt
+ * index = iQuadPt*iParam*iValue
*/
- std::vector<std::vector<double_array> > _paramsCell;
+ double_array _paramsCell;
- /** Single parameter at quadrature points for current cell.
- *
- * size = [numQuadPts*numValues]
- * index = [iQuadPt*numValues+iValue]
- */
- double_array _parameterCell;
-
/** Density value at quadrature points for current cell.
*
- * size = [numQuadPts][1]
- * index = [iQuadPt][0]
+ * size = numQuadPts
+ * index = iQuadPt
*/
- std::vector<double_array> _density;
+ double_array _density;
/** Stress tensor at quadrature points for current cell.
*
- * size = [numQuadPts][tensorSize]
- * index = [iQuadPt][iStress]
+ * size = numQuadPts*tensorSize
+ * index = *iQuadPt*tensorSize + iStress
*/
- std::vector<double_array> _stress;
+ double_array _stress;
/** Elasticity matrix at quadrature points for current cell.
*
- * size = [numQuadPts][numElasticConsts]
- * index = [iQuadPt][iConstant]
+ * size = numQuadPts*numElasticConsts
+ * index = iQuadPt*numElasticConsts+iConstant
*/
- std::vector<double_array> _elasticConsts;
+ double_array _elasticConsts;
}; // class ElasticMaterial
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -14,8 +14,6 @@
#include "ElasticPlaneStrain.hh" // implementation of object methods
-#include "pylith/utils/array.hh" // USES double_array
-
#include "petsc.h" // USES PetscLogFlopsNoCheck
#include <assert.h> // USES assert()
@@ -35,7 +33,6 @@
/// Values expected in spatial database
const int numDBValues = 3;
- const int numParamValues[] = { 1, 1, 1 };
const char* namesDBValues[] = { "density", "vs", "vp" };
/// Indices (order) of database values
@@ -45,12 +42,12 @@
/// Parameters
const int numParameters = 3;
- const char* namesParameters[] = { "density", "mu", "lambda" };
+ const int numParamValues[] = { 1, 1, 1 };
/// Indices (order) of parameters
const int pidDensity = 0;
- const int pidMu = 1;
- const int pidLambda = 2;
+ const int pidMu = pidDensity + 1;
+ const int pidLambda = pidMu + 1;
} // _ElasticPlaneStrain
} // materials
@@ -89,27 +86,17 @@
} // _numDBValues
// ----------------------------------------------------------------------
-// Get names of parameters for physical properties.
-const char**
-pylith::materials::ElasticPlaneStrain::_parameterNames(void) const
-{ // _parameterNames
- return _ElasticPlaneStrain::namesParameters;
-} // _parameterNames
-
-// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
pylith::materials::ElasticPlaneStrain::_dbToParameters(
- std::vector<double_array>* const paramVals,
- const double_array& dbValues) const
+ double* const paramVals,
+ const int numParams,
+ const double_array& dbValues) const
{ // computeParameters
assert(0 != paramVals);
- const int numParams = paramVals->size();
- assert(_ElasticPlaneStrain::numParameters == numParams);
+ assert(_numParamsQuadPt == numParams);
const int numDBValues = dbValues.size();
assert(_ElasticPlaneStrain::numDBValues == numDBValues);
- for (int i=0; i < numParams; ++i)
- assert(_ElasticPlaneStrain::numParamValues[i] == (*paramVals)[i].size());
const double density = dbValues[_ElasticPlaneStrain::didDensity];
const double vs = dbValues[_ElasticPlaneStrain::didVs];
@@ -126,9 +113,9 @@
throw std::runtime_error(msg.str());
} // if
- (*paramVals)[_ElasticPlaneStrain::pidDensity][0] = density;
- (*paramVals)[_ElasticPlaneStrain::pidMu][0] = mu;
- (*paramVals)[_ElasticPlaneStrain::pidLambda][0] = lambda;
+ paramVals[_ElasticPlaneStrain::pidDensity] = density;
+ paramVals[_ElasticPlaneStrain::pidMu] = mu;
+ paramVals[_ElasticPlaneStrain::pidLambda] = lambda;
PetscLogFlopsNoCheck(6);
} // computeParameters
@@ -153,36 +140,38 @@
// Compute density at location from parameters.
void
pylith::materials::ElasticPlaneStrain::_calcDensity(
- double_array* const density,
- const std::vector<double_array>& parameters)
+ double* const density,
+ const double* parameters,
+ const int numParams)
{ // calcDensity
assert(0 != density);
- assert(1 == density->size());
- assert(_ElasticPlaneStrain::numParameters == parameters.size());
- assert(1 == parameters[_ElasticPlaneStrain::pidDensity].size());
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
- (*density)[0] = parameters[_ElasticPlaneStrain::pidDensity][0];
+ density[0] = parameters[_ElasticPlaneStrain::pidDensity];
} // calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor at location from parameters.
void
pylith::materials::ElasticPlaneStrain::_calcStress(
- double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+ double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _calcStress
assert(0 != stress);
- assert(_ElasticPlaneStrain::tensorSize == stress->size());
- assert(_ElasticPlaneStrain::numParameters == parameters.size());
- assert(_ElasticPlaneStrain::tensorSize == totalStrain.size());
- assert(1 == parameters[_ElasticPlaneStrain::pidDensity].size());
- assert(1 == parameters[_ElasticPlaneStrain::pidMu].size());
- assert(1 == parameters[_ElasticPlaneStrain::pidLambda].size());
+ assert(_ElasticPlaneStrain::tensorSize == stressSize);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_ElasticPlaneStrain::tensorSize == strainSize);
- const double density = parameters[_ElasticPlaneStrain::pidDensity][0];
- const double mu = parameters[_ElasticPlaneStrain::pidMu][0];
- const double lambda = parameters[_ElasticPlaneStrain::pidLambda][0];
+ const double density = parameters[_ElasticPlaneStrain::pidDensity];
+ const double mu = parameters[_ElasticPlaneStrain::pidMu];
+ const double lambda = parameters[_ElasticPlaneStrain::pidLambda];
const double mu2 = 2.0 * mu;
@@ -192,9 +181,9 @@
const double s12 = lambda * (e11 + e22);
- (*stress)[0] = s12 + mu2*e11;
- (*stress)[1] = s12 + mu2*e22;
- (*stress)[2] = mu2 * e12;
+ stress[0] = s12 + mu2*e11;
+ stress[1] = s12 + mu2*e22;
+ stress[2] = mu2 * e12;
PetscLogFlopsNoCheck(8);
} // _calcStress
@@ -203,31 +192,33 @@
// Compute density at location from parameters.
void
pylith::materials::ElasticPlaneStrain::_calcElasticConsts(
- double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // calcElasticConsts
assert(0 != elasticConsts);
- assert(_ElasticPlaneStrain::numElasticConsts == elasticConsts->size());
- assert(_ElasticPlaneStrain::numParameters == parameters.size());
- assert(_ElasticPlaneStrain::tensorSize == totalStrain.size());
- assert(1 == parameters[_ElasticPlaneStrain::pidDensity].size());
- assert(1 == parameters[_ElasticPlaneStrain::pidMu].size());
- assert(1 == parameters[_ElasticPlaneStrain::pidLambda].size());
+ assert(_ElasticPlaneStrain::numElasticConsts == numElasticConsts);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_ElasticPlaneStrain::tensorSize == strainSize);
- const double density = parameters[_ElasticPlaneStrain::pidDensity][0];
- const double mu = parameters[_ElasticPlaneStrain::pidMu][0];
- const double lambda = parameters[_ElasticPlaneStrain::pidLambda][0];
+ const double density = parameters[_ElasticPlaneStrain::pidDensity];
+ const double mu = parameters[_ElasticPlaneStrain::pidMu];
+ const double lambda = parameters[_ElasticPlaneStrain::pidLambda];
const double mu2 = 2.0 * mu;
const double lambda2mu = lambda + mu2;
- (*elasticConsts)[0] = lambda2mu; // C1111
- (*elasticConsts)[1] = lambda; // C1122
- (*elasticConsts)[2] = 0; // C1112
- (*elasticConsts)[3] = lambda2mu; // C2222
- (*elasticConsts)[4] = 0; // C2212
- (*elasticConsts)[5] = mu2; // C1212
+ elasticConsts[0] = lambda2mu; // C1111
+ elasticConsts[1] = lambda; // C1122
+ elasticConsts[2] = 0; // C1112
+ elasticConsts[3] = lambda2mu; // C2222
+ elasticConsts[4] = 0; // C2212
+ elasticConsts[5] = mu2; // C1212
PetscLogFlopsNoCheck(2);
} // calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -66,12 +66,6 @@
*/
int _numDBValues(void) const;
- /** Get names of parameters for physical properties.
- *
- * @returns Names of parameters
- */
- const char** _parameterNames(void) const;
-
/** Compute parameters from values in spatial database.
*
* Order of values in arrays matches order used in dbValues() and
@@ -80,7 +74,8 @@
* @param paramVals Array of parameters
* @param dbValues Array of database values
*/
- void _dbToParameters(std::vector<double_array>* const paramVals,
+ void _dbToParameters(double* const paramVals,
+ const int numParams,
const double_array& dbValues) const;
/** Get number of entries in stress tensor.
@@ -105,31 +100,45 @@
/** Compute density from parameters.
*
- * @param density Array for density
- * @param parameters Parameters at location
+ * @param density Array for density.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
*/
- void _calcDensity(double_array* const density,
- const std::vector<double_array>& parameters);
+ void _calcDensity(double* const density,
+ const double* parameters,
+ const int numParams);
/** Compute stress tensor from parameters.
*
- * @param stress Array for stress tensor
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcStress(double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
/** Compute derivatives of elasticity matrix from parameters.
*
- * @param elasticConsts Array for elastic constants
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcElasticConsts(double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcElasticConsts(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -14,8 +14,6 @@
#include "ElasticPlaneStress.hh" // implementation of object methods
-#include "pylith/utils/array.hh" // USES double_array
-
#include "petsc.h" // USES PetscLogFlopsNoCheck
#include <assert.h> // USES assert()
@@ -45,12 +43,11 @@
/// Parameters
const int numParameters = 3;
const int numParamValues[] = { 1, 1, 1 };
- const char* namesParameters[] = { "density", "mu", "lambda" };
/// Indices (order) of parameters
const int pidDensity = 0;
- const int pidMu = 1;
- const int pidLambda = 2;
+ const int pidMu = pidDensity + 1;
+ const int pidLambda = pidMu + 1;
} // _ElasticPlaneStress
} // materials
@@ -88,27 +85,17 @@
} // _numDBValues
// ----------------------------------------------------------------------
-// Get names of parameters for physical properties.
-const char**
-pylith::materials::ElasticPlaneStress::_parameterNames(void) const
-{ // _parameterNames
- return _ElasticPlaneStress::namesParameters;
-} // _parameterNames
-
-// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
pylith::materials::ElasticPlaneStress::_dbToParameters(
- std::vector<double_array>* paramVals,
+ double* paramVals,
+ const int numParams,
const double_array& dbValues) const
{ // _dbToParameters
assert(0 != paramVals);
- const int numParams = paramVals->size();
- assert(_ElasticPlaneStress::numParameters == numParams);
+ assert(_numParamsQuadPt == numParams);
const int numDBValues = dbValues.size();
assert(_ElasticPlaneStress::numDBValues == numDBValues);
- for (int i=0; i < numParams; ++i)
- assert(_ElasticPlaneStress::numParamValues[i] == (*paramVals)[i].size());
const double density = dbValues[_ElasticPlaneStress::didDensity];
const double vs = dbValues[_ElasticPlaneStress::didVs];
@@ -125,9 +112,9 @@
throw std::runtime_error(msg.str());
} // if
- (*paramVals)[_ElasticPlaneStress::pidDensity][0] = density;
- (*paramVals)[_ElasticPlaneStress::pidMu][0] = mu;
- (*paramVals)[_ElasticPlaneStress::pidLambda][0] = lambda;
+ paramVals[_ElasticPlaneStress::pidDensity] = density;
+ paramVals[_ElasticPlaneStress::pidMu] = mu;
+ paramVals[_ElasticPlaneStress::pidLambda] = lambda;
PetscLogFlopsNoCheck(6);
} // _dbToParameters
@@ -151,36 +138,37 @@
// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::ElasticPlaneStress::_calcDensity(double_array* const density,
- const std::vector<double_array>& parameters)
+pylith::materials::ElasticPlaneStress::_calcDensity(double* const density,
+ const double* parameters,
+ const int numParams)
{ // calcDensity
assert(0 != density);
- assert(1 == density->size());
- assert(_ElasticPlaneStress::numParameters == parameters.size());
- assert(1 == parameters[_ElasticPlaneStress::pidDensity].size());
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
- (*density)[0] = parameters[_ElasticPlaneStress::pidDensity][0];
+ density[0] = parameters[_ElasticPlaneStress::pidDensity];
} // calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor at location from parameters.
void
-pylith::materials::ElasticPlaneStress::_calcStress(
- double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+pylith::materials::ElasticPlaneStress::_calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _calcStress
assert(0 != stress);
- assert(_ElasticPlaneStress::tensorSize == stress->size());
- assert(_ElasticPlaneStress::numParameters == parameters.size());
- assert(_ElasticPlaneStress::tensorSize == totalStrain.size());
- assert(1 == parameters[_ElasticPlaneStress::pidDensity].size());
- assert(1 == parameters[_ElasticPlaneStress::pidMu].size());
- assert(1 == parameters[_ElasticPlaneStress::pidLambda].size());
+ assert(_ElasticPlaneStress::tensorSize == stressSize);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_ElasticPlaneStress::tensorSize == strainSize);
- const double density = parameters[_ElasticPlaneStress::pidDensity][0];
- const double mu = parameters[_ElasticPlaneStress::pidMu][0];
- const double lambda = parameters[_ElasticPlaneStress::pidLambda][0];
+ const double density = parameters[_ElasticPlaneStress::pidDensity];
+ const double mu = parameters[_ElasticPlaneStress::pidMu];
+ const double lambda = parameters[_ElasticPlaneStress::pidLambda];
const double mu2 = 2.0 * mu;
const double lambda2mu = lambda + mu2;
@@ -190,9 +178,9 @@
const double e22 = totalStrain[1];
const double e12 = totalStrain[2];
- (*stress)[0] = (2.0*mu2*lambdamu * e11 + mu2*lambda * e22)/lambda2mu;
- (*stress)[1] = (mu2*lambda * e11 + 2.0*mu2*lambdamu * e22)/lambda2mu;
- (*stress)[2] = mu2 * e12;
+ stress[0] = (2.0*mu2*lambdamu * e11 + mu2*lambda * e22) / lambda2mu;
+ stress[1] = (mu2*lambda * e11 + 2.0*mu2*lambdamu * e22) / lambda2mu;
+ stress[2] = mu2 * e12;
PetscLogFlopsNoCheck(18);
} // _calcStress
@@ -201,32 +189,34 @@
// Compute density at location from parameters.
void
pylith::materials::ElasticPlaneStress::_calcElasticConsts(
- double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // calcElasticConsts
assert(0 != elasticConsts);
- assert(_ElasticPlaneStress::numElasticConsts == elasticConsts->size());
- assert(_ElasticPlaneStress::numParameters == parameters.size());
- assert(_ElasticPlaneStress::tensorSize == totalStrain.size());
- assert(1 == parameters[_ElasticPlaneStress::pidDensity].size());
- assert(1 == parameters[_ElasticPlaneStress::pidMu].size());
- assert(1 == parameters[_ElasticPlaneStress::pidLambda].size());
+ assert(_ElasticPlaneStress::numElasticConsts == numElasticConsts);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_ElasticPlaneStress::tensorSize == strainSize);
- const double density = parameters[_ElasticPlaneStress::pidDensity][0];
- const double mu = parameters[_ElasticPlaneStress::pidMu][0];
- const double lambda = parameters[_ElasticPlaneStress::pidLambda][0];
+ const double density = parameters[_ElasticPlaneStress::pidDensity];
+ const double mu = parameters[_ElasticPlaneStress::pidMu];
+ const double lambda = parameters[_ElasticPlaneStress::pidLambda];
const double mu2 = 2.0 * mu;
const double lambda2mu = lambda + mu2;
const double c11 = 2.0 * mu2 * (lambda + mu) / lambda2mu;
- (*elasticConsts)[0] = c11; // C1111
- (*elasticConsts)[1] = mu2 * lambda / lambda2mu; // C1122
- (*elasticConsts)[2] = 0; // C1112
- (*elasticConsts)[3] = c11; // C2222
- (*elasticConsts)[4] = 0; // C2212
- (*elasticConsts)[5] = mu2; // C1212
+ elasticConsts[0] = c11; // C1111
+ elasticConsts[1] = mu2 * lambda / lambda2mu; // C1122
+ elasticConsts[2] = 0; // C1112
+ elasticConsts[3] = c11; // C2222
+ elasticConsts[4] = 0; // C2212
+ elasticConsts[5] = mu2; // C1212
PetscLogFlopsNoCheck(8);
} // calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -66,21 +66,17 @@
*/
int _numDBValues(void) const;
- /** Get names of parameters for physical properties.
- *
- * @returns Names of parameters
- */
- const char** _parameterNames(void) const;
-
/** Compute parameters from values in spatial database.
*
* Order of values in arrays matches order used in dbValues() and
* parameterNames().
*
* @param paramVals Array of parameters
+ * @param numParams Number of parameters at quadrature point.
* @param dbValues Array of database values
*/
- void _dbToParameters(std::vector<double_array>* const paramVals,
+ void _dbToParameters(double* const paramVals,
+ const int numParams,
const double_array& dbValues) const;
/** Get number of entries in stress tensor.
@@ -105,31 +101,45 @@
/** Compute density from parameters.
*
- * @param density Array for density
- * @param parameters Parameters at location
+ * @param density Array for density.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
*/
- void _calcDensity(double_array* const density,
- const std::vector<double_array>& parameters);
+ void _calcDensity(double* const density,
+ const double* parameters,
+ const int numParams);
/** Compute stress tensor from parameters.
*
- * @param stress Array for stress tensor
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcStress(double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
/** Compute derivatives of elasticity matrix from parameters.
*
- * @param elasticConsts Array for elastic constants
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcElasticConsts(double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcElasticConsts(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -14,8 +14,6 @@
#include "ElasticStrain1D.hh" // implementation of object methods
-#include "pylith/utils/array.hh" // USES double_array
-
#include "petsc.h" // USES PetscLogFlopsNoCheck
#include <assert.h> // USES assert()
@@ -35,18 +33,17 @@
const int numDBValues = 2;
const char* namesDBValues[] = { "density", "vp" };
- /// Indices (order) of database values
+ /// Indices of database values
const int didDensity = 0;
const int didVp = 1;
/// Parameters
const int numParameters = 2;
const int numParamValues[] = { 1, 1 };
- const char* namesParameters[] = { "density", "lambda2mu" };
- /// Indices (order) of parameters
+ /// Indices of parameters
const int pidDensity = 0;
- const int pidLambda2mu = 1;
+ const int pidLambda2mu = pidDensity + 1;
} // _ElasticStrain1D
} // materials
} // pylith
@@ -84,33 +81,24 @@
} // _numDBValues
// ----------------------------------------------------------------------
-// Get names of parameters for physical properties.
-const char**
-pylith::materials::ElasticStrain1D::_parameterNames(void) const
-{ // _parameterNames
- return _ElasticStrain1D::namesParameters;
-} // _parameterNames
-
-// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
-pylith::materials::ElasticStrain1D::_dbToParameters(std::vector<double_array>* const paramVals,
- const double_array& dbValues) const
+pylith::materials::ElasticStrain1D::_dbToParameters(
+ double* const paramVals,
+ const int numParams,
+ const double_array& dbValues) const
{ // _dbToParameters
assert(0 != paramVals);
- const int numParams = paramVals->size();
- assert(_ElasticStrain1D::numParameters == numParams);
+ assert(_numParamsQuadPt == numParams);
const int numDBValues = dbValues.size();
assert(_ElasticStrain1D::numDBValues == numDBValues);
- for (int i=0; i < numParams; ++i)
- assert(_ElasticStrain1D::numParamValues[i] == (*paramVals)[i].size());
const double density = dbValues[_ElasticStrain1D::didDensity];
const double vp = dbValues[_ElasticStrain1D::didVp];
const double lambda2mu = density * vp*vp;
- (*paramVals)[_ElasticStrain1D::pidDensity][0] = density;
- (*paramVals)[_ElasticStrain1D::pidLambda2mu][0] = lambda2mu;
+ paramVals[_ElasticStrain1D::pidDensity] = density;
+ paramVals[_ElasticStrain1D::pidLambda2mu] = lambda2mu;
PetscLogFlopsNoCheck(2);
} // _dbToParameters
@@ -134,37 +122,40 @@
// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::ElasticStrain1D::_calcDensity(double_array* const density,
- const std::vector<double_array>& parameters)
+pylith::materials::ElasticStrain1D::_calcDensity(double* const density,
+ const double* parameters,
+ const int numParams)
{ // _calcDensity
assert(0 != density);
- assert(1 == density->size());
- assert(_ElasticStrain1D::numParameters == parameters.size());
- assert(1 == parameters[_ElasticStrain1D::pidDensity].size());
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
- (*density)[0] = parameters[_ElasticStrain1D::pidDensity][0];
+ density[0] = parameters[_ElasticStrain1D::pidDensity];
} // _calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor at location from parameters.
void
pylith::materials::ElasticStrain1D::_calcStress(
- double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+ double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _calcStress
assert(0 != stress);
- assert(_ElasticStrain1D::tensorSize == stress->size());
- assert(_ElasticStrain1D::numParameters == parameters.size());
- assert(_ElasticStrain1D::tensorSize == totalStrain.size());
- assert(1 == parameters[_ElasticStrain1D::pidDensity].size());
- assert(1 == parameters[_ElasticStrain1D::pidLambda2mu].size());
+ assert(_ElasticStrain1D::tensorSize == stressSize);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_ElasticStrain1D::tensorSize == strainSize);
- const double density = parameters[_ElasticStrain1D::pidDensity][0];
- const double lambda2mu = parameters[_ElasticStrain1D::pidLambda2mu][0];
+ 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;
PetscLogFlopsNoCheck(1);
} // _calcStress
@@ -173,21 +164,24 @@
// Compute derivative of elasticity matrix at location from parameters.
void
pylith::materials::ElasticStrain1D::_calcElasticConsts(
- double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _calcElasticConsts
assert(0 != elasticConsts);
- assert(_ElasticStrain1D::numElasticConsts == elasticConsts->size());
- assert(_ElasticStrain1D::numParameters == parameters.size());
- assert(_ElasticStrain1D::tensorSize == totalStrain.size());
- assert(1 == parameters[_ElasticStrain1D::pidDensity].size());
- assert(1 == parameters[_ElasticStrain1D::pidLambda2mu].size());
+ assert(_ElasticStrain1D::numElasticConsts == numElasticConsts);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_ElasticStrain1D::tensorSize == strainSize);
- const double density = parameters[_ElasticStrain1D::pidDensity][0];
- const double lambda2mu = parameters[_ElasticStrain1D::pidLambda2mu][0];
+ 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 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -65,21 +65,17 @@
*/
int _numDBValues(void) const;
- /** Get names of parameters for physical properties.
- *
- * @returns Names of parameters
- */
- const char** _parameterNames(void) const;
-
/** Compute parameters from values in spatial database.
*
* Order of values in arrays matches order used in dbValues() and
* parameterNames().
*
* @param paramVals Array of parameters
+ * @param numParams Number of parameters at quadrature point.
* @param dbValues Array of database values
*/
- void _dbToParameters(std::vector<double_array>* const paramVals,
+ void _dbToParameters(double* const paramVals,
+ const int numParams,
const double_array& dbValues) const;
/** Get number of entries in stress tensor.
@@ -104,31 +100,45 @@
/** Compute density from parameters.
*
- * @param density Array for density
- * @param parameters Parameters at location
+ * @param density Array for density.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
*/
- void _calcDensity(double_array* const density,
- const std::vector<double_array>& parameters);
+ void _calcDensity(double* const density,
+ const double* parameters,
+ const int numParams);
/** Compute stress tensor from parameters.
*
- * @param stress Array for stress tensor
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcStress(double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
/** Compute derivatives of elasticity matrix from parameters.
*
- * @param elasticConsts Array for elastic constants
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcElasticConsts(double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcElasticConsts(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -14,8 +14,6 @@
#include "ElasticStress1D.hh" // implementation of object methods
-#include "pylith/utils/array.hh" // USES double_array
-
#include "petsc.h" // USES PetscLogFlopsNoCheck
#include <assert.h> // USES assert()
@@ -43,13 +41,12 @@
/// Parameters
const int numParameters = 3;
- const char* namesParameters[] = {"density", "mu", "lambda" };
const int numParamValues[] = { 1, 1, 1 };
-
- /// Indices (order) of parameters
+
+ /// Indices of parameters
const int pidDensity = 0;
- const int pidMu = 1;
- const int pidLambda = 2;
+ const int pidMu = pidDensity + 1;
+ const int pidLambda = pidMu + 1;
} // _ElasticStress1D
} // materials
@@ -87,26 +84,17 @@
} // _numDBValues
// ----------------------------------------------------------------------
-// Get names of parameters for physical properties.
-const char**
-pylith::materials::ElasticStress1D::_parameterNames(void) const
-{ // _parameterNames
- return _ElasticStress1D::namesParameters;
-} // _parameterNames
-
-// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
-pylith::materials::ElasticStress1D::_dbToParameters(std::vector<double_array>* const paramVals,
- const double_array& dbValues) const
+pylith::materials::ElasticStress1D::_dbToParameters(
+ double* const paramVals,
+ const int numParams,
+ const double_array& dbValues) const
{ // computeParameters
assert(0 != paramVals);
- const int numParams = paramVals->size();
- assert(_ElasticStress1D::numParameters == numParams);
+ assert(_numParamsQuadPt == numParams);
const int numDBValues = dbValues.size();
assert(_ElasticStress1D::numDBValues == numDBValues);
- for (int i=0; i < numParams; ++i)
- assert(_ElasticStress1D::numParamValues[i] == (*paramVals)[i].size());
const double density = dbValues[_ElasticStress1D::didDensity];
const double vs = dbValues[_ElasticStress1D::didVs];
@@ -115,9 +103,9 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- (*paramVals)[_ElasticStress1D::pidDensity][0] = density;
- (*paramVals)[_ElasticStress1D::pidMu][0] = mu;
- (*paramVals)[_ElasticStress1D::pidLambda][0] = lambda;
+ paramVals[_ElasticStress1D::pidDensity] = density;
+ paramVals[_ElasticStress1D::pidMu] = mu;
+ paramVals[_ElasticStress1D::pidLambda] = lambda;
PetscLogFlopsNoCheck(6);
} // computeParameters
@@ -141,35 +129,40 @@
// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::ElasticStress1D::_calcDensity(double_array* const density,
- const std::vector<double_array>& parameters)
+pylith::materials::ElasticStress1D::_calcDensity(double* const density,
+ const double* parameters,
+ const int numParams)
{ // _calcDensity
assert(0 != density);
- assert(1 == density->size());
- assert(_ElasticStress1D::numParameters == parameters.size());
- assert(1 == parameters[_ElasticStress1D::pidDensity].size());
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
- (*density)[0] = parameters[_ElasticStress1D::pidDensity][0];
+ density[0] = parameters[_ElasticStress1D::pidDensity];
} // _calcDensity
// ----------------------------------------------------------------------
// Compute stress tensor at location from parameters.
void
-pylith::materials::ElasticStress1D::_calcStress(double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+pylith::materials::ElasticStress1D::_calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _calcStress
assert(0 != stress);
- assert(_ElasticStress1D::tensorSize == stress->size());
- assert(_ElasticStress1D::numParameters == parameters.size());
- assert(_ElasticStress1D::tensorSize == totalStrain.size());
+ assert(_ElasticStress1D::tensorSize == stressSize);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_ElasticStress1D::tensorSize == strainSize);
- const double density = parameters[_ElasticStress1D::pidDensity][0];
- const double mu = parameters[_ElasticStress1D::pidMu][0];
- const double lambda = parameters[_ElasticStress1D::pidLambda][0];
+ 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;
PetscLogFlopsNoCheck(7);
} // _calcStress
@@ -178,20 +171,25 @@
// Compute derivative of elasticity matrix at location from parameters.
void
pylith::materials::ElasticStress1D::_calcElasticConsts(
- double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _calcElasticConsts
assert(0 != elasticConsts);
- assert(_ElasticStress1D::numElasticConsts == elasticConsts->size());
- assert(_ElasticStress1D::numParameters == parameters.size());
- assert(_ElasticStress1D::tensorSize == totalStrain.size());
+ assert(_ElasticStress1D::numElasticConsts == numElasticConsts);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_ElasticStress1D::tensorSize == strainSize);
- const double density = parameters[_ElasticStress1D::pidDensity][0];
- const double mu = parameters[_ElasticStress1D::pidMu][0];
- const double lambda = parameters[_ElasticStress1D::pidLambda][0];
+ 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);
PetscLogFlopsNoCheck(6);
} // _calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -66,21 +66,17 @@
*/
int _numDBValues(void) const;
- /** Get names of parameters for physical properties.
- *
- * @returns Names of parameters
- */
- const char** _parameterNames(void) const;
-
/** Compute parameters from values in spatial database.
*
* Order of values in arrays matches order used in dbValues() and
* parameterNames().
*
* @param paramVals Array of parameters
+ * @param numParams Number of parameters at quadrature point.
* @param dbValues Array of database values
*/
- void _dbToParameters(std::vector<double_array>* const paramVals,
+ void _dbToParameters(double* const paramVals,
+ const int numParams,
const double_array& dbValues) const;
/** Get number of entries in stress tensor.
@@ -105,31 +101,45 @@
/** Compute density from parameters.
*
- * @param density Array for density
- * @param parameters Parameters at location
+ * @param density Array for density.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
*/
- void _calcDensity(double_array* const density,
- const std::vector<double_array>& parameters);
+ void _calcDensity(double* const density,
+ const double* parameters,
+ const int numParams);
/** Compute stress tensor from parameters.
*
- * @param stress Array for stress tensor
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcStress(double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
/** Compute derivatives of elasticity matrix from parameters.
*
- * @param elasticConsts Array for elastic constants
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcElasticConsts(double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcElasticConsts(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -31,7 +31,7 @@
pylith::materials::Material::Material(const int* numParamValues,
const int size) :
_dt(0.0),
- _parameters(0),
+ _numParamsQuadPt(0),
_dimension(0),
_needNewJacobian(false),
_db(0),
@@ -39,6 +39,9 @@
_label(""),
_numParamValues(numParamValues, size)
{ // constructor
+ for (int i=0; i < size; ++i)
+ _numParamsQuadPt += numParamValues[i];
+ assert(_numParamsQuadPt >= 0);
} // constructor
// ----------------------------------------------------------------------
@@ -47,8 +50,6 @@
{ // destructor
// Python db object owns database, so just set pointer to null
_db = 0;
-
- if (_parameters) {delete _parameters;}; _parameters = 0;
} // destructor
// ----------------------------------------------------------------------
@@ -71,102 +72,37 @@
assert(!cells.isNull());
const ALE::Mesh::label_sequence::iterator cellsEnd = cells->end();
- // Check to make sure we have cells
-#if 0
- if (0 == cells->size()) {
- std::ostringstream msg;
- msg << "Could not find any cells for material '" << _label << "'.";
- throw std::runtime_error(msg.str());
- } // if
-#endif
-
// Create sections to hold parameters for physical properties
- if (0 != _parameters) // Can't delete NULL pointer that holds reference
- delete _parameters;
- _parameters = new topology::FieldsManager(mesh);
- assert(0 != _parameters);
+ _parameters = new real_section_type(mesh->comm(), mesh->debug());
+ assert(!_parameters.isNull());
const int numQuadPts = quadrature->numQuadPts();
const int spaceDim = quadrature->spaceDim();
- const int tensorSize = (3 == spaceDim) ? 6 : ((2 == spaceDim) ? 4 : 1);
const int numParams = _numParamValues.size();
- const char** paramNames = _parameterNames();
- std::vector<ALE::Obj<real_section_type> > paramSections(numParams);
-
- // Reduce memory storage by reusing layouts from basic sections we expect
- int indexScalarLayout = -1;
- int indexVectorLayout = -1;
- int indexTensorLayout = -1;
- for (int iParam=0; iParam < numParams; ++iParam) {
- _parameters->addReal(paramNames[iParam]);
- paramSections[iParam] = _parameters->getReal(paramNames[iParam]);
- assert(!paramSections[iParam].isNull());
+ // Fiber dimension is number of quadrature points times number of
+ // values per parameter
+ const int numParamsQuadPt = _numParamsQuadPt;
+ const int fiberDim = numParamsQuadPt * numQuadPts;
+ _parameters->setFiberDimension(cells, fiberDim);
+ mesh->allocate(_parameters);
- // Fiber dimension is number of quadrature points times number of
- // values per parameter
- const int fiberDim = numQuadPts * _numParamValues[iParam];
- if (1 == _numParamValues[iParam])
- if (-1 != indexScalarLayout) {
- paramSections[iParam]->setAtlas(paramSections[indexScalarLayout]->getAtlas());
- paramSections[iParam]->allocateStorage();
- } else {
- indexScalarLayout = iParam;
- paramSections[iParam]->setFiberDimension(cells, fiberDim);
- mesh->allocate(paramSections[iParam]);
- } // if/else
- else if (spaceDim == _numParamValues[iParam]) {
- if (-1 != indexVectorLayout) {
- paramSections[iParam]->setAtlas(paramSections[indexVectorLayout]->getAtlas());
- paramSections[iParam]->allocateStorage();
- } else {
- indexVectorLayout = iParam;
- paramSections[iParam]->setFiberDimension(cells, fiberDim);
- mesh->allocate(paramSections[iParam]);
- } // if/else
- } // if/else
- else if (tensorSize == _numParamValues[iParam]) {
- if (-1 != indexTensorLayout) {
- paramSections[iParam]->setAtlas(paramSections[indexTensorLayout]->getAtlas());
- paramSections[iParam]->allocateStorage();
- } else {
- indexTensorLayout = iParam;
- paramSections[iParam]->setFiberDimension(cells, fiberDim);
- mesh->allocate(paramSections[iParam]);
- } // if/else
- } else {
- paramSections[iParam]->setFiberDimension(cells, fiberDim);
- mesh->allocate(paramSections[iParam]);
- } // if/else
- } // for
-
// Setup database for querying
const int numValues = _numDBValues();
_db->open();
_db->queryVals(_dbValues(), numValues);
- // Loop over cells
-
// Container for data returned in query of database
double_array queryData(numValues);
- // Container for parameters at a quadrature point
- std::vector<double_array> paramData(numParams);
+ // Container of parameters at cell's quadrature points.
+ double_array cellData(fiberDim);
- // Container of parameter data for a given cell (quadpts + parameters)
- std::vector<double_array> cellData(numParams);
-
- for (int iParam = 0; iParam < numParams; ++iParam) {
- const int fiberDim = numQuadPts * _numParamValues[iParam];
- cellData[iParam].resize(fiberDim);
- paramData[iParam].resize(_numParamValues[iParam]);
- } // for
-
- for (ALE::Mesh::label_sequence::iterator cellIter=cells->begin();
- cellIter != cellsEnd;
- ++cellIter) {
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
// Compute geometry information for current cell
- quadrature->computeGeometry(mesh, coordinates, *cellIter);
+ quadrature->computeGeometry(mesh, coordinates, *c_iter);
const double_array& quadPts = quadrature->quadPts();
@@ -187,18 +123,12 @@
<< "using spatial database '" << _db->label() << "'.";
throw std::runtime_error(msg.str());
} // if
- _dbToParameters(¶mData, queryData);
+ _dbToParameters(&cellData[numParamsQuadPt*iQuadPt], numParamsQuadPt,
+ queryData);
- for (int iParam=0; iParam < numParams; ++iParam) {
- const int numValues = _numParamValues[iParam];
- for (int iValue=0; iValue < numValues; ++iValue)
- cellData[iParam][iQuadPt*numValues+iValue] =
- paramData[iParam][iValue];
- } // for
} // for
// Insert cell contribution into fields
- for (int iParam=0; iParam < numParams; ++iParam)
- paramSections[iParam]->updatePoint(*cellIter, &cellData[iParam][0]);
+ _parameters->updatePoint(*c_iter, &cellData[0]);
} // for
// Close database
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -22,8 +22,9 @@
#if !defined(pylith_materials_material_hh)
#define pylith_materials_material_hh
-#include "pylith/utils/array.hh" // USES std::vector, double_array
+#include "pylith/utils/array.hh" // USES double_array
#include <string> // HASA std::string
+#include "pylith/utils/sievetypes.hh" // USES real_section_type
/// Namespace for pylith package
namespace pylith {
@@ -51,12 +52,6 @@
} // geocoords
} // spatialdata
-/// Namespace for Sieve package.
-namespace ALE {
- class Mesh;
- template<class T> class Obj;
-} // ALE
-
/// C++ abstract base class for Material object.
class pylith::materials::Material
{ // class Material
@@ -167,13 +162,6 @@
virtual
int _numDBValues(void) const = 0;
- /** Get names of parameters for physical properties.
- *
- * @returns Names of parameters
- */
- virtual
- const char** _parameterNames(void) const = 0;
-
/** Get number of values for each parameter for physical properties.
*
* @returns Array of number of values for each parameter.
@@ -182,13 +170,13 @@
/** Compute parameters from values in spatial database.
*
- * @param paramVals Array of parameters
- * @param numParams Number of parameters
+ * @param paramVals Array of parameters.
+ * @param numParams Number of parameters at quadrature point.
* @param dbValues Array of database values
- * @param numValues Number of database values
*/
virtual
- void _dbToParameters(std::vector<double_array>* const paramVals,
+ void _dbToParameters(double* const paramVals,
+ const int numParams,
const double_array& dbValues) const = 0;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
@@ -205,9 +193,10 @@
double _dt; ///< Current time step
- ///< Manager of parameters for physical properties of material
- topology::FieldsManager* _parameters;
-
+ ///< Section containing parameters for physical properties of material.
+ ALE::Obj<real_section_type> _parameters;
+
+ int _numParamsQuadPt; ///< Number of parameters per quadrature point.
int _dimension; ///< Spatial dimension associated with material
bool _needNewJacobian; ///< True if need to reform Jacobian, false otherwise
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -46,16 +46,14 @@
/// Parameters
const int numParameters = 6;
const int numParamValues[] = { 1, 1, 1, 1, 6, 6};
- const char* namesParameters[] =
- {"density", "mu", "lambda" , "maxwellTime", "strainT", "visStrain"};
/// Indices (order) of parameters
const int pidDensity = 0;
- const int pidMu = 1;
- const int pidLambda = 2;
- const int pidMaxwellTime = 3;
- const int pidStrainT = 4;
- const int pidVisStrain = 5;
+ const int pidMu = pidDensity + 1;
+ const int pidLambda = pidMu + 1;
+ const int pidMaxwellTime = pidLambda + 1;
+ const int pidStrainT = pidMaxwellTime + 1;
+ const int pidVisStrain = pidStrainT + 6;
} // _MaxwellIsotropic3D
} // materials
} // pylith
@@ -63,7 +61,7 @@
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::MaxwellIsotropic3D::MaxwellIsotropic3D(void) :
- ElasticMaterial(_MaxwellIsotropic3D::numParamValues,
+ ElasticMaterial(_MaxwellIsotropic3D::numParamValues,
_MaxwellIsotropic3D::numParameters),
_calcElasticConstsFn(&pylith::materials::MaxwellIsotropic3D::_calcElasticConstsElastic),
_calcStressFn(&pylith::materials::MaxwellIsotropic3D::_calcStressElastic),
@@ -95,26 +93,17 @@
} // _numDBValues
// ----------------------------------------------------------------------
-// Get names of parameters for physical properties.
-const char**
-pylith::materials::MaxwellIsotropic3D::_parameterNames(void) const
-{ // _parameterNames
- return _MaxwellIsotropic3D::namesParameters;
-} // _parameterNames
-
-// ----------------------------------------------------------------------
// Compute parameters from values in spatial database.
void
-pylith::materials::MaxwellIsotropic3D::_dbToParameters(std::vector<double_array>* paramVals,
- const double_array& dbValues) const
+pylith::materials::MaxwellIsotropic3D::_dbToParameters(
+ double* const paramVals,
+ const int numParams,
+ const double_array& dbValues) const
{ // _dbToParameters
assert(0 != paramVals);
- const int numParams = paramVals->size();
- assert(_MaxwellIsotropic3D::numParameters == numParams);
+ assert(_numParamsQuadPt == numParams);
const int numDBValues = dbValues.size();
assert(_MaxwellIsotropic3D::numDBValues == numDBValues);
- for (int i=0; i < numParams; ++i)
- assert(_MaxwellIsotropic3D::numParamValues[i] == (*paramVals)[i].size());
const double density = dbValues[_MaxwellIsotropic3D::didDensity];
const double vs = dbValues[_MaxwellIsotropic3D::didVs];
@@ -135,10 +124,10 @@
assert(mu > 0);
const double maxwelltime = viscosity / mu;
- (*paramVals)[_MaxwellIsotropic3D::pidDensity][0] = density;
- (*paramVals)[_MaxwellIsotropic3D::pidMu][0] = mu;
- (*paramVals)[_MaxwellIsotropic3D::pidLambda][0] = lambda;
- (*paramVals)[_MaxwellIsotropic3D::pidMaxwellTime][0] = maxwelltime;
+ paramVals[_MaxwellIsotropic3D::pidDensity] = density;
+ paramVals[_MaxwellIsotropic3D::pidMu] = mu;
+ paramVals[_MaxwellIsotropic3D::pidLambda] = lambda;
+ paramVals[_MaxwellIsotropic3D::pidMaxwellTime] = maxwelltime;
PetscLogFlopsNoCheck(7);
} // _dbToParameters
@@ -162,16 +151,15 @@
// ----------------------------------------------------------------------
// Compute density at location from parameters.
void
-pylith::materials::MaxwellIsotropic3D::_calcDensity(
- double_array* const density,
- const std::vector<double_array>& parameters)
+pylith::materials::MaxwellIsotropic3D::_calcDensity(double* const density,
+ const double* parameters,
+ const int numParams)
{ // _calcDensity
assert(0 != density);
- assert(1 == density->size());
- assert(_MaxwellIsotropic3D::numParameters == parameters.size());
- assert(1 == parameters[_MaxwellIsotropic3D::pidDensity].size());
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
- (*density)[0] = parameters[_MaxwellIsotropic3D::pidDensity][0];
+ density[0] = parameters[_MaxwellIsotropic3D::pidDensity];
} // _calcDensity
// ----------------------------------------------------------------------
@@ -179,24 +167,24 @@
// material.
void
pylith::materials::MaxwellIsotropic3D::_calcStressElastic(
- double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+ double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _calcStressElastic
assert(0 != stress);
- assert(_MaxwellIsotropic3D::tensorSize == stress->size());
- assert(_MaxwellIsotropic3D::numParameters == parameters.size());
- assert(_MaxwellIsotropic3D::tensorSize == totalStrain.size());
- assert(1 == parameters[_MaxwellIsotropic3D::pidDensity].size());
- assert(1 == parameters[_MaxwellIsotropic3D::pidMu].size());
- assert(1 == parameters[_MaxwellIsotropic3D::pidLambda].size());
- assert(6 == parameters[_MaxwellIsotropic3D::pidStrainT].size());
- assert(6 == parameters[_MaxwellIsotropic3D::pidVisStrain].size());
+ assert(_MaxwellIsotropic3D::tensorSize == stressSize);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- const double density = parameters[_MaxwellIsotropic3D::pidDensity][0];
- const double mu = parameters[_MaxwellIsotropic3D::pidMu][0];
- const double lambda = parameters[_MaxwellIsotropic3D::pidLambda][0];
- const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime][0];
+ const double density = parameters[_MaxwellIsotropic3D::pidDensity];
+ const double mu = parameters[_MaxwellIsotropic3D::pidMu];
+ const double lambda = parameters[_MaxwellIsotropic3D::pidLambda];
+ const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime];
const double mu2 = 2.0 * mu;
@@ -210,12 +198,12 @@
const double traceStrainTpdt = e11 + e22 + e33;
const double s123 = lambda * traceStrainTpdt;
- (*stress)[0] = s123 + mu2*e11;
- (*stress)[1] = s123 + mu2*e22;
- (*stress)[2] = s123 + mu2*e33;
- (*stress)[3] = mu2 * e12;
- (*stress)[4] = mu2 * e23;
- (*stress)[5] = mu2 * e13;
+ stress[0] = s123 + mu2*e11;
+ stress[1] = s123 + mu2*e22;
+ stress[2] = s123 + mu2*e33;
+ stress[3] = mu2 * e12;
+ stress[4] = mu2 * e23;
+ stress[5] = mu2 * e13;
// std::cout << " _calcStressElastic: " << std::endl;
// std::cout << " totalStrain: " << std::endl;
// for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp)
@@ -233,24 +221,24 @@
// material.
void
pylith::materials::MaxwellIsotropic3D::_calcStressViscoelastic(
- double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
-{ // _calcStressViscoelastic
+ double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
+{ // _calcStressElastic
assert(0 != stress);
- assert(_MaxwellIsotropic3D::tensorSize == stress->size());
- assert(_MaxwellIsotropic3D::numParameters == parameters.size());
- assert(_MaxwellIsotropic3D::tensorSize == totalStrain.size());
- assert(1 == parameters[_MaxwellIsotropic3D::pidDensity].size());
- assert(1 == parameters[_MaxwellIsotropic3D::pidMu].size());
- assert(1 == parameters[_MaxwellIsotropic3D::pidLambda].size());
- assert(6 == parameters[_MaxwellIsotropic3D::pidStrainT].size());
- assert(6 == parameters[_MaxwellIsotropic3D::pidVisStrain].size());
+ assert(_MaxwellIsotropic3D::tensorSize == stressSize);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- const double density = parameters[_MaxwellIsotropic3D::pidDensity][0];
- const double mu = parameters[_MaxwellIsotropic3D::pidMu][0];
- const double lambda = parameters[_MaxwellIsotropic3D::pidLambda][0];
- const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime][0];
+ const double density = parameters[_MaxwellIsotropic3D::pidDensity];
+ const double mu = parameters[_MaxwellIsotropic3D::pidMu];
+ const double lambda = parameters[_MaxwellIsotropic3D::pidLambda];
+ const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime];
const double mu2 = 2.0 * mu;
const double bulkmodulus = lambda + mu2/3.0;
@@ -275,9 +263,9 @@
// << " " << parameters[_MaxwellIsotropic3D::pidVisStrain][iComp]
// << std::endl;
- const double meanStrainT = (parameters[_MaxwellIsotropic3D::pidStrainT][0] +
- parameters[_MaxwellIsotropic3D::pidStrainT][1] +
- parameters[_MaxwellIsotropic3D::pidStrainT][2])/3.0;
+ const double meanStrainT = (parameters[_MaxwellIsotropic3D::pidStrainT+0] +
+ parameters[_MaxwellIsotropic3D::pidStrainT+1] +
+ parameters[_MaxwellIsotropic3D::pidStrainT+2])/3.0;
PetscLogFlopsNoCheck(11);
// The code below should probably be in a separate function since it
@@ -314,17 +302,18 @@
// std::cout << " stress totalStrain visStrain: " << std::endl;
for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp) {
devStrainTpdt = totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
- devStrainT = parameters[_MaxwellIsotropic3D::pidStrainT][iComp] -
+ devStrainT = parameters[_MaxwellIsotropic3D::pidStrainT+iComp] -
diag[iComp]*meanStrainT;
- visStrain = expFac*parameters[_MaxwellIsotropic3D::pidVisStrain][iComp] +
+ visStrain = expFac*parameters[_MaxwellIsotropic3D::pidVisStrain+iComp] +
dq*(devStrainTpdt - devStrainT);
devStressTpdt = elasFac*visStrain;
// Later I will want to put in initial stresses.
- (*stress)[iComp] =diag[iComp]*meanStressTpdt+devStressTpdt;
+ stress[iComp] =diag[iComp]*meanStressTpdt+devStressTpdt;
// Temporary to get stresses and strains.
- // std::cout << " " << (*stress)[iComp] << " " << totalStrain[iComp] << " " << visStrain << std:: endl;
+ // std::cout << " " << stress[iComp] << " " << totalStrain[iComp] << " " << visStrain << std:: endl;
} // for
+
PetscLogFlopsNoCheck(11 * _MaxwellIsotropic3D::tensorSize);
} // _calcStress
@@ -332,45 +321,50 @@
// Compute derivative of elasticity matrix at location from parameters.
void
pylith::materials::MaxwellIsotropic3D::_calcElasticConstsElastic(
- double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _calcElasticConstsElastic
assert(0 != elasticConsts);
- assert(_MaxwellIsotropic3D::numElasticConsts == elasticConsts->size());
- assert(_MaxwellIsotropic3D::numParameters == parameters.size());
- assert(_MaxwellIsotropic3D::tensorSize == totalStrain.size());
+ assert(_MaxwellIsotropic3D::numElasticConsts == numElasticConsts);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- const double density = parameters[_MaxwellIsotropic3D::pidDensity][0];
- const double mu = parameters[_MaxwellIsotropic3D::pidMu][0];
- const double lambda = parameters[_MaxwellIsotropic3D::pidLambda][0];
- const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime][0];
+ const double density = parameters[_MaxwellIsotropic3D::pidDensity];
+ const double mu = parameters[_MaxwellIsotropic3D::pidMu];
+ const double lambda = parameters[_MaxwellIsotropic3D::pidLambda];
+ const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime];
const double mu2 = 2.0 * mu;
const double lambda2mu = lambda + mu2;
- const double bulkmodulus = lambda + mu2/3.0;
+ const double bulkmodulus = lambda + mu2 / 3.0;
- (*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] = mu2; // C1212
- (*elasticConsts)[16] = 0; // C1223
- (*elasticConsts)[17] = 0; // C1213
- (*elasticConsts)[18] = mu2; // C2323
- (*elasticConsts)[19] = 0; // C2313
- (*elasticConsts)[20] = mu2; // 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] = mu2; // C1212
+ elasticConsts[16] = 0; // C1223
+ elasticConsts[17] = 0; // C1213
+ elasticConsts[18] = mu2; // C2323
+ elasticConsts[19] = 0; // C2313
+ elasticConsts[20] = mu2; // C1313
PetscLogFlopsNoCheck(4);
} // _calcElasticConstsElastic
@@ -380,19 +374,24 @@
// as an elastic material.
void
pylith::materials::MaxwellIsotropic3D::_calcElasticConstsViscoelastic(
- double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain)
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _calcElasticConstsViscoelastic
assert(0 != elasticConsts);
- assert(_MaxwellIsotropic3D::numElasticConsts == elasticConsts->size());
- assert(_MaxwellIsotropic3D::numParameters == parameters.size());
- assert(_MaxwellIsotropic3D::tensorSize == totalStrain.size());
+ assert(_MaxwellIsotropic3D::numElasticConsts == numElasticConsts);
+ assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- const double density = parameters[_MaxwellIsotropic3D::pidDensity][0];
- const double mu = parameters[_MaxwellIsotropic3D::pidMu][0];
- const double lambda = parameters[_MaxwellIsotropic3D::pidLambda][0];
- const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime][0];
+ const double density = parameters[_MaxwellIsotropic3D::pidDensity];
+ const double mu = parameters[_MaxwellIsotropic3D::pidMu];
+ const double lambda = parameters[_MaxwellIsotropic3D::pidLambda];
+ const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime];
const double mu2 = 2.0 * mu;
const double bulkmodulus = lambda + mu2/3.0;
@@ -420,27 +419,27 @@
} // else
const double visFac = mu*dq/3.0;
- (*elasticConsts)[ 0] = bulkmodulus + 4.0*visFac; // C1111
- (*elasticConsts)[ 1] = bulkmodulus - 2.0*visFac; // C1122
- (*elasticConsts)[ 2] = (*elasticConsts)[1]; // C1133
- (*elasticConsts)[ 3] = 0; // C1112
- (*elasticConsts)[ 4] = 0; // C1123
- (*elasticConsts)[ 5] = 0; // C1113
- (*elasticConsts)[ 6] = (*elasticConsts)[0]; // C2222
- (*elasticConsts)[ 7] = (*elasticConsts)[1]; // C2233
- (*elasticConsts)[ 8] = 0; // C2212
- (*elasticConsts)[ 9] = 0; // C2223
- (*elasticConsts)[10] = 0; // C2213
- (*elasticConsts)[11] = (*elasticConsts)[0]; // C3333
- (*elasticConsts)[12] = 0; // C3312
- (*elasticConsts)[13] = 0; // C3323
- (*elasticConsts)[14] = 0; // C3313
- (*elasticConsts)[15] = 6.0 * visFac; // C1212
- (*elasticConsts)[16] = 0; // C1223
- (*elasticConsts)[17] = 0; // C1213
- (*elasticConsts)[18] = (*elasticConsts)[15]; // C2323
- (*elasticConsts)[19] = 0; // C2313
- (*elasticConsts)[20] = (*elasticConsts)[15]; // C1313
+ elasticConsts[ 0] = bulkmodulus + 4.0*visFac; // C1111
+ elasticConsts[ 1] = bulkmodulus - 2.0*visFac; // C1122
+ elasticConsts[ 2] = elasticConsts[1]; // C1133
+ elasticConsts[ 3] = 0; // C1112
+ elasticConsts[ 4] = 0; // C1123
+ elasticConsts[ 5] = 0; // C1113
+ elasticConsts[ 6] = elasticConsts[0]; // C2222
+ elasticConsts[ 7] = elasticConsts[1]; // C2233
+ elasticConsts[ 8] = 0; // C2212
+ elasticConsts[ 9] = 0; // C2223
+ elasticConsts[10] = 0; // C2213
+ elasticConsts[11] = elasticConsts[0]; // C3333
+ elasticConsts[12] = 0; // C3312
+ elasticConsts[13] = 0; // C3323
+ elasticConsts[14] = 0; // C3313
+ elasticConsts[15] = 6.0 * visFac; // C1212
+ elasticConsts[16] = 0; // C1223
+ elasticConsts[17] = 0; // C1213
+ elasticConsts[18] = elasticConsts[15]; // C2323
+ elasticConsts[19] = 0; // C2313
+ elasticConsts[20] = elasticConsts[15]; // C1313
PetscLogFlopsNoCheck(7);
} // _calcElasticConstsViscoelastic
@@ -449,21 +448,18 @@
// Update state variables.
void
pylith::materials::MaxwellIsotropic3D::_updateStateElastic(
- std::vector<double_array>* parameters,
- const double_array& totalStrain)
+ double* const parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _updateStateElastic
assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- assert(_MaxwellIsotropic3D::numParameters == parameters->size());
- assert(_MaxwellIsotropic3D::tensorSize == totalStrain.size());
- assert(1 == (*parameters)[_MaxwellIsotropic3D::pidDensity].size());
- assert(1 == (*parameters)[_MaxwellIsotropic3D::pidMu].size());
- assert(1 == (*parameters)[_MaxwellIsotropic3D::pidLambda].size());
- assert(6 == (*parameters)[_MaxwellIsotropic3D::pidStrainT].size());
- assert(6 == (*parameters)[_MaxwellIsotropic3D::pidVisStrain].size());
+ const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime];
- const double maxwelltime = (*parameters)[_MaxwellIsotropic3D::pidMaxwellTime][0];
-
const double e11 = totalStrain[0];
const double e22 = totalStrain[1];
const double e33 = totalStrain[2];
@@ -478,8 +474,8 @@
// _calcStressElastic(&stress, (*parameters), totalStrain);
for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp) {
- (*parameters)[_MaxwellIsotropic3D::pidStrainT][iComp] = totalStrain[iComp];
- (*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp] =
+ parameters[_MaxwellIsotropic3D::pidStrainT+iComp] = totalStrain[iComp];
+ parameters[_MaxwellIsotropic3D::pidVisStrain+iComp] =
totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
} // for
PetscLogFlopsNoCheck(5 * _MaxwellIsotropic3D::tensorSize);
@@ -487,59 +483,59 @@
// std::cout << " updateStateElastic: "<< std::endl;
// std::cout << " StrainT VisStrain Stress: " << std::endl;
// for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp)
-// std::cout << " " << (*parameters)[_MaxwellIsotropic3D::pidStrainT][iComp]
-// << " " << (*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp]
+// std::cout << " " << parameters[_MaxwellIsotropic3D::pidStrainT+iComp]
+// << " " << parameters[_MaxwellIsotropic3D::pidVisStrain+iComp]
// << " " << stress[iComp]
// << std::endl;
+
_needNewJacobian = true;
} // _updateStateElastic
-// ************** Finish adding PETSc logging from here *******************
-
// ----------------------------------------------------------------------
// Update state variables.
void
pylith::materials::MaxwellIsotropic3D::_updateStateViscoelastic(
- std::vector<double_array>* parameters,
- const double_array& totalStrain)
+ double* const parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize)
{ // _updateStateViscoelastic
assert(0 != parameters);
+ assert(_numParamsQuadPt == numParams);
+ assert(0 != totalStrain);
+ assert(_MaxwellIsotropic3D::tensorSize == strainSize);
- assert(_MaxwellIsotropic3D::numParameters == parameters->size());
- assert(_MaxwellIsotropic3D::tensorSize == totalStrain.size());
- assert(1 == (*parameters)[_MaxwellIsotropic3D::pidDensity].size());
- assert(1 == (*parameters)[_MaxwellIsotropic3D::pidMu].size());
- assert(1 == (*parameters)[_MaxwellIsotropic3D::pidLambda].size());
- assert(6 == (*parameters)[_MaxwellIsotropic3D::pidStrainT].size());
- assert(6 == (*parameters)[_MaxwellIsotropic3D::pidVisStrain].size());
+ const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime];
- const double maxwelltime = (*parameters)[_MaxwellIsotropic3D::pidMaxwellTime][0];
-
const double e11 = totalStrain[0];
const double e22 = totalStrain[1];
const double e33 = totalStrain[2];
- const double meanStrainTpdt = (e11 + e22 + e33)/3.0;
+ const double meanStrainTpdt = (e11 + e22 + e33) / 3.0;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
// Temporary to get stresses.
- double_array stress(6);
- _calcStressViscoelastic(&stress, (*parameters), totalStrain);
+ double stress[6];
+ const int stressSize = 6;
+ _calcStressViscoelastic(stress, stressSize,
+ parameters, numParams,
+ totalStrain, strainSize);
const double meanStrainT =
- ((*parameters)[_MaxwellIsotropic3D::pidStrainT][0] +
- (*parameters)[_MaxwellIsotropic3D::pidStrainT][1] +
- (*parameters)[_MaxwellIsotropic3D::pidStrainT][2])/3.0;
+ (parameters[_MaxwellIsotropic3D::pidStrainT+0] +
+ parameters[_MaxwellIsotropic3D::pidStrainT+1] +
+ parameters[_MaxwellIsotropic3D::pidStrainT+2]) / 3.0;
PetscLogFlopsNoCheck(6);
+
// The code below should probably be in a separate function since it
// is used more than once. I should also probably cover the possibility
// that Maxwell time is zero (although this should never happen).
const double timeFrac = 1.0e-5;
const int numTerms = 5;
double dq = 0.0;
- if(maxwelltime < timeFrac*_dt) {
+ if (maxwelltime < timeFrac*_dt) {
double fSign = 1.0;
double factorial = 1.0;
double fraction = 1.0;
@@ -563,13 +559,13 @@
PetscLogFlopsNoCheck(3);
for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp) {
devStrainTpdt = totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
- devStrainT = (*parameters)[_MaxwellIsotropic3D::pidStrainT][iComp] -
+ devStrainT = parameters[_MaxwellIsotropic3D::pidStrainT+iComp] -
diag[iComp] * meanStrainT;
visStrain = expFac *
- (*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp] +
+ parameters[_MaxwellIsotropic3D::pidVisStrain+iComp] +
dq * (devStrainTpdt - devStrainT);
- (*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp] = visStrain;
- (*parameters)[_MaxwellIsotropic3D::pidStrainT][iComp] = totalStrain[iComp];
+ parameters[_MaxwellIsotropic3D::pidVisStrain+iComp] = visStrain;
+ parameters[_MaxwellIsotropic3D::pidStrainT+iComp] = totalStrain[iComp];
} // for
PetscLogFlopsNoCheck(8 * _MaxwellIsotropic3D::tensorSize);
@@ -579,8 +575,8 @@
// std::cout << " updateStateViscoelastic: "<< std::endl;
// std::cout << " StrainT VisStrain Stress: " << std::endl;
// for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp)
-// std::cout << " " << (*parameters)[_MaxwellIsotropic3D::pidStrainT][iComp]
-// << " " << (*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp]
+// std::cout << " " << parameters[_MaxwellIsotropic3D::pidStrainT+iComp]
+// << " " << parameters[_MaxwellIsotropic3D::pidVisStrain+iComp]
// << " " << stress[iComp]
// << std::endl;
} // _updateStateViscoelastic
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -85,12 +85,6 @@
*/
int _numDBValues(void) const;
- /** Get names of parameters for physical properties.
- *
- * @returns Names of parameters
- */
- const char** _parameterNames(void) const;
-
/** Compute parameters from values in spatial database.
*
* Order of values in arrays matches order used in dbValues() and
@@ -99,7 +93,8 @@
* @param paramVals Array of parameters
* @param dbValues Array of database values
*/
- void _dbToParameters(std::vector<double_array>* paramVals,
+ void _dbToParameters(double* const paramVals,
+ const int numParams,
const double_array& dbValues) const;
/** Get number of entries in stress/strain tensors.
@@ -124,120 +119,178 @@
/** Compute density from parameters.
*
- * @param density Array for density
- * @param parameters Parameters at location
+ * @param density Array for density.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
*/
- void _calcDensity(double_array* const density,
- const std::vector<double_array>& parameters);
+ void _calcDensity(double* const density,
+ const double* parameters,
+ const int numParams);
/** Compute stress tensor from parameters.
*
- * @param stress Array for stress tensor
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcStress(double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
/** Compute derivatives of elasticity matrix from parameters.
*
- * @param elasticConsts Array for elastic constants
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcElasticConsts(double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcElasticConsts(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
- /** Update state variables after solve.
+ /** Update parameters (for next time step).
*
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _updateState(std::vector<double_array>* parameters,
- const double_array& totalStrain);
+ void _updateState(double* const parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
// PRIVATE TYPEDEFS ///////////////////////////////////////////////////
private :
/// Member prototype for _calcStress()
typedef void (pylith::materials::MaxwellIsotropic3D::*calcStress_fn_type)
- (double_array* const,
- const std::vector<double_array>&,
- const double_array&);
+ (double* const,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int);
/// Member prototype for _calcElasticConsts()
typedef void (pylith::materials::MaxwellIsotropic3D::*calcElasticConsts_fn_type)
- (double_array* const,
- const std::vector<double_array>&,
- const double_array&);
+ (double* const,
+ const int,
+ const double*,
+ const int,
+ const double*,
+ const int);
/// Member prototype for _updateState()
typedef void (pylith::materials::MaxwellIsotropic3D::*updateState_fn_type)
- (std::vector<double_array>* const,
- const double_array&);
+ (double* const,
+ const int,
+ const double*,
+ const int);
// PRIVATE METHODS ////////////////////////////////////////////////////
private :
/** Compute stress tensor from parameters as an elastic material.
*
- * @param stress Array for stress tensor
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
* @param parameters Parameters at locations.
+ * @param numParams Number of parameters.
* @param totalStrain Total strain at locations.
+ * @param strainSize Size of strain tensor.
*/
- void _calcStressElastic(double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcStressElastic(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
/** Compute stress tensor from parameters as an viscoelastic material.
*
- * @param stress Array for stress tensor
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
* @param parameters Parameters at locations.
+ * @param numParams Number of parameters.
* @param totalStrain Total strain at locations.
+ * @param strainSize Size of strain tensor.
*/
- void _calcStressViscoelastic(double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcStressViscoelastic(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
/** Compute derivatives of elasticity matrix from parameters as an
* elastic material.
*
- * @param elasticConsts Array for elastic constants
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcElasticConstsElastic(double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcElasticConstsElastic(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
/** Compute derivatives of elasticity matrix from parameters as a
* viscoelastic material.
*
- * @param elasticConsts Array for elastic constants
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _calcElasticConstsViscoelastic(double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain);
+ void _calcElasticConstsViscoelastic(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
/** Update state variables after solve as an elastic material.
*
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _updateStateElastic(std::vector<double_array>* parameters,
- const double_array& totalStrain);
+ void _updateStateElastic(double* const parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
/** Update state variables after solve as a viscoelastic material.
*
- * @param parameters Parameters at locations.
- * @param totalStrain Total strain at locations.
+ * @param parameters Parameters at location.
+ * @param numParams Number of parameters.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
*/
- void _updateStateViscoelastic(std::vector<double_array>* parameters,
- const double_array& totalStrain);
+ void _updateStateViscoelastic(double* const parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.icc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -58,34 +58,44 @@
// Compute stress tensor from parameters.
inline
void
-pylith::materials::MaxwellIsotropic3D::_calcStress(
- double_array* const stress,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain) {
+pylith::materials::MaxwellIsotropic3D::_calcStress(double* const stress,
+ const int stressSize,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize) {
assert(0 != _calcStressFn);
- CALL_MEMBER_FN(*this, _calcStressFn)(stress, parameters, totalStrain);
+ CALL_MEMBER_FN(*this, _calcStressFn)(stress, stressSize,
+ parameters, numParams,
+ totalStrain, strainSize);
} // _calcStress
// Compute derivatives of elasticity matrix from parameters.
inline
void
pylith::materials::MaxwellIsotropic3D::_calcElasticConsts(
- double_array* const elasticConsts,
- const std::vector<double_array>& parameters,
- const double_array& totalStrain) {
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize) {
assert(0 != _calcElasticConstsFn);
- CALL_MEMBER_FN(*this, _calcElasticConstsFn)
- (elasticConsts, parameters, totalStrain);
+ CALL_MEMBER_FN(*this, _calcElasticConstsFn)(elasticConsts, numElasticConsts,
+ parameters, numParams,
+ totalStrain, strainSize);
} // _calcElasticConsts
// Update state variables after solve.
inline
void
-pylith::materials::MaxwellIsotropic3D::_updateState(
- std::vector<double_array>* parameters,
- const double_array& totalStrain) {
+pylith::materials::MaxwellIsotropic3D::_updateState(double* const parameters,
+ const int numParams,
+ const double* totalStrain,
+ const int strainSize) {
assert(0 != _updateStateFn);
- CALL_MEMBER_FN(*this, _updateStateFn)(parameters, totalStrain);
+ CALL_MEMBER_FN(*this, _updateStateFn)(parameters, numParams,
+ totalStrain, strainSize);
} // _updateState
// End of file
Modified: short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -123,7 +123,12 @@
} else if (0 != fields) {
// If field is not in mesh, try to get it from fields manager
// ADD STUFF HERE
- } // else
+ } else {
+ std::ostringstream msg;
+ msg << "Could not find field '" << _vertexFields[i]
+ << "' in fields available for output.";
+ throw std::runtime_error(msg.str());
+ } // if/else
// Extract values from section
// ADD STUFF HERE
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestIntegratorElasticity.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -44,24 +44,19 @@
const double dispVals[] = { 0.5, 1.5 };
const double strainE[] = { 0.5, 0.25 };
- std::vector<double_array> strain(numQuadPts);
- for (int i=0; i < numQuadPts; ++i)
- strain[i].resize(tensorSize);
+ const int size = numQuadPts * tensorSize;
+ double_array strain(size);
double_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
double_array disp(dispVals, numBasis*dim);
- IntegratorElasticity::_calcTotalStrain1D(&strain,
- basisDeriv, disp, numBasis);
+ IntegratorElasticity::_calcTotalStrain1D(&strain,
+ basisDeriv, disp, numBasis, numQuadPts);
const double tolerance = 1.0e-06;
- CPPUNIT_ASSERT_EQUAL(numQuadPts, int(strain.size()));
- for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
- CPPUNIT_ASSERT_EQUAL(tensorSize, int(strain[iQuad].size()));
- for (int iStrain=0; iStrain < tensorSize; ++iStrain)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i++], strain[iQuad][iStrain],
- tolerance);
- } // for
+ CPPUNIT_ASSERT_EQUAL(size, int(strain.size()));
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i], strain[i], tolerance);
} // testCalcTotalStrain1D
// ----------------------------------------------------------------------
@@ -97,24 +92,19 @@
0.6, -0.4, 1.3
};
- std::vector<double_array> strain(numQuadPts);
- for (int i=0; i < numQuadPts; ++i)
- strain[i].resize(tensorSize);
+ const int size = numQuadPts * tensorSize;
+ double_array strain(size);
double_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
double_array disp(dispVals, numBasis*dim);
- IntegratorElasticity::_calcTotalStrain2D(&strain,
- basisDeriv, disp, numBasis);
+ IntegratorElasticity::_calcTotalStrain2D(&strain, basisDeriv, disp,
+ numBasis, numQuadPts);
const double tolerance = 1.0e-06;
- CPPUNIT_ASSERT_EQUAL(numQuadPts, int(strain.size()));
- for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
- CPPUNIT_ASSERT_EQUAL(tensorSize, int(strain[iQuad].size()));
- for (int iStrain=0; iStrain < tensorSize; ++iStrain)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i++], strain[iQuad][iStrain],
- tolerance);
- } // for
+ CPPUNIT_ASSERT_EQUAL(size, int(strain.size()));
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i], strain[i], tolerance);
} // testCalcTotalStrain2D
// ----------------------------------------------------------------------
@@ -160,24 +150,19 @@
0.9, -0.6, -0.9, 1.95, 0.75, 0.9
};
- std::vector<double_array> strain(numQuadPts);
- for (int i=0; i < numQuadPts; ++i)
- strain[i].resize(tensorSize);
+ const int size = numQuadPts * tensorSize;
+ double_array strain(size);
double_array basisDeriv(basisDerivVals, numQuadPts*numBasis*dim);
double_array disp(dispVals, numBasis*dim);
- IntegratorElasticity::_calcTotalStrain3D(&strain,
- basisDeriv, disp, numBasis);
+ IntegratorElasticity::_calcTotalStrain3D(&strain, basisDeriv, disp,
+ numBasis, numQuadPts);
const double tolerance = 1.0e-06;
- CPPUNIT_ASSERT_EQUAL(numQuadPts, int(strain.size()));
- for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
- CPPUNIT_ASSERT_EQUAL(tensorSize, int(strain[iQuad].size()));
- for (int iStrain=0; iStrain < tensorSize; ++iStrain)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i++], strain[iQuad][iStrain],
- tolerance);
- } // for
+ CPPUNIT_ASSERT_EQUAL(size, int(strain.size()));
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(strainE[i], strain[i], tolerance);
} // testCalcTotalStrain3D
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticIsotropic3D.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticIsotropic3D.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -100,12 +100,11 @@
const int numParams = data.numParameters;
- std::vector<double_array> parameters(numParams);
+ double_array parameters(numParams);
const int paramsSize = 1;
for (int i=0; i < numParams; ++i) {
- parameters[i].resize(numParams);
for (int j=0; j < paramsSize; ++j)
- parameters[i][j] = i+j;
+ parameters[i*paramsSize+j] = i+j;
} // for
const int tensorSize = 9;
@@ -113,12 +112,13 @@
for (int i=0; i < tensorSize; ++i)
totalStrain[i] = i;
- material._updateState(¶meters, totalStrain);
+ material._updateState(¶meters[0], numParams, &totalStrain[0], tensorSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParams; ++i)
for (int j=0; j < paramsSize; ++j)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i+j), parameters[i][j], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i+j), parameters[i*paramsSize+j],
+ tolerance);
for (int i=0; i < tensorSize; ++i)
CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i), totalStrain[i], tolerance);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -59,50 +59,29 @@
ElasticIsotropic3D material;
ElasticIsotropic3DData data;
- delete material._parameters;
- material._parameters = new topology::FieldsManager(mesh);
const int numQuadPts = 2;
- const int fiberDim = numQuadPts; // number of values in field per cell
+ const int_array& numParamValues = material._getNumParamValues();
+ const int numParams = data.numParameters;
+ const int numParamsQuadPt = data.numParamsQuadPt;
+
+ Mesh::label_sequence::iterator c_iter = cells->begin();
- Mesh::label_sequence::iterator cellIter=cells->begin();
+ const int fiberDim = numQuadPts * numParamsQuadPt;
+ material._parameters = new real_section_type(mesh->comm(), mesh->debug());
+ material._parameters->setFiberDimension(cells, fiberDim);
+ mesh->allocate(material._parameters);
- material._parameters->addReal("density");
- const ALE::Obj<real_section_type>& parameterDensity =
- material._parameters->getReal("density");
- parameterDensity->setFiberDimension(cells, fiberDim);
- mesh->allocate(parameterDensity);
- double cellData[numQuadPts];
- cellData[0] = data.parameterData[0];
- cellData[1] = data.parameterData[3];
- parameterDensity->updateAddPoint(*cellIter, cellData);
+ material._parameters->updatePoint(*c_iter, data.parameterData);
- material._parameters->addReal("mu");
- const ALE::Obj<real_section_type>& parameterMu =
- material._parameters->getReal("mu");
- parameterMu->setFiberDimension(cells, fiberDim);
- mesh->allocate(parameterMu);
- cellData[0] = data.parameterData[1];
- cellData[1] = data.parameterData[4];
- parameterMu->updateAddPoint(*cellIter, cellData);
+ material.getStateVarsCell(*c_iter, numQuadPts);
+ const double_array& density = material.calcDensity();
- material._parameters->addReal("lambda");
- const ALE::Obj<real_section_type>& parameterLambda =
- material._parameters->getReal("lambda");
- parameterLambda->setFiberDimension(cells, fiberDim);
- mesh->allocate(parameterLambda);
- cellData[0] = data.parameterData[2];
- cellData[1] = data.parameterData[5];
- parameterLambda->updateAddPoint(*cellIter, cellData);
-
- material.getStateVarsCell(*cellIter, numQuadPts);
- const std::vector<double_array>& density = material.calcDensity();
-
const double tolerance = 1.0e-06;
const double* densityE = data.density;
CPPUNIT_ASSERT(0 != densityE);
const double size = numQuadPts;
for (int i=0; i < size; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, density[i][0]/densityE[i], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, density[i]/densityE[i], tolerance);
} // testCalcDensity
// ----------------------------------------------------------------------
@@ -140,41 +119,20 @@
ElasticIsotropic3D material;
ElasticIsotropic3DData data;
- delete material._parameters;
- material._parameters = new topology::FieldsManager(mesh);
const int numQuadPts = 2;
- const int fiberDim = numQuadPts; // number of values in field per cell
+ const int_array& numParamValues = material._getNumParamValues();
+ const int numParams = data.numParameters;
+ const int numParamsQuadPt = data.numParamsQuadPt;
- Mesh::label_sequence::iterator cellIter=cells->begin();
+ Mesh::label_sequence::iterator c_iter = cells->begin();
- material._parameters->addReal("density");
- const ALE::Obj<real_section_type>& parameterDensity =
- material._parameters->getReal("density");
- parameterDensity->setFiberDimension(cells, fiberDim);
- mesh->allocate(parameterDensity);
- double cellData[numQuadPts];
- cellData[0] = data.parameterData[0];
- cellData[1] = data.parameterData[3];
- parameterDensity->updateAddPoint(*cellIter, cellData);
+ const int fiberDim = numQuadPts * numParamsQuadPt;
+ material._parameters = new real_section_type(mesh->comm(), mesh->debug());
+ material._parameters->setFiberDimension(cells, fiberDim);
+ mesh->allocate(material._parameters);
- material._parameters->addReal("mu");
- const ALE::Obj<real_section_type>& parameterMu =
- material._parameters->getReal("mu");
- parameterMu->setFiberDimension(cells, fiberDim);
- mesh->allocate(parameterMu);
- cellData[0] = data.parameterData[1];
- cellData[1] = data.parameterData[4];
- parameterMu->updateAddPoint(*cellIter, cellData);
+ material._parameters->updatePoint(*c_iter, data.parameterData);
- material._parameters->addReal("lambda");
- const ALE::Obj<real_section_type>& parameterLambda =
- material._parameters->getReal("lambda");
- parameterLambda->setFiberDimension(cells, fiberDim);
- mesh->allocate(parameterLambda);
- cellData[0] = data.parameterData[2];
- cellData[1] = data.parameterData[5];
- parameterLambda->updateAddPoint(*cellIter, cellData);
-
int strainSize = 0;
switch(data.dimension)
{ // switch
@@ -188,31 +146,22 @@
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
+ double_array strain(data.strain, numQuadPts*strainSize);
- material.getStateVarsCell(*cellIter, numQuadPts);
- const std::vector<double_array>& stress = material.calcStress(strain);
+ material.getStateVarsCell(*c_iter, numQuadPts);
+ const double_array& stress = material.calcStress(strain);
const double tolerance = 1.0e-06;
const double* stressE = data.stress;
CPPUNIT_ASSERT(0 != stressE);
- 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);
+ const int size = numQuadPts * strainSize;;
+ CPPUNIT_ASSERT_EQUAL(size, int(stress.size()));
+ for (int i=0; i < size; ++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[iQuad][iStress],
- tolerance);
- } // for
- } // for
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(stressE[i], stress[i], tolerance);
} // testCalcStress
// ----------------------------------------------------------------------
@@ -250,41 +199,20 @@
ElasticIsotropic3D material;
ElasticIsotropic3DData data;
- delete material._parameters;
- material._parameters = new topology::FieldsManager(mesh);
const int numQuadPts = 2;
- const int fiberDim = numQuadPts; // number of values in field per cell
+ const int_array& numParamValues = material._getNumParamValues();
+ const int numParams = data.numParameters;
+ const int numParamsQuadPt = data.numParamsQuadPt;
- Mesh::label_sequence::iterator cellIter=cells->begin();
+ Mesh::label_sequence::iterator c_iter = cells->begin();
- material._parameters->addReal("density");
- const ALE::Obj<real_section_type>& parameterDensity =
- material._parameters->getReal("density");
- parameterDensity->setFiberDimension(cells, fiberDim);
- mesh->allocate(parameterDensity);
- double cellData[numQuadPts];
- cellData[0] = data.parameterData[0];
- cellData[1] = data.parameterData[3];
- parameterDensity->updateAddPoint(*cellIter, cellData);
+ const int fiberDim = numQuadPts * numParamsQuadPt;
+ material._parameters = new real_section_type(mesh->comm(), mesh->debug());
+ material._parameters->setFiberDimension(cells, fiberDim);
+ mesh->allocate(material._parameters);
- material._parameters->addReal("mu");
- const ALE::Obj<real_section_type>& parameterMu =
- material._parameters->getReal("mu");
- parameterMu->setFiberDimension(cells, fiberDim);
- mesh->allocate(parameterMu);
- cellData[0] = data.parameterData[1];
- cellData[1] = data.parameterData[4];
- parameterMu->updateAddPoint(*cellIter, cellData);
+ material._parameters->updatePoint(*c_iter, data.parameterData);
- material._parameters->addReal("lambda");
- const ALE::Obj<real_section_type>& parameterLambda =
- material._parameters->getReal("lambda");
- parameterLambda->setFiberDimension(cells, fiberDim);
- mesh->allocate(parameterLambda);
- cellData[0] = data.parameterData[2];
- cellData[1] = data.parameterData[5];
- parameterLambda->updateAddPoint(*cellIter, cellData);
-
int strainSize = 0;
switch(data.dimension)
{ // switch
@@ -298,34 +226,25 @@
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
+ double_array strain(data.strain, numQuadPts*strainSize);
- material.getStateVarsCell(*cellIter, numQuadPts);
- const std::vector<double_array>& elasticConsts =
- material.calcDerivElastic(strain);
+ material.getStateVarsCell(*c_iter, numQuadPts);
+ const double_array& elasticConsts = material.calcDerivElastic(strain);
const double tolerance = 1.0e-06;
const double* elasticConstsE = data.elasticConsts;
CPPUNIT_ASSERT(0 != elasticConstsE);
- for (int iQuad=0, i=0; iQuad < numQuadPts; ++iQuad) {
- const int numElasticConsts = elasticConsts[iQuad].size();
- for (int iConst=0; iConst < numElasticConsts; ++iConst, ++i) {
+ const int size = elasticConsts.size();
+ for (int i=0; i < size; ++i)
if (fabs(elasticConstsE[i]) > tolerance)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,
- elasticConsts[iQuad][iConst]/elasticConstsE[i],
+ elasticConsts[i]/elasticConstsE[i],
tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(elasticConstsE[i],
- elasticConsts[iQuad][iConst],
+ elasticConsts[i],
tolerance);
- } // for
- } // for
} // testCalcDerivElastic
// ----------------------------------------------------------------------
@@ -335,19 +254,10 @@
ElasticMaterial* material,
const ElasticMaterialData& data) const
{ // _testCalcDensity
- const int numParameters = data.numParameters;
- std::vector<double_array> parameters(numParameters);
- for (int iParam=0, i=0; iParam < numParameters; ++iParam) {
- const int numValues = data.numParamValues[iParam];
- parameters[iParam].resize(numValues);
- for (int iValue=0; iValue < numValues; ++iValue)
- parameters[iParam][iValue] = data.parameterData[i++];
- } // for
-
double_array density(1);
- material->_calcDensity(&density, parameters);
+ material->_calcDensity(&density[0], data.parameterData, data.numParamsQuadPt);
+
const double* densityE = data.density;
-
CPPUNIT_ASSERT(0 != densityE);
const double tolerance = 1.0e-06;
@@ -361,15 +271,6 @@
ElasticMaterial* material,
const ElasticMaterialData& data) const
{ // _testCalcElasticConsts
- const int numParameters = data.numParameters;
- std::vector<double_array> parameters(numParameters);
- for (int iParam=0, i=0; iParam < numParameters; ++iParam) {
- const int numValues = data.numParamValues[iParam];
- parameters[iParam].resize(numValues);
- for (int iValue=0; iValue < numValues; ++iValue)
- parameters[iParam][iValue] = data.parameterData[i++];
- } // for
-
int stressSize = 0;
switch(data.dimension)
{ // switch
@@ -384,12 +285,11 @@
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);
+ material->_calcStress(&stress[0], stress.size(),
+ data.parameterData, data.numParamsQuadPt,
+ data.strain, stressSize);
+
const double* stressE = data.stress;
-
CPPUNIT_ASSERT(0 != stressE);
const double tolerance = 1.0e-06;
@@ -409,15 +309,6 @@
ElasticMaterial* material,
const ElasticMaterialData& data) const
{ // _testCalcElasticConsts
- const int numParameters = data.numParameters;
- std::vector<double_array> parameters(numParameters);
- for (int iParam=0, i=0; iParam < numParameters; ++iParam) {
- const int numValues = data.numParamValues[iParam];
- parameters[iParam].resize(numValues);
- for (int iValue=0; iValue < numValues; ++iValue)
- parameters[iParam][iValue] = data.parameterData[i++];
- } // for
-
int numConsts = 0;
int strainSize = 0;
switch(data.dimension)
@@ -436,12 +327,11 @@
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);
+ material->_calcElasticConsts(&elasticConsts[0], numConsts,
+ data.parameterData, data.numParamsQuadPt,
+ data.strain, strainSize);
+
const double* elasticConstsE = data.elasticConsts;
-
CPPUNIT_ASSERT(0 != elasticConstsE);
const double tolerance = 1.0e-06;
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStrain.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStrain.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -100,12 +100,11 @@
const int numParams = data.numParameters;
- std::vector<double_array> parameters(numParams);
+ double_array parameters(numParams);
const int paramsSize = 1;
for (int i=0; i < numParams; ++i) {
- parameters[i].resize(numParams);
for (int j=0; j < paramsSize; ++j)
- parameters[i][j] = i+j;
+ parameters[i*paramsSize+j] = i+j;
} // for
const int tensorSize = 9;
@@ -113,12 +112,13 @@
for (int i=0; i < tensorSize; ++i)
totalStrain[i] = i;
- material._updateState(¶meters, totalStrain);
+ material._updateState(¶meters[0], numParams, &totalStrain[0], tensorSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParams; ++i)
for (int j=0; j < paramsSize; ++j)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i+j), parameters[i][j], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i+j), parameters[i*paramsSize+j],
+ tolerance);
for (int i=0; i < tensorSize; ++i)
CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i), totalStrain[i], tolerance);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStress.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStress.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticPlaneStress.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -100,12 +100,11 @@
const int numParams = data.numParameters;
- std::vector<double_array> parameters(numParams);
+ double_array parameters(numParams);
const int paramsSize = 1;
for (int i=0; i < numParams; ++i) {
- parameters[i].resize(numParams);
for (int j=0; j < paramsSize; ++j)
- parameters[i][j] = i+j;
+ parameters[i*paramsSize+j] = i+j;
} // for
const int tensorSize = 9;
@@ -113,12 +112,13 @@
for (int i=0; i < tensorSize; ++i)
totalStrain[i] = i;
- material._updateState(¶meters, totalStrain);
+ material._updateState(¶meters[0], numParams, &totalStrain[0], tensorSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParams; ++i)
for (int j=0; j < paramsSize; ++j)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i+j), parameters[i][j], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i+j), parameters[i*paramsSize+j],
+ tolerance);
for (int i=0; i < tensorSize; ++i)
CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i), totalStrain[i], tolerance);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStrain1D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStrain1D.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStrain1D.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -100,12 +100,11 @@
const int numParams = data.numParameters;
- std::vector<double_array> parameters(numParams);
+ double_array parameters(numParams);
const int paramsSize = 1;
for (int i=0; i < numParams; ++i) {
- parameters[i].resize(numParams);
for (int j=0; j < paramsSize; ++j)
- parameters[i][j] = i+j;
+ parameters[i*paramsSize+j] = i+j;
} // for
const int tensorSize = 9;
@@ -113,12 +112,13 @@
for (int i=0; i < tensorSize; ++i)
totalStrain[i] = i;
- material._updateState(¶meters, totalStrain);
+ material._updateState(¶meters[0], numParams, &totalStrain[0], tensorSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParams; ++i)
for (int j=0; j < paramsSize; ++j)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i+j), parameters[i][j], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i+j), parameters[i*paramsSize+j],
+ tolerance);
for (int i=0; i < tensorSize; ++i)
CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i), totalStrain[i], tolerance);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStress1D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStress1D.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticStress1D.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -100,12 +100,11 @@
const int numParams = data.numParameters;
- std::vector<double_array> parameters(numParams);
+ double_array parameters(numParams);
const int paramsSize = 1;
for (int i=0; i < numParams; ++i) {
- parameters[i].resize(numParams);
for (int j=0; j < paramsSize; ++j)
- parameters[i][j] = i+j;
+ parameters[i*paramsSize+j] = i+j;
} // for
const int tensorSize = 9;
@@ -113,12 +112,13 @@
for (int i=0; i < tensorSize; ++i)
totalStrain[i] = i;
- material._updateState(¶meters, totalStrain);
+ material._updateState(¶meters[0], numParams, &totalStrain[0], tensorSize);
const double tolerance = 1.0e-06;
for (int i=0; i < numParams; ++i)
for (int j=0; j < paramsSize; ++j)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i+j), parameters[i][j], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i+j), parameters[i*paramsSize+j],
+ tolerance);
for (int i=0; i < tensorSize; ++i)
CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i), totalStrain[i], tolerance);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -189,32 +189,34 @@
// Get cells associated with material
const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
- Mesh::label_sequence::iterator cellIter=cells->begin();
+ Mesh::label_sequence::iterator c_iter = cells->begin();
const double tolerance = 1.0e-06;
- const ALE::Obj<real_section_type>& parameterDensity =
- material._parameters->getReal("density");
- const real_section_type::value_type* densityCell =
- parameterDensity->restrictPoint(*cellIter);
- CPPUNIT_ASSERT(0 != densityCell);
- for (int i=0; i < numQuadPts; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, densityCell[i]/densityE[i], tolerance);
+ const real_section_type::value_type* paramsCell =
+ material._parameters->restrictPoint(*c_iter);
+ CPPUNIT_ASSERT(0 != paramsCell);
+
+ const int pidDensity = 0;
+ const int pidMu = 1;
+ const int pidLambda = 2;
+
+ // density
+ for (int i=0; i < numQuadPts; ++i) {
+ const int index = i*material._numParamsQuadPt + pidDensity;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, paramsCell[index]/densityE[i], tolerance);
+ } // for
- const ALE::Obj<real_section_type>& parameterMu =
- material._parameters->getReal("mu");
- const real_section_type::value_type* muCell =
- parameterMu->restrictPoint(*cellIter);
- CPPUNIT_ASSERT(0 != muCell);
- for (int i=0; i < numQuadPts; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, muCell[i]/muE[i], tolerance);
+ // mu
+ for (int i=0; i < numQuadPts; ++i) {
+ const int index = i*material._numParamsQuadPt + pidMu;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, paramsCell[index]/muE[i], tolerance);
+ } // for
- const ALE::Obj<real_section_type>& parameterLambda =
- material._parameters->getReal("lambda");
- const real_section_type::value_type* lambdaCell =
- parameterLambda->restrictPoint(*cellIter);
- CPPUNIT_ASSERT(0 != lambdaCell);
- for (int i=0; i < numQuadPts; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, lambdaCell[i]/lambdaE[i], tolerance);
+ // lambda
+ for (int i=0; i < numQuadPts; ++i) {
+ const int index = i*material._numParamsQuadPt + pidLambda;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, paramsCell[index]/lambdaE[i], tolerance);
+ } // for
} // testInitialize
// ----------------------------------------------------------------------
@@ -235,28 +237,23 @@
const int numParameters = data.numParameters;
int numParamEntries = 0;
-
- std::vector<double_array> parameterData(numParameters);
- for (int iParam=0; iParam < numParameters; ++iParam) {
- parameterData[iParam].resize(data.numParamValues[iParam]);
+ for (int iParam=0; iParam < numParameters; ++iParam)
numParamEntries += data.numParamValues[iParam];
- } // for
+ double_array parameterData(numParamEntries);
+
double* const parameterDataE = &data.parameterData[iLoc*numParamEntries];
- material->_dbToParameters(¶meterData, dbData);
+ material->_dbToParameters(¶meterData[0], numParamEntries, dbData);
const double tolerance = 1.0e-06;
- for (int iParam=0, i=0; iParam < numParameters; ++iParam) {
- const int numParamValues = data.numParamValues[iParam];
- CPPUNIT_ASSERT_EQUAL(numParamValues, int(parameterData[iParam].size()));
- for (int iValue=0; iValue < numParamValues; ++iValue)
- if (fabs(parameterDataE[i]) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,
- parameterData[iParam][iValue]/parameterDataE[i++],
- tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(parameterDataE[i++], parameterData[iParam][iValue],
- tolerance);
+ for (int i=0; i < numParamEntries; ++i) {
+ if (fabs(parameterDataE[i]) > tolerance)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,
+ parameterData[i]/parameterDataE[i],
+ tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(parameterDataE[i], parameterData[i],
+ tolerance);
} // for
} // for
} // _testDBToParameters
@@ -279,7 +276,7 @@
} // _testDBValues
// ----------------------------------------------------------------------
-// Test _numParameters() and _parameterNames()
+// Test _numParameters().
void
pylith::materials::TestMaterial::_testParameters(Material* material,
const MaterialData& data) const
@@ -291,12 +288,6 @@
const int_array& numParamValues = material->_getNumParamValues();
CPPUNIT_ASSERT_EQUAL(numParameters, int(numParamValues.size()));
- char** const namesE = data.parameterNames;
- const char** names = material->_parameterNames();
- for (int i=0; i < numParameters; ++i) {
- CPPUNIT_ASSERT_EQUAL(data.numParamValues[i], numParamValues[i]);
- CPPUNIT_ASSERT(0 == strcmp(namesE[i], names[i]));
- } // for
} // _testParameters
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -19,7 +19,6 @@
#include "pylith/materials/MaxwellIsotropic3D.hh" // USES MaxwellIsotropic3D
-#include <stdexcept> // TEMPORARY
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestMaxwellIsotropic3D );
@@ -150,46 +149,45 @@
MaxwellIsotropic3D material;
MaxwellIsotropic3DElasticData data;
- const int numParams = data.numParameters;
-
+ const int_array& numParamValues = material._getNumParamValues();
+ const int numParams = numParamValues.size();
+ const int numParamsQuadPt = material._numParamsQuadPt;;
+
const int tensorSize = 6;
double_array totalStrain(tensorSize);
for (int i=0; i < tensorSize; ++i)
totalStrain[i] = i;
- const double meanStrain = (totalStrain[0] + totalStrain[1] + totalStrain[2])/3.0;
+ const double meanStrain =
+ (totalStrain[0] + totalStrain[1] + totalStrain[2]) / 3.0;
- std::vector<double_array> parameters(numParams);
- std::vector<double_array> paramdata(numParams);
- const int paramsSize [] = { 1, 1, 1, 1, 6, 6};
- for (int i=0; i < numParams; ++i) {
- parameters[i].resize(paramsSize[i]);
- paramdata[i].resize(paramsSize[i]);
- for (int j=0; j < paramsSize[i]; ++j)
- parameters[i][j] = i+j;
- } // for
+ double_array parameters(numParamsQuadPt);
+ double_array parametersE(numParamsQuadPt);
+ for (int i=0, index=0; i < numParams; ++i)
+ for (int j=0; j < numParamValues[i]; ++j, ++index) {
+ parametersE[index] = i+j;
+ parameters[index] = i+j;
+ } // for
// Set up vector parameters, which are the only ones updated.
- for (int i=0; i< 3; ++i) {
- paramdata[4][i] = totalStrain[i];
- paramdata[5][i] = totalStrain[i] - meanStrain;
+ const int pidStrainT = 4;
+ const int pidVisStrain = pidStrainT + 6;
+ for (int i=0; i < 3; ++i) {
+ parametersE[pidStrainT+i] = totalStrain[i];
+ parametersE[pidVisStrain+i] = totalStrain[i] - meanStrain;
} // for
- for (int i=3; i< 6; ++i) {
- paramdata[4][i] = totalStrain[i];
- paramdata[5][i] = totalStrain[i];
+ for (int i=3; i < 6; ++i) {
+ parametersE[pidStrainT+i] = totalStrain[i];
+ parametersE[pidVisStrain+i] = totalStrain[i];
} // for
- material._updateState(¶meters, totalStrain);
+ material._updateState(¶meters[0], numParamsQuadPt,
+ &totalStrain[0], tensorSize);
const double tolerance = 1.0e-06;
- // Only test vector parameters
- for (int i=4; i < numParams; ++i)
- for (int j=0; j < paramsSize[i]; ++j)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(paramdata[i][j], parameters[i][j], tolerance);
-
- for (int i=0; i < tensorSize; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i), totalStrain[i], tolerance);
+ for (int i=0; i < numParamsQuadPt; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(parametersE[i], parameters[i], tolerance);
} // testUpdateStateElastic
// ----------------------------------------------------------------------
@@ -200,7 +198,7 @@
MaxwellIsotropic3D material;
MaxwellIsotropic3DTimeDepData data;
material.useElasticBehavior(false);
- double dt = 2.0e5;
+ double dt = 2.0e+5;
material.timeStep(dt);
_testCalcStress(&material, data);
} // testCalcStressTimeDep
@@ -213,7 +211,7 @@
MaxwellIsotropic3D material;
MaxwellIsotropic3DTimeDepData data;
material.useElasticBehavior(false);
- double dt = 2.0e5;
+ double dt = 2.0e+5;
material.timeStep(dt);
_testCalcElasticConsts(&material, data);
} // testElasticConstsTimeDep
@@ -226,14 +224,16 @@
MaxwellIsotropic3D material;
MaxwellIsotropic3DTimeDepData data;
- const int numParams = data.numParameters;
+ const int_array& numParamValues = material._getNumParamValues();
+ const int numParams = numParamValues.size();
+ const int numParamsQuadPt = material._numParamsQuadPt;;
material.useElasticBehavior(false);
- const double dt = 2.0e5;
+ const double dt = 2.0e+5;
material.timeStep(dt);
- const double viscosity = 1.0e18;
- const double mu = 3.0e10;
- const double maxwelltime = viscosity/mu;
+ const double viscosity = 1.0e+18;
+ const double mu = 3.0e+10;
+ const double maxwelltime = viscosity / mu;
const int tensorSize = 6;
double_array totalStrainTpdt(tensorSize);
@@ -241,28 +241,32 @@
double_array visStrainT(tensorSize);
for (int i=0; i < tensorSize; ++i) {
totalStrainTpdt[i] = i;
- totalStrainT[i] = totalStrainTpdt[i]/2.0;
- visStrainT[i] = totalStrainTpdt[i]/4.0;
+ totalStrainT[i] = totalStrainTpdt[i] / 2.0;
+ visStrainT[i] = totalStrainTpdt[i] / 4.0;
} // for
- const double meanStrainTpdt = (totalStrainTpdt[0] + totalStrainTpdt[1] + totalStrainTpdt[2])/3.0;
- const double meanStrainT = (totalStrainT[0] + totalStrainT[1] + totalStrainT[2])/3.0;
+ const double meanStrainTpdt =
+ (totalStrainTpdt[0] + totalStrainTpdt[1] + totalStrainTpdt[2]) / 3.0;
+ const double meanStrainT =
+ (totalStrainT[0] + totalStrainT[1] + totalStrainT[2]) / 3.0;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
- std::vector<double_array> parameters(numParams);
- std::vector<double_array> paramdata(numParams);
- const int paramsSize []= { 1, 1, 1, 1, 6, 6};
- for (int i=0; i < numParams; ++i) {
- parameters[i].resize(paramsSize[i]);
- paramdata[i].resize(paramsSize[i]);
- for (int j=0; j < paramsSize[i]; ++j)
- parameters[i][j] = i+j;
- } // for
+ double_array parameters(numParamsQuadPt);
+ double_array parametersE(numParamsQuadPt);
+ for (int i=0, index=0; i < numParams; ++i)
+ for (int j=0; j < numParamValues[i]; ++j, ++index) {
+ parametersE[index] = i+j;
+ parameters[index] = i+j;
+ } // for
- parameters[3][0] = maxwelltime;
- paramdata[3][0] = maxwelltime;
+ const int pidMaxwellTime = 3;
+ const int pidStrainT = pidMaxwellTime + 1;
+ const int pidVisStrain = pidStrainT + 6;
+ parameters[pidMaxwellTime] = maxwelltime;
+ parametersE[pidMaxwellTime] = maxwelltime;
+
const double dq = maxwelltime*(1.0-exp(-dt/maxwelltime))/dt;
const double expFac = exp(-dt/maxwelltime);
double devStrainTpdt = 0.0;
@@ -271,22 +275,19 @@
for (int i=0; i < tensorSize; ++i) {
devStrainTpdt = totalStrainTpdt[i] - diag[i]*meanStrainTpdt;
devStrainT = totalStrainT[i] - diag[i]*meanStrainT;
- parameters[4][i] = totalStrainT[i];
- parameters[5][i] = visStrainT[i];
- paramdata[4][i] = totalStrainTpdt[i];
- paramdata[5][i] = expFac * visStrainT[i] + dq * (devStrainTpdt - devStrainT);
+ parameters[pidStrainT+i] = totalStrainT[i];
+ parameters[pidVisStrain+i] = visStrainT[i];
+ parametersE[pidStrainT+i] = totalStrainTpdt[i];
+ parametersE[pidVisStrain+i] =
+ expFac * visStrainT[i] + dq * (devStrainTpdt - devStrainT);
} //for
- material._updateState(¶meters, totalStrainTpdt);
+ material._updateState(¶meters[0], numParamsQuadPt,
+ &totalStrainTpdt[0], tensorSize);
const double tolerance = 1.0e-06;
- // Test vector parameters and Maxwell time.
- for (int i=3; i < numParams; ++i)
- for (int j=0; j < paramsSize[i]; ++j)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(paramdata[i][j], parameters[i][j], tolerance);
-
- for (int i=0; i < tensorSize; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i), totalStrainTpdt[i], tolerance);
+ for (int i=0; i < numParamsQuadPt; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(parametersE[i], parameters[i], tolerance);
} // testUpdateStateTimeDep
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -21,6 +21,8 @@
const int pylith::materials::ElasticIsotropic3DData::_numParameters = 3;
+const int pylith::materials::ElasticIsotropic3DData::_numParamsQuadPt = 3;
+
const int pylith::materials::ElasticIsotropic3DData::_numLocs = 2;
const int pylith::materials::ElasticIsotropic3DData::_numParamValues[] = {
@@ -35,12 +37,6 @@
"vp",
};
-const char* pylith::materials::ElasticIsotropic3DData::_parameterNames[] = {
-"density",
-"mu",
-"lambda",
-};
-
const double pylith::materials::ElasticIsotropic3DData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -144,10 +140,10 @@
dimension = _dimension;
numDBValues = _numDBValues;
numParameters = _numParameters;
+ numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
- parameterNames = const_cast<char**>(_parameterNames);
dbData = const_cast<double*>(_dbData);
parameterData = const_cast<double*>(_parameterData);
density = const_cast<double*>(_density);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticIsotropic3DData.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -43,14 +43,14 @@
static const int _numParameters;
+ static const int _numParamsQuadPt;
+
static const int _numLocs;
static const int _numParamValues[];
static const char* _dbValues[];
- static const char* _parameterNames[];
-
static const double _dbData[];
static const double _parameterData[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticMaterialApp.py 2008-01-03 23:00:07 UTC (rev 8984)
@@ -61,9 +61,9 @@
self.dimension = None
self.numDBValues = None
self.numParameters = None
+ self.numParamsQuadPt = None
self.numParamValues = None
self.dbValues = None
- self.parameterNames = None
self.dbData = None
self.parameterData = None
@@ -97,6 +97,8 @@
def _initData(self):
+ self.numParamsQuadPt = numpy.sum(self.numParamValues)
+
self.data.addScalar(vtype="int", name="_dimension",
value=self.dimension,
format="%d")
@@ -106,14 +108,14 @@
self.data.addScalar(vtype="int", name="_numParameters",
value=self.numParameters,
format="%d")
+ self.data.addScalar(vtype="int", name="_numParamsQuadPt",
+ value=self.numParamsQuadPt,
+ format="%d")
self.data.addArray(vtype="int", name="_numParamValues",
values=self.numParamValues,
format="%d", ncols=1)
self.data.addArray(vtype="char*", name="_dbValues", values=self.dbValues,
format="\"%s\"", ncols=1)
- self.data.addArray(vtype="char*", name="_parameterNames",
- values=self.parameterNames,
- format="\"%s\"", ncols=1)
self.data.addArray(vtype="double", name="_dbData", values=self.dbData,
format="%16.8e", ncols=1)
self.data.addArray(vtype="double", name="_parameterData",
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -21,6 +21,8 @@
const int pylith::materials::ElasticPlaneStrainData::_numParameters = 3;
+const int pylith::materials::ElasticPlaneStrainData::_numParamsQuadPt = 3;
+
const int pylith::materials::ElasticPlaneStrainData::_numLocs = 2;
const int pylith::materials::ElasticPlaneStrainData::_numParamValues[] = {
@@ -35,12 +37,6 @@
"vp",
};
-const char* pylith::materials::ElasticPlaneStrainData::_parameterNames[] = {
-"density",
-"mu",
-"lambda",
-};
-
const double pylith::materials::ElasticPlaneStrainData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -102,10 +98,10 @@
dimension = _dimension;
numDBValues = _numDBValues;
numParameters = _numParameters;
+ numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
- parameterNames = const_cast<char**>(_parameterNames);
dbData = const_cast<double*>(_dbData);
parameterData = const_cast<double*>(_parameterData);
density = const_cast<double*>(_density);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStrainData.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -43,14 +43,14 @@
static const int _numParameters;
+ static const int _numParamsQuadPt;
+
static const int _numLocs;
static const int _numParamValues[];
static const char* _dbValues[];
- static const char* _parameterNames[];
-
static const double _dbData[];
static const double _parameterData[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -21,6 +21,8 @@
const int pylith::materials::ElasticPlaneStressData::_numParameters = 3;
+const int pylith::materials::ElasticPlaneStressData::_numParamsQuadPt = 3;
+
const int pylith::materials::ElasticPlaneStressData::_numLocs = 2;
const int pylith::materials::ElasticPlaneStressData::_numParamValues[] = {
@@ -35,12 +37,6 @@
"vp",
};
-const char* pylith::materials::ElasticPlaneStressData::_parameterNames[] = {
-"density",
-"mu",
-"lambda",
-};
-
const double pylith::materials::ElasticPlaneStressData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -102,10 +98,10 @@
dimension = _dimension;
numDBValues = _numDBValues;
numParameters = _numParameters;
+ numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
- parameterNames = const_cast<char**>(_parameterNames);
dbData = const_cast<double*>(_dbData);
parameterData = const_cast<double*>(_parameterData);
density = const_cast<double*>(_density);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticPlaneStressData.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -43,14 +43,14 @@
static const int _numParameters;
+ static const int _numParamsQuadPt;
+
static const int _numLocs;
static const int _numParamValues[];
static const char* _dbValues[];
- static const char* _parameterNames[];
-
static const double _dbData[];
static const double _parameterData[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -21,6 +21,8 @@
const int pylith::materials::ElasticStrain1DData::_numParameters = 2;
+const int pylith::materials::ElasticStrain1DData::_numParamsQuadPt = 2;
+
const int pylith::materials::ElasticStrain1DData::_numLocs = 2;
const int pylith::materials::ElasticStrain1DData::_numParamValues[] = {
@@ -33,11 +35,6 @@
"vp",
};
-const char* pylith::materials::ElasticStrain1DData::_parameterNames[] = {
-"density",
-"lambda2mu",
-};
-
const double pylith::materials::ElasticStrain1DData::_dbData[] = {
2.50000000e+03,
5.00000000e+03,
@@ -77,10 +74,10 @@
dimension = _dimension;
numDBValues = _numDBValues;
numParameters = _numParameters;
+ numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
- parameterNames = const_cast<char**>(_parameterNames);
dbData = const_cast<double*>(_dbData);
parameterData = const_cast<double*>(_parameterData);
density = const_cast<double*>(_density);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStrain1DData.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -43,14 +43,14 @@
static const int _numParameters;
+ static const int _numParamsQuadPt;
+
static const int _numLocs;
static const int _numParamValues[];
static const char* _dbValues[];
- static const char* _parameterNames[];
-
static const double _dbData[];
static const double _parameterData[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -21,6 +21,8 @@
const int pylith::materials::ElasticStress1DData::_numParameters = 3;
+const int pylith::materials::ElasticStress1DData::_numParamsQuadPt = 3;
+
const int pylith::materials::ElasticStress1DData::_numLocs = 2;
const int pylith::materials::ElasticStress1DData::_numParamValues[] = {
@@ -35,12 +37,6 @@
"vp",
};
-const char* pylith::materials::ElasticStress1DData::_parameterNames[] = {
-"density",
-"mu",
-"lambda",
-};
-
const double pylith::materials::ElasticStress1DData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -84,10 +80,10 @@
dimension = _dimension;
numDBValues = _numDBValues;
numParameters = _numParameters;
+ numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
- parameterNames = const_cast<char**>(_parameterNames);
dbData = const_cast<double*>(_dbData);
parameterData = const_cast<double*>(_parameterData);
density = const_cast<double*>(_density);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/ElasticStress1DData.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -43,14 +43,14 @@
static const int _numParameters;
+ static const int _numParamsQuadPt;
+
static const int _numLocs;
static const int _numParamValues[];
static const char* _dbValues[];
- static const char* _parameterNames[];
-
static const double _dbData[];
static const double _parameterData[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -19,9 +19,9 @@
numLocs(0),
numDBValues(0),
numParameters(0),
+ numParamsQuadPt(0),
numParamValues(0),
- dbValues(0),
- parameterNames(0)
+ dbValues(0)
{ // constructor
} // constructor
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaterialData.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -39,10 +39,10 @@
int numDBValues; ///< Number of database values
int numParameters; ///< Number of parameters
+ int numParamsQuadPt; ///< Number of parameters per quadrature point.
int* numParamValues; ///< Number of values for each parameter
char** dbValues; ///< Aray of names of database values;
- char** parameterNames; //< Array of names of parameters
double* dbData; ///< Array of database values at locations
double* parameterData; ///< Array of parameter values at locations
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -21,6 +21,8 @@
const int pylith::materials::MaxwellIsotropic3DElasticData::_numParameters = 6;
+const int pylith::materials::MaxwellIsotropic3DElasticData::_numParamsQuadPt = 16;
+
const int pylith::materials::MaxwellIsotropic3DElasticData::_numLocs = 2;
const int pylith::materials::MaxwellIsotropic3DElasticData::_numParamValues[] = {
@@ -39,15 +41,6 @@
"viscosity",
};
-const char* pylith::materials::MaxwellIsotropic3DElasticData::_parameterNames[] = {
-"density",
-"mu",
-"lambda",
-"maxwellTime",
-"strainT",
-"visStrain",
-};
-
const double pylith::materials::MaxwellIsotropic3DElasticData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -179,10 +172,10 @@
dimension = _dimension;
numDBValues = _numDBValues;
numParameters = _numParameters;
+ numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
- parameterNames = const_cast<char**>(_parameterNames);
dbData = const_cast<double*>(_dbData);
parameterData = const_cast<double*>(_parameterData);
density = const_cast<double*>(_density);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -43,14 +43,14 @@
static const int _numParameters;
+ static const int _numParamsQuadPt;
+
static const int _numLocs;
static const int _numParamValues[];
static const char* _dbValues[];
- static const char* _parameterNames[];
-
static const double _dbData[];
static const double _parameterData[];
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc 2008-01-03 23:00:07 UTC (rev 8984)
@@ -21,6 +21,8 @@
const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numParameters = 6;
+const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numParamsQuadPt = 16;
+
const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numLocs = 2;
const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numParamValues[] = {
@@ -39,15 +41,6 @@
"viscosity",
};
-const char* pylith::materials::MaxwellIsotropic3DTimeDepData::_parameterNames[] = {
-"density",
-"mu",
-"lambda",
-"maxwellTime",
-"strainT",
-"visStrain",
-};
-
const double pylith::materials::MaxwellIsotropic3DTimeDepData::_dbData[] = {
2.50000000e+03,
3.00000000e+03,
@@ -179,10 +172,10 @@
dimension = _dimension;
numDBValues = _numDBValues;
numParameters = _numParameters;
+ numParamsQuadPt = _numParamsQuadPt;
numLocs = _numLocs;
numParamValues = const_cast<int*>(_numParamValues);
dbValues = const_cast<char**>(_dbValues);
- parameterNames = const_cast<char**>(_parameterNames);
dbData = const_cast<double*>(_dbData);
parameterData = const_cast<double*>(_parameterData);
density = const_cast<double*>(_density);
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh 2008-01-03 22:44:59 UTC (rev 8983)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh 2008-01-03 23:00:07 UTC (rev 8984)
@@ -43,14 +43,14 @@
static const int _numParameters;
+ static const int _numParamsQuadPt;
+
static const int _numLocs;
static const int _numParamValues[];
static const char* _dbValues[];
- static const char* _parameterNames[];
-
static const double _dbData[];
static const double _parameterData[];
More information about the cig-commits
mailing list