[cig-commits] r19509 - in short/3D/PyLith/trunk: libsrc/pylith/materials tests/3d/plasticity tests/3d/refine

brad at geodynamics.org brad at geodynamics.org
Sat Jan 28 13:33:24 PST 2012


Author: brad
Date: 2012-01-28 13:33:23 -0800 (Sat, 28 Jan 2012)
New Revision: 19509

Added:
   short/3D/PyLith/trunk/tests/3d/plasticity/dynamic/
Modified:
   short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.hh
   short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.icc
   short/3D/PyLith/trunk/tests/3d/refine/tet4.exo.gz
   short/3D/PyLith/trunk/tests/3d/refine/tet4.jou
Log:
Merge from stable.

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc	2012-01-28 21:32:37 UTC (rev 19508)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/DruckerPrager3D.cc	2012-01-28 21:33:23 UTC (rev 19509)
@@ -65,10 +65,14 @@
 
       // Values expected in properties spatial database
       const int numDBProperties = 6;
-      const char* dbProperties[] = {"density", "vs", "vp" ,
-				    "friction-angle",
-				    "cohesion",
-				    "dilatation-angle"};
+      const char* dbProperties[6] = {
+	"density", 
+	"vs", 
+	"vp" ,
+	"friction-angle",
+	"cohesion",
+	"dilatation-angle",
+      };
 
       /// Number of state variables.
       const int numStateVars = 1;
@@ -80,12 +84,13 @@
 
       // Values expected in state variables spatial database.
       const int numDBStateVars = 6;
-      const char* dbStateVars[] = { "plastic-strain-xx",
-				    "plastic-strain-yy",
-				    "plastic-strain-zz",
-				    "plastic-strain-xy",
-				    "plastic-strain-yz",
-				    "plastic-strain-xz"
+      const char* dbStateVars[6] = {
+	"plastic-strain-xx",
+	"plastic-strain-yy",
+	"plastic-strain-zz",
+	"plastic-strain-xy",
+	"plastic-strain-yz",
+	"plastic-strain-xz",
       };
 
     } // _DruckerPrager3D
@@ -200,7 +205,7 @@
 				PylithScalar* const propValues,
 				const scalar_array& dbValues)
 { // _dbToProperties
-  assert(0 != propValues);
+  assert(propValues);
   const int numDBValues = dbValues.size();
   assert(_DruckerPrager3D::numDBProperties == numDBValues);
 
@@ -211,9 +216,12 @@
   const PylithScalar cohesion = dbValues[db_cohesion];
   const PylithScalar dilatationAngle = dbValues[db_dilatationAngle];
  
-  if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || frictionAngle < 0.0
-      || cohesion < 0.0 || dilatationAngle < 0.0
-      || frictionAngle < dilatationAngle) {
+  const double pi = M_PI;
+
+  if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || cohesion < 0.0 ||
+      frictionAngle < 0.0 || frictionAngle > pi/2 ||
+      dilatationAngle < 0.0 || dilatationAngle > pi/2 ||
+      frictionAngle < dilatationAngle) {
     std::ostringstream msg;
     msg << "Spatial database returned illegal value for physical "
 	<< "properties.\n"
@@ -290,8 +298,8 @@
 pylith::materials::DruckerPrager3D::_nondimProperties(PylithScalar* const values,
 					         const int nvalues) const
 { // _nondimProperties
-  assert(0 != _normalizer);
-  assert(0 != values);
+  assert(_normalizer);
+  assert(values);
   assert(nvalues == _numPropsQuadPt);
 
   const PylithScalar densityScale = _normalizer->densityScale();
@@ -316,8 +324,8 @@
 pylith::materials::DruckerPrager3D::_dimProperties(PylithScalar* const values,
 						      const int nvalues) const
 { // _dimProperties
-  assert(0 != _normalizer);
-  assert(0 != values);
+  assert(_normalizer);
+  assert(values);
   assert(nvalues == _numPropsQuadPt);
 
   const PylithScalar densityScale = _normalizer->densityScale();
@@ -342,7 +350,7 @@
 				PylithScalar* const stateValues,
 				const scalar_array& dbValues)
 { // _dbToStateVars
-  assert(0 != stateValues);
+  assert(stateValues);
   const int numDBValues = dbValues.size();
   assert(_DruckerPrager3D::numDBStateVars == numDBValues);
 
@@ -351,8 +359,6 @@
   assert(totalSize == numDBValues);
   memcpy(&stateValues[s_plasticStrain], &dbValues[db_plasticStrain],
 	 _tensorSize*sizeof(PylithScalar));
-
-  PetscLogFlops(0);
 } // _dbToStateVars
 
 // ----------------------------------------------------------------------
@@ -361,11 +367,9 @@
 pylith::materials::DruckerPrager3D::_nondimStateVars(PylithScalar* const values,
 						const int nvalues) const
 { // _nondimStateVars
-  assert(0 != _normalizer);
-  assert(0 != values);
+  assert(_normalizer);
+  assert(values);
   assert(nvalues == _numVarsQuadPt);
-
-  PetscLogFlops(0);
 } // _nondimStateVars
 
 // ----------------------------------------------------------------------
@@ -374,11 +378,9 @@
 pylith::materials::DruckerPrager3D::_dimStateVars(PylithScalar* const values,
 					     const int nvalues) const
 { // _dimStateVars
-  assert(0 != _normalizer);
-  assert(0 != values);
+  assert(_normalizer);
+  assert(values);
   assert(nvalues == _numVarsQuadPt);
-
-  PetscLogFlops(0);
 } // _dimStateVars
 
 // ----------------------------------------------------------------------
@@ -390,8 +392,8 @@
 					    const PylithScalar* stateVars,
 					    const int numStateVars)
 { // _calcDensity
-  assert(0 != density);
-  assert(0 != properties);
+  assert(density);
+  assert(properties);
   assert(_numPropsQuadPt == numProperties);
 
   density[0] = properties[p_density];
@@ -406,14 +408,15 @@
 				  const PylithScalar* stateVars,
 				  const int numStateVars) const
 { // _stableTimeStepImplicit
-  assert(0 != properties);
+  assert(properties);
   assert(_numPropsQuadPt == numProperties);
-  assert(0 != stateVars);
+  assert(stateVars);
   assert(_numVarsQuadPt == numStateVars);
+
   // It's unclear what to do for an elasto-plastic material, which has no
   // inherent time scale. For now, just set dtStable to a large value.
   const PylithScalar dtStable = pylith::PYLITH_MAXDOUBLE;
-  PetscLogFlops(0);
+
   return dtStable;
 } // _stableTimeStepImplicit
 
@@ -436,17 +439,17 @@
 					 const int initialStrainSize,
 					 const bool computeStateVars)
 { // _calcStressElastic
-  assert(0 != stress);
+  assert(stress);
   assert(_DruckerPrager3D::tensorSize == stressSize);
-  assert(0 != properties);
+  assert(properties);
   assert(_numPropsQuadPt == numProperties);
-  assert(0 != stateVars);
+  assert(stateVars);
   assert(_numVarsQuadPt == numStateVars);
-  assert(0 != totalStrain);
+  assert(totalStrain);
   assert(_DruckerPrager3D::tensorSize == strainSize);
-  assert(0 != initialStress);
+  assert(initialStress);
   assert(_DruckerPrager3D::tensorSize == initialStressSize);
-  assert(0 != initialStrain);
+  assert(initialStrain);
   assert(_DruckerPrager3D::tensorSize == initialStrainSize);
 
   const PylithScalar mu = properties[p_mu];
@@ -492,20 +495,21 @@
 					const int initialStrainSize,
 					const bool computeStateVars)
 { // _calcStressElastoplastic
-  assert(0 != stress);
+  assert(stress);
   assert(_DruckerPrager3D::tensorSize == stressSize);
-  assert(0 != properties);
+  assert(properties);
   assert(_numPropsQuadPt == numProperties);
-  assert(0 != stateVars);
+  assert(stateVars);
   assert(_numVarsQuadPt == numStateVars);
-  assert(0 != totalStrain);
+  assert(totalStrain);
   assert(_DruckerPrager3D::tensorSize == strainSize);
-  assert(0 != initialStress);
+  assert(initialStress);
   assert(_DruckerPrager3D::tensorSize == initialStressSize);
-  assert(0 != initialStrain);
+  assert(initialStrain);
   assert(_DruckerPrager3D::tensorSize == initialStrainSize);
 
-  const int tensorSize = _tensorSize;
+  const int tensorSize = 6;
+  assert(_tensorSize == tensorSize);
   const PylithScalar mu = properties[p_mu];
   const PylithScalar lambda = properties[p_lambda];
     
@@ -521,39 +525,35 @@
     const PylithScalar ae = 1.0/mu2;
     const PylithScalar am = 1.0/(3.0 * bulkModulus);
 
-    const PylithScalar plasticStrainT[] = {stateVars[s_plasticStrain],
-				     stateVars[s_plasticStrain + 1],
-				     stateVars[s_plasticStrain + 2],
-				     stateVars[s_plasticStrain + 3],
-				     stateVars[s_plasticStrain + 4],
-				     stateVars[s_plasticStrain + 5]};
+    const PylithScalar plasticStrainT[tensorSize] = {
+      stateVars[s_plasticStrain  ],
+      stateVars[s_plasticStrain+1],
+      stateVars[s_plasticStrain+2],
+      stateVars[s_plasticStrain+3],
+      stateVars[s_plasticStrain+4],
+      stateVars[s_plasticStrain+5],
+    };
     const PylithScalar meanPlasticStrainT = (plasticStrainT[0] +
 				       plasticStrainT[1] +
 				       plasticStrainT[2])/3.0;
     PylithScalar devPlasticStrainT[tensorSize];
-    pylith::materials::ElasticMaterial::calcDeviatoric3D(devPlasticStrainT,
-							 plasticStrainT,
-							 meanPlasticStrainT);
+    calcDeviatoric3D(devPlasticStrainT, plasticStrainT, meanPlasticStrainT);
 
-    const PylithScalar diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+    const PylithScalar diag[tensorSize] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
 
     // Initial stress values
     const PylithScalar meanStressInitial = (initialStress[0] +
 				      initialStress[1] +
 				      initialStress[2])/3.0;
     PylithScalar devStressInitial[tensorSize];
-    pylith::materials::ElasticMaterial::calcDeviatoric3D(devStressInitial,
-							 initialStress,
-							 meanStressInitial);
+    calcDeviatoric3D(devStressInitial, initialStress, meanStressInitial);
 
     // Initial strain values
     const PylithScalar meanStrainInitial = (initialStrain[0] +
 				      initialStrain[1] +
 				      initialStrain[2])/3.0;
     PylithScalar devStrainInitial[tensorSize];
-    pylith::materials::ElasticMaterial::calcDeviatoric3D(devStrainInitial,
-							 initialStrain,
-							 meanStrainInitial);
+    calcDeviatoric3D(devStrainInitial, initialStrain, meanStrainInitial);
 
     // Values for current time step
     const PylithScalar e11 = totalStrain[0];
@@ -563,30 +563,28 @@
     const PylithScalar meanStrainPPTpdt = meanStrainTpdt - meanPlasticStrainT -
       meanStrainInitial;
 
-    const PylithScalar strainPPTpdt[] =
-      { totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
-	devStrainInitial[0],
-	totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
-	devStrainInitial[1],
-	totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] -
-	devStrainInitial[2],
-	totalStrain[3] - devPlasticStrainT[3] - devStrainInitial[3],
-	totalStrain[4] - devPlasticStrainT[4] - devStrainInitial[4],
-	totalStrain[5] - devPlasticStrainT[5] - devStrainInitial[5] };
+    const PylithScalar strainPPTpdt[tensorSize] = {
+      totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
+      totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
+      totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+      totalStrain[3] - devPlasticStrainT[3] - devStrainInitial[3],
+      totalStrain[4] - devPlasticStrainT[4] - devStrainInitial[4],
+      totalStrain[5] - devPlasticStrainT[5] - devStrainInitial[5],
+    };
 
     // Compute trial elastic stresses and yield function to see if yield should
     // occur.
-    const PylithScalar trialDevStress[] = { strainPPTpdt[0]/ae + devStressInitial[0],
-				      strainPPTpdt[1]/ae + devStressInitial[1],
-				      strainPPTpdt[2]/ae + devStressInitial[2],
-				      strainPPTpdt[3]/ae + devStressInitial[3],
-				      strainPPTpdt[4]/ae + devStressInitial[4],
-				      strainPPTpdt[5]/ae + devStressInitial[5]};
+    const PylithScalar trialDevStress[tensorSize] = { 
+      strainPPTpdt[0]/ae + devStressInitial[0],
+      strainPPTpdt[1]/ae + devStressInitial[1],
+      strainPPTpdt[2]/ae + devStressInitial[2],
+      strainPPTpdt[3]/ae + devStressInitial[3],
+      strainPPTpdt[4]/ae + devStressInitial[4],
+      strainPPTpdt[5]/ae + devStressInitial[5],
+    };
     const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
     const PylithScalar stressInvar2 =
-      sqrt(0.5 *
-	   pylith::materials::ElasticMaterial::scalarProduct3D(trialDevStress,
-							       trialDevStress));
+      sqrt(0.5 * scalarProduct3D(trialDevStress, trialDevStress));
     const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress +
       stressInvar2 - beta;
 #if 0 // DEBUGGING
@@ -602,16 +600,12 @@
     // If yield function is greater than zero, compute elastoplastic stress.
     if (yieldFunction >= 0.0) {
       const PylithScalar devStressInitialProd = 
-	pylith::materials::ElasticMaterial::scalarProduct3D(devStressInitial,
-							    devStressInitial);
-      const PylithScalar strainPPTpdtProd =
-	pylith::materials::ElasticMaterial::scalarProduct3D(strainPPTpdt,
-							    strainPPTpdt);
+	scalarProduct3D(devStressInitial, devStressInitial);
+      const PylithScalar strainPPTpdtProd = 
+	scalarProduct3D(strainPPTpdt, strainPPTpdt);
       const PylithScalar d =
 	sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
-	   pylith::materials::ElasticMaterial::scalarProduct3D(devStressInitial,
-							       strainPPTpdt) +
-	     strainPPTpdtProd);
+	   scalarProduct3D(devStressInitial, strainPPTpdt) + strainPPTpdtProd);
       const PylithScalar plasticMult = 2.0 * ae * am *
 	(3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
 	(6.0 * alphaYield * alphaFlow * ae + am);
@@ -641,16 +635,18 @@
       stress[5] = strainPPTpdt[5]/ae + devStressInitial[5]; 
     } // if
 
+  } else {
     // If state variables have already been updated, the plastic strain for the
     // time step has already been computed.
-  } else {
     const PylithScalar mu2 = 2.0 * mu;
-    const PylithScalar plasticStrainTpdt[] = {stateVars[s_plasticStrain],
-					stateVars[s_plasticStrain + 1],
-					stateVars[s_plasticStrain + 2],
-					stateVars[s_plasticStrain + 3],
-					stateVars[s_plasticStrain + 4],
-					stateVars[s_plasticStrain + 5]};
+    const PylithScalar plasticStrainTpdt[tensorSize] = {
+      stateVars[s_plasticStrain  ],
+      stateVars[s_plasticStrain+1],
+      stateVars[s_plasticStrain+2],
+      stateVars[s_plasticStrain+3],
+      stateVars[s_plasticStrain+4],
+      stateVars[s_plasticStrain+5],
+    };
 
     const PylithScalar e11 = totalStrain[0] - plasticStrainTpdt[0] - initialStrain[0];
     const PylithScalar e22 = totalStrain[1] - plasticStrainTpdt[1] - initialStrain[1];
@@ -671,7 +667,7 @@
 
     PetscLogFlops(31);
 
-  } // else
+  } // if/else
 #if 0 // DEBUGGING
   const PylithScalar alphaYield = properties[p_alphaYield];
   const PylithScalar beta = properties[p_beta];
@@ -684,9 +680,7 @@
 				   stress[4],
 				   stress[5]};
   const PylithScalar stressInvar2Test =
-    sqrt(0.5 *
-	 pylith::materials::ElasticMaterial::scalarProduct3D(devStressTest,
-							     devStressTest));
+    sqrt(0.5 * scalarProduct3D(devStressTest, devStressTest));
   
   const PylithScalar yieldFunctionTest = 3.0 * alphaYield * meanStressTest +
       stressInvar2Test - beta;
@@ -809,7 +803,7 @@
 
   // Duplicate functionality of _calcStressElastoplastic
   // Get properties
-  const int tensorSize = _tensorSize;
+  const int tensorSize = 6;
   const PylithScalar mu = properties[p_mu];
   const PylithScalar lambda = properties[p_lambda];
   const PylithScalar alphaYield = properties[p_alphaYield];
@@ -821,19 +815,19 @@
   const PylithScalar am = 1.0/(3.0 * bulkModulus);
   
   // Get state variables from previous time step
-  const PylithScalar plasticStrainT[] = {stateVars[s_plasticStrain],
-				   stateVars[s_plasticStrain + 1],
-				   stateVars[s_plasticStrain + 2],
-				   stateVars[s_plasticStrain + 3],
-				   stateVars[s_plasticStrain + 4],
-				   stateVars[s_plasticStrain + 5]};
+  const PylithScalar plasticStrainT[tensorSize] = {
+    stateVars[s_plasticStrain  ],
+    stateVars[s_plasticStrain+1],
+    stateVars[s_plasticStrain+2],
+    stateVars[s_plasticStrain+3],
+    stateVars[s_plasticStrain+4],
+    stateVars[s_plasticStrain+5],
+  };
   const PylithScalar meanPlasticStrainT = (plasticStrainT[0] +
 				     plasticStrainT[1] +
 				     plasticStrainT[2])/3.0;
   PylithScalar devPlasticStrainT[tensorSize];
-  pylith::materials::ElasticMaterial::calcDeviatoric3D(devPlasticStrainT,
-						       plasticStrainT,
-						       meanPlasticStrainT);
+  calcDeviatoric3D(devPlasticStrainT, plasticStrainT, meanPlasticStrainT);
 
   const PylithScalar diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
 
@@ -842,18 +836,14 @@
 				    initialStress[1] +
 				    initialStress[2])/3.0;
   PylithScalar devStressInitial[tensorSize];
-  pylith::materials::ElasticMaterial::calcDeviatoric3D(devStressInitial,
-						       initialStress,
-						       meanStressInitial);
+  calcDeviatoric3D(devStressInitial, initialStress, meanStressInitial);
 
   // Initial strain values
   const PylithScalar meanStrainInitial = (initialStrain[0] +
 				    initialStrain[1] +
 				    initialStrain[2])/3.0;
   PylithScalar devStrainInitial[tensorSize];
-  pylith::materials::ElasticMaterial::calcDeviatoric3D(devStrainInitial,
-						       initialStrain,
-						       meanStrainInitial);
+  calcDeviatoric3D(devStrainInitial, initialStrain, meanStrainInitial);
 
   // Values for current time step
   const PylithScalar meanStrainTpdt = (totalStrain[0] +
@@ -883,9 +873,7 @@
 				    strainPPTpdt[5]/ae + devStressInitial[5]};
   const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
   const PylithScalar stressInvar2 =
-    sqrt(0.5 *
-	 pylith::materials::ElasticMaterial::scalarProduct3D(trialDevStress,
-							     trialDevStress));
+    sqrt(0.5 * scalarProduct3D(trialDevStress, trialDevStress));
   const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress +
     stressInvar2 - beta;
 #if 0 // DEBUGGING
@@ -902,16 +890,12 @@
   // corresponding tangent matrix.
   if (yieldFunction >= 0.0) {
     const PylithScalar devStressInitialProd = 
-      pylith::materials::ElasticMaterial::scalarProduct3D(devStressInitial,
-							  devStressInitial);
+      scalarProduct3D(devStressInitial, devStressInitial);
     const PylithScalar strainPPTpdtProd =
-      pylith::materials::ElasticMaterial::scalarProduct3D(strainPPTpdt,
-							  strainPPTpdt);
+      scalarProduct3D(strainPPTpdt, strainPPTpdt);
     const PylithScalar d = 
       sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
-	   pylith::materials::ElasticMaterial::scalarProduct3D(devStressInitial,
-							       strainPPTpdt) +
-	   strainPPTpdtProd);
+	   scalarProduct3D(devStressInitial, strainPPTpdt) + strainPPTpdtProd);
     const PylithScalar plasticFac = 2.0 * ae * am/
       (6.0 * alphaYield * alphaFlow * ae + am);
     const PylithScalar meanStrainFac = 3.0 * alphaYield/am;
@@ -1098,9 +1082,7 @@
 				     plasticStrainT[1] +
 				     plasticStrainT[2])/3.0;
   PylithScalar devPlasticStrainT[tensorSize];
-  pylith::materials::ElasticMaterial::calcDeviatoric3D(devPlasticStrainT,
-						       plasticStrainT,
-						       meanPlasticStrainT);
+  calcDeviatoric3D(devPlasticStrainT, plasticStrainT, meanPlasticStrainT);
 
   const PylithScalar diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
 
@@ -1109,18 +1091,14 @@
 				    initialStress[1] +
 				    initialStress[2])/3.0;
   PylithScalar devStressInitial[tensorSize];
-  pylith::materials::ElasticMaterial::calcDeviatoric3D(devStressInitial,
-						       initialStress,
-						       meanStressInitial);
+  calcDeviatoric3D(devStressInitial, initialStress, meanStressInitial);
 
   // Initial strain values
   const PylithScalar meanStrainInitial = (initialStrain[0] +
 				    initialStrain[1] +
 				    initialStrain[2])/3.0;
   PylithScalar devStrainInitial[tensorSize];
-  pylith::materials::ElasticMaterial::calcDeviatoric3D(devStrainInitial,
-						       initialStrain,
-						       meanStrainInitial);
+  calcDeviatoric3D(devStrainInitial, initialStrain, meanStrainInitial);
 
   // Values for current time step
   const PylithScalar e11 = totalStrain[0];
@@ -1151,9 +1129,7 @@
 				    strainPPTpdt[5]/ae + devStressInitial[5]};
   const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
   const PylithScalar stressInvar2 =
-    sqrt(0.5 *
-	 pylith::materials::ElasticMaterial::scalarProduct3D(trialDevStress,
-							     trialDevStress));
+    sqrt(0.5 * scalarProduct3D(trialDevStress, trialDevStress));
   const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress +
     stressInvar2 - beta;
 #if 0 // DEBUGGING
@@ -1170,16 +1146,12 @@
   // Otherwise, plastic strains remain the same.
   if (yieldFunction >= 0.0) {
     const PylithScalar devStressInitialProd = 
-      pylith::materials::ElasticMaterial::scalarProduct3D(devStressInitial,
-							  devStressInitial);
+      scalarProduct3D(devStressInitial, devStressInitial);
     const PylithScalar strainPPTpdtProd =
-      pylith::materials::ElasticMaterial::scalarProduct3D(strainPPTpdt,
-							  strainPPTpdt);
+      scalarProduct3D(strainPPTpdt, strainPPTpdt);
     const PylithScalar d =
       sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
-	   pylith::materials::ElasticMaterial::scalarProduct3D(devStressInitial,
-							       strainPPTpdt) +
-	   strainPPTpdtProd);
+	   scalarProduct3D(devStressInitial, strainPPTpdt) + strainPPTpdtProd);
     const PylithScalar plasticMult = 2.0 * ae * am *
       (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
       (6.0 * alphaYield * alphaFlow * ae + am);

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc	2012-01-28 21:32:37 UTC (rev 19508)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticIsotropic3D.cc	2012-01-28 21:33:23 UTC (rev 19509)
@@ -211,16 +211,16 @@
 						   const int initialStrainSize,
 						   const bool computeStateVars)
 { // _calcStress
-  assert(0 != stress);
+  assert(stress);
   assert(_ElasticIsotropic3D::tensorSize == stressSize);
-  assert(0 != properties);
+  assert(properties);
   assert(_numPropsQuadPt == numProperties);
   assert(0 == numStateVars);
-  assert(0 != totalStrain);
+  assert(totalStrain);
   assert(_ElasticIsotropic3D::tensorSize == strainSize);
-  assert(0 != initialStress);
+  assert(initialStress);
   assert(_ElasticIsotropic3D::tensorSize == initialStressSize);
-  assert(0 != initialStrain);
+  assert(initialStrain);
   assert(_ElasticIsotropic3D::tensorSize == initialStrainSize);
 
   const PylithScalar mu = properties[p_mu];

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.hh	2012-01-28 21:32:37 UTC (rev 19508)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.hh	2012-01-28 21:33:23 UTC (rev 19509)
@@ -59,42 +59,6 @@
   virtual
   void deallocate(void);
   
-  /** Compute 2D deviatoric stress/strain from vector and mean value.
-   *
-   * @param deviatoric Array for deviatoric tensor.
-   * @param vec Input tensor (as vector).
-   * @param vecMean Tensor trace divided by spatial_dimension.
-   */
-  void calcDeviatoric2D(PylithScalar* const deviatoric,
-			const PylithScalar* vec,
-			const PylithScalar vecMean);
-  
-  /** Compute 3D deviatoric stress/strain from vector and mean value.
-   *
-   * @param deviatoric Array for deviatoric tensor.
-   * @param vec Input tensor (as vector).
-   * @param vecMean Tensor trace divided by spatial_dimension.
-   */
-  void calcDeviatoric3D(PylithScalar* const deviatoric,
-			const PylithScalar* vec,
-			const PylithScalar vecMean);
-  
-  /** Compute 2D scalar product of two tensors represented as vectors.
-   *
-   * @param tensor1 First tensor.
-   * @param tensor2 Second tensor.
-   */
-  PylithScalar scalarProduct2D(const PylithScalar* tensor1,
-			 const PylithScalar* tensor2) const;
-  
-  /** Compute 3D scalar product of two tensors represented as vectors.
-   *
-   * @param tensor1 First tensor.
-   * @param tensor2 Second tensor.
-   */
-  PylithScalar scalarProduct3D(const PylithScalar* tensor1,
-			 const PylithScalar* tensor2) const;
-  
   /** Set database for initial stress state.
    *
    * @param db Spatial database.
@@ -355,6 +319,46 @@
 				 const PylithScalar* stateVars,
 				 const int numStateVars) const = 0;
 
+  /** Compute 2D deviatoric stress/strain from vector and mean value.
+   *
+   * @param deviatoric Array for deviatoric tensor.
+   * @param vec Input tensor (as vector).
+   * @param vecMean Tensor trace divided by spatial_dimension.
+   */
+  static
+  void calcDeviatoric2D(PylithScalar* const deviatoric,
+			const PylithScalar* vec,
+			const PylithScalar vecMean);
+  
+  /** Compute 3D deviatoric stress/strain from vector and mean value.
+   *
+   * @param deviatoric Array for deviatoric tensor.
+   * @param vec Input tensor (as vector).
+   * @param vecMean Tensor trace divided by spatial_dimension.
+   */
+  static
+  void calcDeviatoric3D(PylithScalar* const deviatoric,
+			const PylithScalar* vec,
+			const PylithScalar vecMean);
+  
+  /** Compute 2D scalar product of two tensors represented as vectors.
+   *
+   * @param tensor1 First tensor.
+   * @param tensor2 Second tensor.
+   */
+  static
+  PylithScalar scalarProduct2D(const PylithScalar* tensor1,
+			       const PylithScalar* tensor2);
+  
+  /** Compute 3D scalar product of two tensors represented as vectors.
+   *
+   * @param tensor1 First tensor.
+   * @param tensor2 Second tensor.
+   */
+  static
+  PylithScalar scalarProduct3D(const PylithScalar* tensor1,
+			       const PylithScalar* tensor2);
+  
   // PRIVATE METHODS ////////////////////////////////////////////////////
 private :
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.icc	2012-01-28 21:32:37 UTC (rev 19508)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.icc	2012-01-28 21:33:23 UTC (rev 19509)
@@ -22,6 +22,41 @@
 
 #include <cassert>
 
+// Set database for initial stress state.
+inline
+void
+pylith::materials::ElasticMaterial::dbInitialStress(spatialdata::spatialdb::SpatialDB* db) {
+  _dbInitialStress = db;
+}
+
+// Set database for initial strain state.
+inline
+void
+pylith::materials::ElasticMaterial::dbInitialStrain(spatialdata::spatialdb::SpatialDB* db) {
+  _dbInitialStrain = db;
+}
+
+// Set whether elastic or inelastic constitutive relations are used.
+inline
+void
+pylith::materials::ElasticMaterial::useElasticBehavior(const bool flag) {
+} // useElasticBehavior
+
+// Get flag indicating whether material implements an empty
+// _updateProperties() method.
+inline
+bool
+pylith::materials::ElasticMaterial::hasStateVars(void) const {
+  return _numVarsQuadPt > 0;
+} // usesUpdateProperties
+
+// Get initial stress/strain fields.
+inline
+const pylith::topology::Fields<pylith::topology::Field<pylith::topology::Mesh> >*
+pylith::materials::ElasticMaterial::initialFields(void) const {
+  return _initialFields;
+} // initialFields
+
 // Compute 2D deviatoric stress/strain from vector and mean value.
 // 2 FLOPs per call
 inline
@@ -45,9 +80,9 @@
 					const PylithScalar* vec,
 					const PylithScalar vecMean)
 {
-  deviatoric[0] = vec[0] -vecMean;
-  deviatoric[1] = vec[1] -vecMean;
-  deviatoric[2] = vec[2] -vecMean;
+  deviatoric[0] = vec[0] - vecMean;
+  deviatoric[1] = vec[1] - vecMean;
+  deviatoric[2] = vec[2] - vecMean;
   deviatoric[3] = vec[3];
   deviatoric[4] = vec[4];
   deviatoric[5] = vec[5];
@@ -59,12 +94,14 @@
 inline
 PylithScalar
 pylith::materials::ElasticMaterial::scalarProduct2D(const PylithScalar* tensor1,
-						    const PylithScalar* tensor2) const
+						    const PylithScalar* tensor2)
 { // scalarProduct2D
-  const PylithScalar scalarProduct = tensor1[0] * tensor2[0] +
-    tensor1[1] * tensor2[1] + 2.0 * tensor1[2] * tensor2[2];
+  const PylithScalar scalarProduct = 
+    tensor1[0] * tensor2[0] +
+    tensor1[1] * tensor2[1] + 
+    2.0 * tensor1[2] * tensor2[2];
+
   return scalarProduct;
-
 } // scalarProduct2D
 
 
@@ -73,54 +110,20 @@
 inline
 PylithScalar
 pylith::materials::ElasticMaterial::scalarProduct3D(const PylithScalar* tensor1,
-						    const PylithScalar* tensor2) const
+						    const PylithScalar* tensor2)
 { // scalarProduct3D
-  const PylithScalar scalarProduct = tensor1[0] * tensor2[0] +
+  const PylithScalar scalarProduct = 
+    tensor1[0] * tensor2[0] +
     tensor1[1] * tensor2[1] +
     tensor1[2] * tensor2[2] +
     2.0 * (tensor1[3] * tensor2[3] +
 	   tensor1[4] * tensor2[4] +
 	   tensor1[5] * tensor2[5]);
+
   return scalarProduct;
-
 } // scalarProduct3D
 
 
-// Set database for initial stress state.
-inline
-void
-pylith::materials::ElasticMaterial::dbInitialStress(spatialdata::spatialdb::SpatialDB* db) {
-  _dbInitialStress = db;
-}
-
-// Set database for initial strain state.
-inline
-void
-pylith::materials::ElasticMaterial::dbInitialStrain(spatialdata::spatialdb::SpatialDB* db) {
-  _dbInitialStrain = db;
-}
-
-// Set whether elastic or inelastic constitutive relations are used.
-inline
-void
-pylith::materials::ElasticMaterial::useElasticBehavior(const bool flag) {
-} // useElasticBehavior
-
-// Get flag indicating whether material implements an empty
-// _updateProperties() method.
-inline
-bool
-pylith::materials::ElasticMaterial::hasStateVars(void) const {
-  return _numVarsQuadPt > 0;
-} // usesUpdateProperties
-
-// Get initial stress/strain fields.
-inline
-const pylith::topology::Fields<pylith::topology::Field<pylith::topology::Mesh> >*
-pylith::materials::ElasticMaterial::initialFields(void) const {
-  return _initialFields;
-} // initialFields
-
 // ----------------------------------------------------------------------
 // Compute density for cell at quadrature points.
 inline

Copied: short/3D/PyLith/trunk/tests/3d/plasticity/dynamic (from rev 19507, short/3D/PyLith/branches/v1.6-stable/tests/3d/plasticity/dynamic)

Modified: short/3D/PyLith/trunk/tests/3d/refine/tet4.exo.gz
===================================================================
(Binary files differ)

Modified: short/3D/PyLith/trunk/tests/3d/refine/tet4.jou
===================================================================
--- short/3D/PyLith/trunk/tests/3d/refine/tet4.jou	2012-01-28 21:32:37 UTC (rev 19508)
+++ short/3D/PyLith/trunk/tests/3d/refine/tet4.jou	2012-01-28 21:33:23 UTC (rev 19509)
@@ -7,7 +7,7 @@
 # Set discretization size
 # ----------------------------------------------------------------------
 volume all scheme tetmesh
-volume all size 500
+volume all size 1500
 
 # ----------------------------------------------------------------------
 # Generate the mesh



More information about the CIG-COMMITS mailing list