[cig-commits] r6521 - in short/3D/PyLith/trunk: libsrc libsrc/materials unittests/libtests/materials

[brad at geodynamics.org]

Author: brad
Date: 2007-04-07 19:00:47 -0700 (Sat, 07 Apr 2007)
New Revision: 6521

Added:
   short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc
   short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh
   short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.icc
   short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc
   short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh
   short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.icc
Removed:
   short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.cc
   short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.hh
   short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.icc
Modified:
   short/3D/PyLith/trunk/libsrc/Makefile.am
   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/Makefile.am
   short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
Log:
Moved quadrature loop from specific material to ElasticMaterial (simplifies specific Material routines). Updated unit tests accordingly. Created meaningful 1-D elastic materials. Still need to setup tests for 1-D and 2-D materials.

Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am	2007-04-08 02:00:47 UTC (rev 6521)
@@ -32,7 +32,8 @@
 	feassemble/Quadrature3D.cc \
 	feassemble/ParameterManager.cc \
 	materials/Material.cc \
-	materials/ElasticIsotropic1D.cc \
+	materials/ElasticStress1D.cc \
+	materials/ElasticStrain1D.cc \
 	materials/ElasticIsotropic3D.cc \
 	materials/ElasticMaterial.cc \
 	materials/ElasticPlaneStrain.cc \

Deleted: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.cc	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.cc	2007-04-08 02:00:47 UTC (rev 6521)
@@ -1,236 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-//                           Brad T. Aagaard
-//                        U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-#include <portinfo>
-
-#include "ElasticIsotropic1D.hh" // implementation of object methods
-
-#include <assert.h> // USES assert()
-
-// ----------------------------------------------------------------------
-// _ElasticIsotropic1D is a helper class for ElasticIsotropic1D. We define
-// it in this implementation file to insulate other objects from these
-// details.
-namespace pylith {
-  namespace materials {
-    class _ElasticIsotropic1D;
-  } // materials
-} // pylith
-
-class pylith::materials::_ElasticIsotropic1D {
-public:
-  // Number of entries in stress tensor.
-  static const int stressSize;
-
-  // Number of entries in derivative of elasticity matrix.
-  static const int numElasticConsts;
-
-  // Values expected in spatial database
-  static const int numDBValues;
-  static const char* namesDBValues[];
-
-  // Indices (order) of database values
-  static const int didDensity;
-  static const int didVs;
-  static const int didVp;
-
-  // Parameters
-  static const int numParameters;
-  static const char* namesParameters[];
-
-  // Indices (order) of parameters
-  static const int pidDensity;
-  static const int pidMu;
-  static const int pidLambda;
-}; // _ElasticIsotropic1D
-
-const int pylith::materials::_ElasticIsotropic1D::stressSize = 1;
-const int pylith::materials::_ElasticIsotropic1D::numElasticConsts = 1;
-const int pylith::materials::_ElasticIsotropic1D::numDBValues = 3;
-const char* pylith::materials::_ElasticIsotropic1D::namesDBValues[] =
-  {"density", "vs", "vp" };
-const int pylith::materials::_ElasticIsotropic1D::numParameters = 3;
-const char* pylith::materials::_ElasticIsotropic1D::namesParameters[] =
-  {"density", "mu", "lambda" };
-const int pylith::materials::_ElasticIsotropic1D::didDensity = 0;
-const int pylith::materials::_ElasticIsotropic1D::didVs = 1;
-const int pylith::materials::_ElasticIsotropic1D::didVp = 2;
-const int pylith::materials::_ElasticIsotropic1D::pidDensity = 0;
-const int pylith::materials::_ElasticIsotropic1D::pidMu = 1;
-const int pylith::materials::_ElasticIsotropic1D::pidLambda = 2;
-
-// ----------------------------------------------------------------------
-// Default constructor.
-pylith::materials::ElasticIsotropic1D::ElasticIsotropic1D(void)
-{ // constructor
-  _dimension = 1;
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor.
-pylith::materials::ElasticIsotropic1D::~ElasticIsotropic1D(void)
-{ // destructor
-} // destructor
-
-// ----------------------------------------------------------------------
-// Copy constructor.
-pylith::materials::ElasticIsotropic1D::ElasticIsotropic1D(
-						const ElasticIsotropic1D& m) :
-  ElasticMaterial(m)
-{ // copy constructor
-} // copy constructor
-
-// ----------------------------------------------------------------------
-// Get number of entries in stress tensor.
-int
-pylith::materials::ElasticIsotropic1D::stressSize(void) const
-{ // stressSize
-  return _ElasticIsotropic1D::stressSize;
-} // stressSize
-
-// ----------------------------------------------------------------------
-// Get number of elastic constants for material.
-int
-pylith::materials::ElasticIsotropic1D::numElasticConsts(void) const
-{ // numElasticConsts
-  return _ElasticIsotropic1D::numElasticConsts;
-} // numElasticConsts
-
-// ----------------------------------------------------------------------
-// Get names of values expected to be in database of parameters for
-const char**
-pylith::materials::ElasticIsotropic1D::_dbValues(void) const
-{ // _dbValues
-  return _ElasticIsotropic1D::namesDBValues;
-} // _dbValues
-
-// ----------------------------------------------------------------------
-// Get number of values expected to be in database of parameters for
-int
-pylith::materials::ElasticIsotropic1D::_numDBValues(void) const
-{ // _numDBValues
-  return _ElasticIsotropic1D::numDBValues;
-} // _numDBValues
-
-// ----------------------------------------------------------------------
-// Get names of parameters for physical properties.
-const char**
-pylith::materials::ElasticIsotropic1D::_parameterNames(void) const
-{ // _parameterNames
-  return _ElasticIsotropic1D::namesParameters;
-} // _parameterNames
-
-// ----------------------------------------------------------------------
-// Get number of parameters for physical properties.
-int
-pylith::materials::ElasticIsotropic1D::_numParameters(void) const
-{ // _numParameters
-  return _ElasticIsotropic1D::numParameters;
-} // _numParameters
-
-// ----------------------------------------------------------------------
-// Compute parameters from values in spatial database.
-void
-pylith::materials::ElasticIsotropic1D::_dbToParameters(double* paramVals,
-						       const int numParams,
-						       const double* dbValues,
-						       const int numValues) const
-{ // computeParameters
-  assert(0 != paramVals);
-  assert(_ElasticIsotropic1D::numParameters == numParams);
-  assert(0 != dbValues);
-  assert(_ElasticIsotropic1D::numDBValues == numValues);
-
-  const double density = dbValues[_ElasticIsotropic1D::didDensity];
-  const double vs = dbValues[_ElasticIsotropic1D::didVs];
-  const double vp = dbValues[_ElasticIsotropic1D::didVp];
- 
-  const double mu = density * vs*vs;
-  const double lambda = density * vp*vp - 2.0*mu;
-
-  paramVals[_ElasticIsotropic1D::pidDensity] = density;
-  paramVals[_ElasticIsotropic1D::pidMu] = mu;
-  paramVals[_ElasticIsotropic1D::pidLambda] = lambda;
-} // computeParameters
-
-// ----------------------------------------------------------------------
-// Compute density at location from parameters.
-void
-pylith::materials::ElasticIsotropic1D::_calcDensity(const double* parameters,
-						    const int numParameters,
-						    const int numLocs)
-{ // _calcDensity
-  assert(0 != _density);
-  assert(0 != parameters);
-
-  for (int iLoc=0, index=0; iLoc < numLocs; ++iLoc, index+=numParameters)
-    _density[iLoc] = 
-      parameters[index+_ElasticIsotropic1D::pidDensity];
-} // _calcDensity
-
-// ----------------------------------------------------------------------
-// Compute stress tensor at location from parameters.
-void
-pylith::materials::ElasticIsotropic1D::_calcStress(const double* parameters,
-						   const int numParameters,
-						   const double* totalStrain,
-						   const int numLocs,
-						   const int spaceDim)
-{ // _calcStress
-  assert(0 != _stress);
-  assert(0 != parameters);
-  assert(_ElasticIsotropic1D::numParameters == numParameters);
-  assert(spaceDim == _dimension);
-
-  for (int iLoc=0, indexP=0, indexS=0;
-       iLoc < numLocs; 
-       ++iLoc, 
-	 indexP+=_ElasticIsotropic1D::numParameters,
-	 indexS+=_ElasticIsotropic1D::stressSize) {
-    const double density = parameters[indexP+_ElasticIsotropic1D::pidDensity];
-    const double mu = parameters[indexP+_ElasticIsotropic1D::pidMu];
-    const double lambda = parameters[indexP+_ElasticIsotropic1D::pidLambda];
-
-    const double e11 = totalStrain[indexS  ];
-    _stress[indexS  ] = mu*(3.0*lambda+2.0*mu)/(lambda + mu) * e11;
-  } // for
-} // _calcStress
-
-// ----------------------------------------------------------------------
-// Compute derivative of elasticity matrix at location from parameters.
-void
-pylith::materials::ElasticIsotropic1D::_calcElasticConsts(const double* parameters,
-							  const int numParameters,
-							  const double* totalStrain,
-							  const int numLocs,
-							  const int spaceDim)
-{ // _calcElasticConsts
-  assert(0 != _elasticConsts);
-  assert(0 != parameters);
-  assert(_ElasticIsotropic1D::numParameters == numParameters);
-  assert(spaceDim == _dimension);
- 
-  for (int iLoc=0, indexP=0, indexC=0;
-       iLoc < numLocs; 
-       ++iLoc, 
-       indexP+=_ElasticIsotropic1D::numParameters,
-       indexC+=_ElasticIsotropic1D::numElasticConsts) {
-    const double density = parameters[indexP+_ElasticIsotropic1D::pidDensity];
-    const double mu = parameters[indexP+_ElasticIsotropic1D::pidMu];
-    const double lambda = parameters[indexP+_ElasticIsotropic1D::pidLambda];
-
-    _elasticConsts[indexC  ] = mu*(3.0*lambda+2.0*mu)/(lambda + mu);
-  } // for
-} // _calcElasticConsts
-
-
-// End of file 

Deleted: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.hh	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.hh	2007-04-08 02:00:47 UTC (rev 6521)
@@ -1,188 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-//                           Brad T. Aagaard
-//                        U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-/** @file libsrc/materials/ElasticIsotropic1D.h
- *
- * @brief C++ ElasticIsotropic1D object
- *
- * 1-D, isotropic, linear elastic material. The physical properties
- * are specified using density and compressional-wave speed. The
- * physical properties are stored internally using density and lambda
- * + 2 mu, which are directly related to the elasticity constants used
- * in the finite-element integration.
- */
-
-#if !defined(pylith_materials_elasticisotropic1d_hh)
-#define pylith_materials_elasticisotropic1d_hh
-
-#include "ElasticMaterial.hh"
-
-/// Namespace for pylith package
-namespace pylith {
-  namespace materials {
-    class ElasticIsotropic1D;
-    class TestElasticIsotropic1D; // unit testing
-  } // materials
-} // pylith
-
-/// 3-D, isotropic, linear elastic material.
-class pylith::materials::ElasticIsotropic1D : public ElasticMaterial
-{ // class ElasticIsotropic1D
-  friend class TestElasticIsotropic1D; // unit testing
-
-  // PUBLIC METHODS /////////////////////////////////////////////////////
-public :
-
-  /// Default constructor
-  ElasticIsotropic1D(void);
-
-  /// Destructor
-  ~ElasticIsotropic1D(void);
-
-  /** Create a pointer to a copy of this.
-   *
-   * @returns Pointer to copy
-   */
-  ElasticMaterial* clone(void) const;
-
-  /** Get number of entries in stress tensor.
-   *
-   * 1-D = 1
-   * 2-D = 3
-   * 3-D = 6
-   *
-   * @returns Number of entries in stress tensor.
-   */
-  int stressSize(void) const;
-
-  /** Get number of elastic constants for material.
-   *
-   * 1-D = 1
-   * 2-D = 6
-   * 3-D = 21
-   *
-   * @returns Number of elastic constants
-   */
-  int numElasticConsts(void) const;
-
-  // PROTECTED METHODS //////////////////////////////////////////////////
-protected :
-
-  /** Copy constructor
-   *
-   * @param m Material to copy
-   */
-  ElasticIsotropic1D(const ElasticIsotropic1D& m);
-
-  /** Get names of values expected to be in database of parameters for
-   *  physical properties.
-   *
-   * @returns Names of values
-   */
-  const char** _dbValues(void) const;
-
-  /** Get number of values expected to be in database of parameters for
-   *  physical properties.
-   *
-   * @returns Number of values
-   */
-  int _numDBValues(void) const;
-
-  /** Get names of parameters for physical properties.
-   *
-   * @returns Names of parameters
-   */
-  const char** _parameterNames(void) const;
-
-  /** Get number of parameters for physical properties.
-   *
-   * @returns Number of parameters
-   */
-  int _numParameters(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
-   * @param dbValues Array of database values
-   * @param numValues Number of database values
-   */
-  void _dbToParameters(double* paramVals,
-		       const int numParams,
-		       const double* dbValues,
-		       const int numValues) const;
-
-  /** Compute density at locations from parameters.
-   *
-   * Results are stored in _density.
-   *
-   * Index into parameters = iLoc*numParameters + iParam
-   *
-   * @param parameters Parameters at location
-   * @param numParameters Number of parameters
-   * @param numLocs Number of locations
-   */
-  void _calcDensity(const double* parameters,
-		    const int numParameters,
-		    const int numLocs);
-
-  /** Compute stress tensor at locations from parameters.
-   *
-   * Results are stored in _stress.
-   *
-   * Index into parameters = iLoc*numParameters + iParam
-   *
-   * @param parameters Parameters at location
-   * @param numParameters Number of parameters
-   * @param totalStrain Total strain at locations
-   * @param numLocs Number of locations
-   * @param spaceDim Spatial dimension for locations.
-   */
-  void _calcStress(const double* parameters,
-		   const int numParameters,
-		   const double* totalStrain,
-		   const int numLocs,
-		   const int spaceDim);
-
-  /** Compute derivatives of elasticity matrix at locations from parameters.
-   *
-   * Results are stored in _elasticConsts.
-   *
-   * @param parameters Parameters at locations.
-   * @param numParameters Number of parameters.
-   * @param totalStrain Total strain at locations.
-   * @param numLocs Number of locations.
-   * @param spaceDim Spatial dimension for locations.
-   */
-  void _calcElasticConsts(const double* parameters,
-			  const int numParameters,
-			  const double* totalStrain,
-			  const int numLocs,
-			  const int spaceDim);
-
-  // NOT IMPLEMENTED ////////////////////////////////////////////////////
-private :
-
-  /// Not implemented
-  const ElasticIsotropic1D& operator=(const ElasticIsotropic1D& m);
-
-}; // class ElasticIsotropic1D
-
-#include "ElasticIsotropic1D.icc" // inline methods
-
-#endif // pylith_materials_elasticisotropic1d_hh
-
-
-// End of file 

Deleted: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.icc	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.icc	2007-04-08 02:00:47 UTC (rev 6521)
@@ -1,25 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-//                           Brad T. Aagaard
-//                        U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-#if !defined(pylith_materials_elasticisotropic1d_hh)
-#error "ElasticIsotropic1D.icc can only be included from ElasticIsotropic1D.hh"
-#endif
-
-// Create a pointer to a copy of this.
-inline
-pylith::materials::ElasticMaterial*
-pylith::materials::ElasticIsotropic1D::clone(void) const {
-  return new ElasticIsotropic1D(*this);
-}
-
-
-// End of file 

Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc	2007-04-08 02:00:47 UTC (rev 6521)
@@ -17,37 +17,38 @@
 #include <assert.h> // USES assert()
 
 // ----------------------------------------------------------------------
-// _ElasticIsotropic is a helper class for ElasticIsotropic. We define
-// it in this implementation file to insulate other objects from these
-// details.
 namespace pylith {
   namespace materials {
     class _ElasticIsotropic3D;
   } // materials
 } // pylith
 
+/** _ElasticIsotropic is a helper class for ElasticIsotropic. We define
+ * it in this implementation file to insulate other objects from these
+ * details.
+ */
 class pylith::materials::_ElasticIsotropic3D {
 public:
-  // Number of entries in stress tensor.
+  /// Number of entries in stress tensor.
   static const int stressSize;
 
-  // Number of entries in derivative of elasticity matrix.
+  /// Number of entries in derivative of elasticity matrix.
   static const int numElasticConsts;
 
-  // Values expected in spatial database
+  /// Values expected in spatial database
   static const int numDBValues;
   static const char* namesDBValues[];
 
-  // Indices (order) of database values
+  /// Indices (order) of database values
   static const int didDensity;
   static const int didVs;
   static const int didVp;
 
-  // Parameters
+  /// Parameters
   static const int numParameters;
   static const char* namesParameters[];
 
-  // Indices (order) of parameters
+  /// Indices (order) of parameters
   static const int pidDensity;
   static const int pidMu;
   static const int pidLambda;
@@ -165,108 +166,100 @@
 // ----------------------------------------------------------------------
 // Compute density at location from parameters.
 void
-pylith::materials::ElasticIsotropic3D::_calcDensity(const double* parameters,
-						    const int numParameters,
-						    const int numLocs)
+pylith::materials::ElasticIsotropic3D::_calcDensity(double* const density,
+						    const int size,
+						    const double* parameters,
+						    const int numParameters)
 { // _calcDensity
-  assert(0 != _density);
+  assert(0 != density);
+  assert(1 == size);
   assert(0 != parameters);
 
-  for (int iLoc=0, index=0; iLoc < numLocs; ++iLoc, index+=numParameters)
-    _density[iLoc] = 
-      parameters[index+_ElasticIsotropic3D::pidDensity];
+  *density = parameters[_ElasticIsotropic3D::pidDensity];
 } // _calcDensity
 
 // ----------------------------------------------------------------------
 // Compute stress tensor at location from parameters.
 void
-pylith::materials::ElasticIsotropic3D::_calcStress(const double* parameters,
+pylith::materials::ElasticIsotropic3D::_calcStress(double* const stress,
+						   const int size,
+						   const double* parameters,
 						   const int numParameters,
 						   const double* totalStrain,
-						   const int numLocs,
 						   const int spaceDim)
 { // _calcStress
-  assert(0 != _stress);
+  assert(0 != stress);
+  assert(_ElasticIsotropic3D::stressSize == size);
   assert(0 != parameters);
   assert(_ElasticIsotropic3D::numParameters == numParameters);
   assert(spaceDim == _dimension);
 
-  for (int iLoc=0, indexP=0, indexS=0;
-       iLoc < numLocs; 
-       ++iLoc, 
-	 indexP+=_ElasticIsotropic3D::numParameters,
-	 indexS+=_ElasticIsotropic3D::stressSize) {
-    const double density = parameters[indexP+_ElasticIsotropic3D::pidDensity];
-    const double mu = parameters[indexP+_ElasticIsotropic3D::pidMu];
-    const double lambda = parameters[indexP+_ElasticIsotropic3D::pidLambda];
+  const double density = parameters[_ElasticIsotropic3D::pidDensity];
+  const double mu = parameters[_ElasticIsotropic3D::pidMu];
+  const double lambda = parameters[_ElasticIsotropic3D::pidLambda];
 
-    const double lambda2mu = lambda + 2.0 * mu;
+  const double lambda2mu = lambda + 2.0 * mu;
   
-    const double e11 = totalStrain[indexS  ];
-    const double e22 = totalStrain[indexS+1];
-    const double e33 = totalStrain[indexS+2];
-    const double e12 = totalStrain[indexS+3];
-    const double e23 = totalStrain[indexS+4];
-    const double e13 = totalStrain[indexS+5];
+  const double e11 = totalStrain[0];
+  const double e22 = totalStrain[1];
+  const double e33 = totalStrain[2];
+  const double e12 = totalStrain[3];
+  const double e23 = totalStrain[4];
+  const double e13 = totalStrain[5];
+  
+  const double s123 = lambda * (e11 + e22 + e33);
 
-    const double s123 = lambda * (e11 + e22 + e33);
-
-    _stress[indexS  ] = s123 + 2.0*mu*e11;
-    _stress[indexS+1] = s123 + 2.0*mu*e22;
-    _stress[indexS+2] = s123 + 2.0*mu*e33;
-    _stress[indexS+3] = 2.0 * mu * e12;
-    _stress[indexS+4] = 2.0 * mu * e23;
-    _stress[indexS+5] = 2.0 * mu * e13;
-  } // for
+  stress[0] = s123 + 2.0*mu*e11;
+  stress[1] = s123 + 2.0*mu*e22;
+  stress[2] = s123 + 2.0*mu*e33;
+  stress[3] = 2.0 * mu * e12;
+  stress[4] = 2.0 * mu * e23;
+  stress[5] = 2.0 * mu * e13;
 } // _calcStress
 
 // ----------------------------------------------------------------------
 // Compute derivative of elasticity matrix at location from parameters.
 void
-pylith::materials::ElasticIsotropic3D::_calcElasticConsts(const double* parameters,
+pylith::materials::ElasticIsotropic3D::_calcElasticConsts(double* const elasticConsts,
+							  const int size,
+							  const double* parameters,
 							  const int numParameters,
 							  const double* totalStrain,
-							  const int numLocs,
 							  const int spaceDim)
 { // _calcElasticConsts
-  assert(0 != _elasticConsts);
+  assert(0 != elasticConsts);
+  assert(_ElasticIsotropic3D::numElasticConsts == size);
   assert(0 != parameters);
   assert(_ElasticIsotropic3D::numParameters == numParameters);
   assert(spaceDim == _dimension);
  
-  for (int iLoc=0, indexP=0, indexC=0;
-       iLoc < numLocs; 
-       ++iLoc, 
-       indexP+=_ElasticIsotropic3D::numParameters,
-       indexC+=_ElasticIsotropic3D::numElasticConsts) {
-    const double density = parameters[indexP+_ElasticIsotropic3D::pidDensity];
-    const double mu = parameters[indexP+_ElasticIsotropic3D::pidMu];
-    const double lambda = parameters[indexP+_ElasticIsotropic3D::pidLambda];
+  const double density = parameters[_ElasticIsotropic3D::pidDensity];
+  const double mu = parameters[_ElasticIsotropic3D::pidMu];
+  const double lambda = parameters[_ElasticIsotropic3D::pidLambda];
 
-    const double lambda2mu = lambda + 2.0 * mu;
+  const double lambda2mu = lambda + 2.0 * mu;
    
-    _elasticConsts[indexC+ 0] = lambda2mu; // C1111
-    _elasticConsts[indexC+ 1] = lambda; // C1122
-    _elasticConsts[indexC+ 2] = lambda; // C1133
-    _elasticConsts[indexC+ 3] = 0; // C1112
-    _elasticConsts[indexC+ 4] = 0; // C1123
-    _elasticConsts[indexC+ 5] = 0; // C1113
-    _elasticConsts[indexC+ 6] = lambda2mu; // C2222
-    _elasticConsts[indexC+ 7] = lambda; // C2233
-    _elasticConsts[indexC+ 8] = 0; // C2212
-    _elasticConsts[indexC+ 9] = 0; // C2223
-    _elasticConsts[indexC+10] = 0; // C2213
-    _elasticConsts[indexC+11] = lambda2mu; // C3333
-    _elasticConsts[indexC+12] = 0; // C3312
-    _elasticConsts[indexC+13] = 0; // C3323
-    _elasticConsts[indexC+14] = 0; // C3313
-    _elasticConsts[indexC+15] = 2.0 * mu; // C1212
-    _elasticConsts[indexC+16] = 0; // C1223
-    _elasticConsts[indexC+17] = 0; // C1213
-    _elasticConsts[indexC+18] = 2.0 * mu; // C2323
-    _elasticConsts[indexC+19] = 0; // C2313
-    _elasticConsts[indexC+20] = 2.0 * mu; // C1313
-  } // for
+  elasticConsts[ 0] = lambda2mu; // C1111
+  elasticConsts[ 1] = lambda; // C1122
+  elasticConsts[ 2] = lambda; // C1133
+  elasticConsts[ 3] = 0; // C1112
+  elasticConsts[ 4] = 0; // C1123
+  elasticConsts[ 5] = 0; // C1113
+  elasticConsts[ 6] = lambda2mu; // C2222
+  elasticConsts[ 7] = lambda; // C2233
+  elasticConsts[ 8] = 0; // C2212
+  elasticConsts[ 9] = 0; // C2223
+  elasticConsts[10] = 0; // C2213
+  elasticConsts[11] = lambda2mu; // C3333
+  elasticConsts[12] = 0; // C3312
+  elasticConsts[13] = 0; // C3323
+  elasticConsts[14] = 0; // C3313
+  elasticConsts[15] = 2.0 * mu; // C1212
+  elasticConsts[16] = 0; // C1223
+  elasticConsts[17] = 0; // C1213
+  elasticConsts[18] = 2.0 * mu; // C2323
+  elasticConsts[19] = 0; // C2313
+  elasticConsts[20] = 2.0 * mu; // C1313
 } // _calcElasticConsts
 
 

Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh	2007-04-08 02:00:47 UTC (rev 6521)
@@ -125,54 +125,51 @@
 		       const double* dbValues,
 		       const int numValues) const;
 
-  /** Compute density at locations from parameters.
+  /** Compute density from parameters.
    *
-   * Results are stored in _density.
-   *
-   * Index into parameters = iLoc*numParameters + iParam
-   *
+   * @param density Array for density
+   * @param size Size of array for density
    * @param parameters Parameters at location
    * @param numParameters Number of parameters
-   * @param numLocs Number of locations
    */
-  void _calcDensity(const double* parameters,
-		    const int numParameters,
-		    const int numLocs);
+  void _calcDensity(double* const density,
+		    const int size,
+		    const double* parameters,
+		    const int numParameters);
 
-  /** Compute stress tensor at locations from parameters.
+  /** Compute stress tensor from parameters.
    *
-   * Results are stored in _stress.
-   *
-   * Index into parameters = iLoc*numParameters + iParam
-   *
-   * @param parameters Parameters at location
-   * @param numParameters Number of parameters
-   * @param totalStrain Total strain at locations
-   * @param numLocs Number of locations
+   * @param stress Array for stress tensor
+   * @param size Size of array for stress tensor
+   * @param parameters Parameters at locations.
+   * @param numParameters Number of parameters.
+   * @param totalStrain Total strain at locations.
    * @param spaceDim Spatial dimension for locations.
    */
-  void _calcStress(const double* parameters,
+  void _calcStress(double* const stress,
+		   const int size,
+		   const double* parameters,
 		   const int numParameters,
 		   const double* totalStrain,
-		   const int numLocs,
 		   const int spaceDim);
 
-  /** Compute derivatives of elasticity matrix at locations from parameters.
+  /** Compute derivatives of elasticity matrix from parameters.
    *
-   * Results are stored in _elasticConsts.
-   *
+   * @param elasticConsts Array for elastic constants
+   * @param size Size of array
    * @param parameters Parameters at locations.
    * @param numParameters Number of parameters.
    * @param totalStrain Total strain at locations.
-   * @param numLocs Number of locations.
    * @param spaceDim Spatial dimension for locations.
    */
-  void _calcElasticConsts(const double* parameters,
+  void _calcElasticConsts(double* const elasticConsts,
+			  const int size,
+			  const double* parameters,
 			  const int numParameters,
 			  const double* totalStrain,
-			  const int numLocs,
 			  const int spaceDim);
 
+
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :
 

Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc	2007-04-08 02:00:47 UTC (rev 6521)
@@ -62,7 +62,8 @@
   double* paramsCell = 0;
   _getParameters(&paramsCell, cell, numQuadPts);
   const int numParams = _numParameters();
-  _calcDensity(paramsCell, numParams, numQuadPts);
+  for (int iQuad=0, indexP=0; iQuad < numQuadPts; ++iQuad, indexP+=numParams)
+    _calcDensity(&_density[iQuad], 1, &paramsCell[indexP], numParams);
   delete[] paramsCell; paramsCell = 0;
 
   return _density;
@@ -71,11 +72,10 @@
 // ----------------------------------------------------------------------
 // Compute stress tensor for cell at quadrature points.
 const double*
-pylith::materials::ElasticMaterial::calcStress(
-				     const Mesh::point_type& cell,
-				     const double* totalStrain,
-				     const int numQuadPts,
-				     const int spaceDim)
+pylith::materials::ElasticMaterial::calcStress(const Mesh::point_type& cell,
+					       const double* totalStrain,
+					       const int numQuadPts,
+					       const int spaceDim)
 { // calcStress
   assert(0 != totalStrain);
   assert(0 != _stress);
@@ -86,7 +86,13 @@
   double* paramsCell = 0;
   _getParameters(&paramsCell, cell, numQuadPts);
   const int numParams = _numParameters();
-  _calcStress(paramsCell, numParams, totalStrain, numQuadPts, spaceDim);
+  const int stressOffset = stressSize();
+  for (int iQuad=0, iStress=0, indexP=0; 
+       iQuad < numQuadPts;
+       ++iQuad, iStress+=stressOffset, indexP+=numParams)
+    _calcStress(&_stress[iStress], stressOffset, 
+		&paramsCell[indexP], numParams, 
+		&totalStrain[iStress], spaceDim);
   delete[] paramsCell; paramsCell = 0;
 
   return _stress;
@@ -110,7 +116,14 @@
   double* paramsCell = 0;
   _getParameters(&paramsCell, cell, numQuadPts);
   const int numParams = _numParameters();
-  _calcElasticConsts(paramsCell, numParams, totalStrain, numQuadPts, spaceDim);
+  const int stressOffset = stressSize();
+  const int constOffset = numElasticConsts();
+  for (int iQuad=0, iStress=0, iConst=0, indexP=0;
+       iQuad < numQuadPts;
+       ++iQuad, iStress+=stressOffset, iConst += constOffset, indexP+=numParams)
+    _calcElasticConsts(&_elasticConsts[iConst], constOffset, 
+		       &paramsCell[indexP], numParams, 
+		       &totalStrain[iStress], spaceDim);
   delete[] paramsCell; paramsCell = 0;
 
   return _elasticConsts;

Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh	2007-04-08 02:00:47 UTC (rev 6521)
@@ -41,7 +41,7 @@
   // PUBLIC TYPEDEFS ////////////////////////////////////////////////////
 public :
 
-  typedef ALE::Mesh               Mesh;
+  typedef ALE::Mesh Mesh;
   typedef Mesh::real_section_type real_section_type;
 
   // PUBLIC METHODS /////////////////////////////////////////////////////
@@ -168,51 +168,51 @@
    */
   ElasticMaterial(const ElasticMaterial& m);
 
-  /** Compute density at locations from parameters.
+  /** Compute density from parameters.
    *
-   * Results are stored in _density.
-   *
+   * @param density Array for density
+   * @param size Size of array for density
    * @param parameters Parameters at location
    * @param numParameters Number of parameters
-   * @param numLocs Number of locations
    */
   virtual
-  void _calcDensity(const double* parameters,
-		    const int numParameters,
-		    const int numLocs) = 0;
+  void _calcDensity(double* const density,
+		    const int size,
+		    const double* parameters,
+		    const int numParameters) = 0;
 
-  /** Compute stress density at locations from parameters.
+  /** Compute stress tensor from parameters.
    *
-   * Results are stored in _stress.
-   *
+   * @param stress Array for stress tensor
+   * @param size Size of array for stress tensor
    * @param parameters Parameters at locations.
    * @param numParameters Number of parameters.
    * @param totalStrain Total strain at locations.
-   * @param numLocs Number of locations.
    * @param spaceDim Spatial dimension for locations.
    */
   virtual
-  void _calcStress(const double* parameters,
+  void _calcStress(double* const stress,
+		   const int size,
+		   const double* parameters,
 		   const int numParameters,
 		   const double* totalStrain,
-		   const int numLocs,
 		   const int spaceDim) = 0;
 
-  /** Compute derivatives of elasticity matrix at locations from parameters.
+  /** Compute derivatives of elasticity matrix from parameters.
    *
-   * Results are stored in _elasticConsts.
-   *
+   * @param elasticConsts Array for elastic constants
+   * @param size Size of array
    * @param parameters Parameters at locations.
    * @param numParameters Number of parameters.
    * @param totalStrain Total strain at locations.
-   * @param numLocs Number of locations.
    * @param spaceDim Spatial dimension for locations.
    */
   virtual
-  void _calcElasticConsts(const double* parameters,
+  void _calcElasticConsts(double* const elasticConsts,
+			  const int size,
+			  const double* parameters,
 			  const int numParameters,
 			  const double* totalStrain,
-			  const int numLocs,
 			  const int spaceDim) = 0;
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
@@ -221,9 +221,27 @@
   /// Not implemented
   const ElasticMaterial& operator=(const ElasticMaterial& m);
 
-  // PROTECTED MEMBERS //////////////////////////////////////////////////
-protected :
+  // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
 
+  /** Get parameters for cell.
+   *
+   * Parameters are returned in paramsCells.
+   *
+   * size = numQuadPts * numParams
+   * index = iQuad*numParams + iParam
+   *
+   * @param paramsCell Array of parameters for cell
+   * @param cell Finite-element cell
+   * @param numQuadPts Number of quadrature points (consistency check)
+   */
+  void _getParameters(double** paramsCells,
+		      const Mesh::point_type& cell,
+		      const int numQuadPts);
+
+  // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
   /** Density value at quadrature points for current cell.
    *
    * size = numQuadPts
@@ -231,33 +249,20 @@
    */
   double* _density;
 
-  /** Array of stress tensor at quadrature points for current cell.
+  /** Stress tensor at quadrature points for current cell.
    *
    * size = stressSize*numQuadPts
    * index = iQuadPt*stressSize+iStress
    */
   double* _stress;
 
-  /** Array of elasticity constants at quadrature points for current cell.
+  /** Elasticity matrix at quadrature points for current cell.
    *
    * size = numElasticConsts*numQuadPts
    * index = iQuadPt*numElasticConsts+iConstant
    */
   double* _elasticConsts;
 
-  // PRIVATE METHODS ////////////////////////////////////////////////////
-private :
-
-  /** Get parameters for cell.
-   *
-   * @param paramsCell Array of parameters for cell
-   * @param cell Finite-element cell
-   * @param numQuadPts Number of quadrature points (consistency check)
-   */
-  void _getParameters(double** paramsCells,
-		      const Mesh::point_type& cell,
-		      const int numQuadPts);
-
 }; // class ElasticMaterial
 
 #endif // pylith_materials_elasticmaterial_hh

Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc	2007-04-08 02:00:47 UTC (rev 6521)
@@ -17,37 +17,38 @@
 #include <assert.h> // USES assert()
 
 // ----------------------------------------------------------------------
-// _ElasticPlaneStrain is a helper class for ElasticPlaneStrain. We define
-// it in this implementation file to insulate other objects from these
-// details.
 namespace pylith {
   namespace materials {
     class _ElasticPlaneStrain;
   } // materials
 } // pylith
 
+/** _ElasticPlaneStrain is a helper class for ElasticPlaneStrain. We define
+ * it in this implementation file to insulate other objects from these
+ * details.
+ */
 class pylith::materials::_ElasticPlaneStrain {
 public:
-  // Number of entries in stress tensor.
+  /// Number of entries in stress tensor.
   static const int stressSize;
 
-  // Number of elastic constants (for general 3-D elastic material)
+  /// Number of elastic constants (for general 3-D elastic material)
   static const int numElasticConsts;
 
-  // Values expected in spatial database
+  /// Values expected in spatial database
   static const int numDBValues;
   static const char* namesDBValues[];
 
-  // Indices (order) of database values
+  /// Indices (order) of database values
   static const int didDensity;
   static const int didVs;
   static const int didVp;
 
-  // Parameters
+  /// Parameters
   static const int numParameters;
   static const char* namesParameters[];
 
-  // Indices (order) of parameters
+  /// Indices (order) of parameters
   static const int pidDensity;
   static const int pidMu;
   static const int pidLambda;
@@ -165,86 +166,78 @@
 // ----------------------------------------------------------------------
 // Compute density at location from parameters.
 void
-pylith::materials::ElasticPlaneStrain::_calcDensity(const double* parameters,
-						    const int numParameters,
-						    const int numLocs)
+pylith::materials::ElasticPlaneStrain::_calcDensity(double* const density,
+						    const int size,
+						    const double* parameters,
+						    const int numParameters)
 { // calcDensity
-  assert(0 != _density);
+  assert(0 != density);
+  assert(1 == size);
   assert(0 != parameters);
 
-  for (int iLoc=0, index=0; iLoc < numLocs; ++iLoc, index+=numParameters)
-    _density[iLoc] = 
-      parameters[index+_ElasticPlaneStrain::pidDensity];
+  *density = parameters[_ElasticPlaneStrain::pidDensity];
 } // calcDensity
 
 // ----------------------------------------------------------------------
 // Compute stress tensor at location from parameters.
 void
-pylith::materials::ElasticPlaneStrain::_calcStress(const double* parameters,
+pylith::materials::ElasticPlaneStrain::_calcStress(double* const stress,
+						   const int size,
+						   const double* parameters,
 						   const int numParameters,
 						   const double* totalStrain,
-						   const int numLocs,
 						   const int spaceDim)
 { // _calcStress
-  assert(0 != _stress);
+  assert(0 != stress);
+  assert(_ElasticPlaneStrain::stressSize == size);
   assert(0 != parameters);
   assert(_ElasticPlaneStrain::numParameters == numParameters);
   assert(spaceDim == _dimension);
 
-  for (int iLoc=0, indexP=0, indexS=0;
-       iLoc < numLocs; 
-       ++iLoc, 
-	 indexP+=_ElasticPlaneStrain::numParameters,
-	 indexS+=_ElasticPlaneStrain::stressSize) {
-    const double density = parameters[indexP+_ElasticPlaneStrain::pidDensity];
-    const double mu = parameters[indexP+_ElasticPlaneStrain::pidMu];
-    const double lambda = parameters[indexP+_ElasticPlaneStrain::pidLambda];
-
-    const double lambda2mu = lambda + 2.0 * mu;
+  const double density = parameters[_ElasticPlaneStrain::pidDensity];
+  const double mu = parameters[_ElasticPlaneStrain::pidMu];
+  const double lambda = parameters[_ElasticPlaneStrain::pidLambda];
   
-    const double e11 = totalStrain[indexS  ];
-    const double e22 = totalStrain[indexS+1];
-    const double e12 = totalStrain[indexS+3];
+  const double lambda2mu = lambda + 2.0 * mu;
+  
+  const double e11 = totalStrain[0];
+  const double e22 = totalStrain[1];
+  const double e12 = totalStrain[2];
 
-    const double s12 = lambda * (e11 + e22);
+  const double s12 = lambda * (e11 + e22);
 
-    _stress[indexS  ] = s12 + 2.0*mu*e11;
-    _stress[indexS+1] = s12 + 2.0*mu*e22;
-    _stress[indexS+2] = 2.0 * mu * e12;
-  } // for
+  stress[0] = s12 + 2.0*mu*e11;
+  stress[1] = s12 + 2.0*mu*e22;
+  stress[2] = 2.0 * mu * e12;
 } // _calcStress
 
 // ----------------------------------------------------------------------
 // Compute density at location from parameters.
 void
-pylith::materials::ElasticPlaneStrain::_calcElasticConsts(const double* parameters,
+pylith::materials::ElasticPlaneStrain::_calcElasticConsts(double* const elasticConsts,
+							  const int size,
+							  const double* parameters,
 							  const int numParameters,
 							  const double* totalStrain,
-							  const int numLocs,
 							  const int spaceDim)
 { // calcElasticConsts
-  assert(0 != _elasticConsts);
+  assert(0 != elasticConsts);
+  assert(_ElasticPlaneStrain::numElasticConsts == size);
   assert(0 != parameters);
   assert(_ElasticPlaneStrain::numParameters == numParameters);
 
-  for (int iLoc=0, indexP=0, indexC=0;
-       iLoc < numLocs; 
-       ++iLoc, 
-	 indexP+=_ElasticPlaneStrain::numParameters,
-	 indexC+=_ElasticPlaneStrain::numElasticConsts) {
-    const double density = parameters[indexP+_ElasticPlaneStrain::pidDensity];
-    const double mu = parameters[indexP+_ElasticPlaneStrain::pidMu];
-    const double lambda = parameters[indexP+_ElasticPlaneStrain::pidLambda];
+  const double density = parameters[_ElasticPlaneStrain::pidDensity];
+  const double mu = parameters[_ElasticPlaneStrain::pidMu];
+  const double lambda = parameters[_ElasticPlaneStrain::pidLambda];
 
-    const double lambda2mu = lambda + 2.0*mu;
+  const double lambda2mu = lambda + 2.0*mu;
   
-    _elasticConsts[indexC+ 0] = lambda2mu; // C1111
-    _elasticConsts[indexC+ 1] = lambda; // C1122
-    _elasticConsts[indexC+ 3] = 0; // C1112
-    _elasticConsts[indexC+ 6] = lambda2mu; // C2222
-    _elasticConsts[indexC+ 8] = 0; // C2212
-    _elasticConsts[indexC+15] = 2.0*mu; // C1212
-  } // for
+  elasticConsts[0] = lambda2mu; // C1111
+  elasticConsts[1] = lambda; // C1122
+  elasticConsts[2] = 0; // C1112
+  elasticConsts[3] = lambda2mu; // C2222
+  elasticConsts[4] = 0; // C2212
+  elasticConsts[5] = 2.0*mu; // C1212
 } // calcElasticConsts
 
 

Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh	2007-04-08 02:00:47 UTC (rev 6521)
@@ -14,7 +14,7 @@
  *
  * @brief C++ ElasticPlaneStrain object
  *
- * 3-D, isotropic, linear elastic material for plane strain. The
+ * 2-D, isotropic, linear elastic material for plane strain. The
  * physical properties are specified using density, shear-wave speed,
  * and compressional-wave speed. The physical properties are stored
  * internally using density, lambda, and mu, which are directly
@@ -125,52 +125,48 @@
 		       const double* dbValues,
 		       const int numValues) const;
 
-  /** Compute density at locations from parameters.
+  /** Compute density from parameters.
    *
-   * Results are stored in _density.
-   *
-   * Index into parameters = iLoc*numParameters + iParam
-   *
+   * @param density Array for density
+   * @param size Size of array for density
    * @param parameters Parameters at location
    * @param numParameters Number of parameters
-   * @param numLocs Number of locations
    */
-  void _calcDensity(const double* parameters,
-		    const int numParameters,
-		    const int numLocs);
+  void _calcDensity(double* const density,
+		    const int size,
+		    const double* parameters,
+		    const int numParameters);
 
-  /** Compute stress tensor at locations from parameters.
+  /** Compute stress tensor from parameters.
    *
-   * Results are stored in _stress.
-   *
-   * Index into parameters = iLoc*numParameters + iParam
-   *
-   * @param parameters Parameters at location
-   * @param numParameters Number of parameters
-   * @param totalStrain Total strain at locations
-   * @param numLocs Number of locations
+   * @param stress Array for stress tensor
+   * @param size Size of array for stress tensor
+   * @param parameters Parameters at locations.
+   * @param numParameters Number of parameters.
+   * @param totalStrain Total strain at locations.
    * @param spaceDim Spatial dimension for locations.
    */
-  void _calcStress(const double* parameters,
+  void _calcStress(double* const stress,
+		   const int size,
+		   const double* parameters,
 		   const int numParameters,
 		   const double* totalStrain,
-		   const int numLocs,
 		   const int spaceDim);
 
-  /** Compute derivatives of elasticity matrix at locations from parameters.
+  /** Compute derivatives of elasticity matrix from parameters.
    *
-   * Results are stored in _elasticConsts.
-   *
+   * @param elasticConsts Array for elastic constants
+   * @param size Size of array
    * @param parameters Parameters at locations.
    * @param numParameters Number of parameters.
    * @param totalStrain Total strain at locations.
-   * @param numLocs Number of locations.
    * @param spaceDim Spatial dimension for locations.
    */
-  void _calcElasticConsts(const double* parameters,
+  void _calcElasticConsts(double* const elasticConsts,
+			  const int size,
+			  const double* parameters,
 			  const int numParameters,
 			  const double* totalStrain,
-			  const int numLocs,
 			  const int spaceDim);
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////

Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc	2007-04-08 02:00:47 UTC (rev 6521)
@@ -17,37 +17,38 @@
 #include <assert.h> // USES assert()
 
 // ----------------------------------------------------------------------
-// _ElasticPlaneStress is a helper class for ElasticPlaneStress. We define
-// it in this implementation file to insulate other objects from these
-// details.
 namespace pylith {
   namespace materials {
     class _ElasticPlaneStress;
   } // materials
 } // pylith
 
+/** _ElasticPlaneStress is a helper class for ElasticPlaneStress. We define
+ * it in this implementation file to insulate other objects from these
+ * details.
+ */
 class pylith::materials::_ElasticPlaneStress {
 public:
-  // Number of entries in stress tensor.
+  /// Number of entries in stress tensor.
   static const int stressSize;
 
-  // Number of elastic constants (for general 3-D elastic material)
+  /// Number of elastic constants (for general 3-D elastic material)
   static const int numElasticConsts;
 
-  // Values expected in spatial database
+  /// Values expected in spatial database
   static const int numDBValues;
   static const char* namesDBValues[];
 
-  // Indices (order) of database values
+  /// Indices (order) of database values
   static const int didDensity;
   static const int didVs;
   static const int didVp;
 
-  // Parameters
+  /// Parameters
   static const int numParameters;
   static const char* namesParameters[];
 
-  // Indices (order) of parameters
+  /// Indices (order) of parameters
   static const int pidDensity;
   static const int pidMu;
   static const int pidLambda;
@@ -165,88 +166,78 @@
 // ----------------------------------------------------------------------
 // Compute density at location from parameters.
 void
-pylith::materials::ElasticPlaneStress::_calcDensity(const double* parameters,
-						    const int numParameters,
-						    const int numLocs)
+pylith::materials::ElasticPlaneStress::_calcDensity(double* const density,
+						    const int size,
+						    const double* parameters,
+						    const int numParameters)
 { // calcDensity
-  assert(0 != _density);
+  assert(0 != density);
+  assert(1 == size);
   assert(0 != parameters);
 
-  for (int iLoc=0, index=0; iLoc < numLocs; ++iLoc, index+=numParameters)
-    _density[iLoc] = 
-      parameters[index+_ElasticPlaneStress::pidDensity];
+  *density = parameters[_ElasticPlaneStress::pidDensity];
 } // calcDensity
 
 // ----------------------------------------------------------------------
 // Compute stress tensor at location from parameters.
 void
-pylith::materials::ElasticPlaneStress::_calcStress(const double* parameters,
+pylith::materials::ElasticPlaneStress::_calcStress(double* const stress,
+						   const int size,
+						   const double* parameters,
 						   const int numParameters,
 						   const double* totalStrain,
-						   const int numLocs,
 						   const int spaceDim)
 { // _calcStress
-  assert(0 != _stress);
+  assert(0 != stress);
+  assert(_ElasticPlaneStress::stressSize == size);
   assert(0 != parameters);
   assert(_ElasticPlaneStress::numParameters == numParameters);
   assert(spaceDim == _dimension);
 
-  for (int iLoc=0, indexP=0, indexS=0;
-       iLoc < numLocs; 
-       ++iLoc, 
-	 indexP+=_ElasticPlaneStress::numParameters,
-	 indexS+=_ElasticPlaneStress::stressSize) {
-    const double density = parameters[indexP+_ElasticPlaneStress::pidDensity];
-    const double mu = parameters[indexP+_ElasticPlaneStress::pidMu];
-    const double lambda = parameters[indexP+_ElasticPlaneStress::pidLambda];
+  const double density = parameters[_ElasticPlaneStress::pidDensity];
+  const double mu = parameters[_ElasticPlaneStress::pidMu];
+  const double lambda = parameters[_ElasticPlaneStress::pidLambda];
 
-    const double lambda2mu = lambda + 2.0 * mu;
-    const double lambdamu = lambda + mu;
+  const double lambda2mu = lambda + 2.0 * mu;
+  const double lambdamu = lambda + mu;
   
-    const double e11 = totalStrain[indexS  ];
-    const double e22 = totalStrain[indexS+1];
-    const double e12 = totalStrain[indexS+3];
+  const double e11 = totalStrain[0];
+  const double e22 = totalStrain[1];
+  const double e12 = totalStrain[2];
 
-    _stress[indexS  ] = 
-      (4.0*mu*lambdamu * e11 + 2.0*mu*lambda * e22)/lambda2mu;
-    _stress[indexS+1] =
-      (2.0*mu*lambda * e11 + 4.0*mu*lambdamu * e22)/lambda2mu;
-    _stress[indexS+2] = 2.0 * mu * e12;
-  } // for
+  stress[0] = (4.0*mu*lambdamu * e11 + 2.0*mu*lambda * e22)/lambda2mu;
+  stress[1] = (2.0*mu*lambda * e11 + 4.0*mu*lambdamu * e22)/lambda2mu;
+  stress[2] = 2.0 * mu * e12;
 } // _calcStress
 
 // ----------------------------------------------------------------------
 // Compute density at location from parameters.
 void
-pylith::materials::ElasticPlaneStress::_calcElasticConsts(const double* parameters,
+pylith::materials::ElasticPlaneStress::_calcElasticConsts(double* const elasticConsts,
+							  const int size,
+							  const double* parameters,
 							  const int numParameters,
 							  const double* totalStrain,
-							  const int numLocs,
 							  const int spaceDim)
 { // calcElasticConsts
-  assert(0 != _elasticConsts);
+  assert(0 != elasticConsts);
+  assert(_ElasticPlaneStress::numElasticConsts == size);
   assert(0 != parameters);
   assert(_ElasticPlaneStress::numParameters == numParameters);
 
-  for (int iLoc=0, indexP=0, indexC=0;
-       iLoc < numLocs; 
-       ++iLoc, 
-	 indexP+=_ElasticPlaneStress::numParameters,
-	 indexC+=_ElasticPlaneStress::numElasticConsts) {
-    const double density = parameters[indexP+_ElasticPlaneStress::pidDensity];
-    const double mu = parameters[indexP+_ElasticPlaneStress::pidMu];
-    const double lambda = parameters[indexP+_ElasticPlaneStress::pidLambda];
+  const double density = parameters[_ElasticPlaneStress::pidDensity];
+  const double mu = parameters[_ElasticPlaneStress::pidMu];
+  const double lambda = parameters[_ElasticPlaneStress::pidLambda];
   
-    const double lambda2mu = lambda + 2.0 * mu;
-    const double c11 = 4.0 * mu * (lambda + mu) / lambda2mu;
+  const double lambda2mu = lambda + 2.0 * mu;
+  const double c11 = 4.0 * mu * (lambda + mu) / lambda2mu;
 
-    _elasticConsts[indexC+0] = c11; // C1111
-    _elasticConsts[indexC+1] = 2.0 * mu * lambda / lambda2mu; // C1122
-    _elasticConsts[indexC+2] = 0; // C1112
-    _elasticConsts[indexC+3] = c11; // C2222
-    _elasticConsts[indexC+4] = 0; // C2212
-    _elasticConsts[indexC+5] = 2.0 * mu; // C1212
-  } // for
+  elasticConsts[0] = c11; // C1111
+  elasticConsts[1] = 2.0 * mu * lambda / lambda2mu; // C1122
+  elasticConsts[2] = 0; // C1112
+  elasticConsts[3] = c11; // C2222
+  elasticConsts[4] = 0; // C2212
+  elasticConsts[5] = 2.0 * mu; // C1212
 } // calcElasticConsts
 
 

Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh	2007-04-08 02:00:47 UTC (rev 6521)
@@ -125,52 +125,48 @@
 		       const double* dbValues,
 		       const int numValues) const;
 
-  /** Compute density at locations from parameters.
+  /** Compute density from parameters.
    *
-   * Results are stored in _density.
-   *
-   * Index into parameters = iLoc*numParameters + iParam
-   *
+   * @param density Array for density
+   * @param size Size of array for density
    * @param parameters Parameters at location
    * @param numParameters Number of parameters
-   * @param numLocs Number of locations
    */
-  void _calcDensity(const double* parameters,
-		    const int numParameters,
-		    const int numLocs);
+  void _calcDensity(double* const density,
+		    const int size,
+		    const double* parameters,
+		    const int numParameters);
 
-  /** Compute stress tensor at locations from parameters.
+  /** Compute stress tensor from parameters.
    *
-   * Results are stored in _stress.
-   *
-   * Index into parameters = iLoc*numParameters + iParam
-   *
-   * @param parameters Parameters at location
-   * @param numParameters Number of parameters
-   * @param totalStrain Total strain at locations
-   * @param numLocs Number of locations
+   * @param stress Array for stress tensor
+   * @param size Size of array for stress tensor
+   * @param parameters Parameters at locations.
+   * @param numParameters Number of parameters.
+   * @param totalStrain Total strain at locations.
    * @param spaceDim Spatial dimension for locations.
    */
-  void _calcStress(const double* parameters,
+  void _calcStress(double* const stress,
+		   const int size,
+		   const double* parameters,
 		   const int numParameters,
 		   const double* totalStrain,
-		   const int numLocs,
 		   const int spaceDim);
 
-  /** Compute derivatives of elasticity matrix at locations from parameters.
+  /** Compute derivatives of elasticity matrix from parameters.
    *
-   * Results are stored in _elasticConsts.
-   *
+   * @param elasticConsts Array for elastic constants
+   * @param size Size of array
    * @param parameters Parameters at locations.
    * @param numParameters Number of parameters.
    * @param totalStrain Total strain at locations.
-   * @param numLocs Number of locations.
    * @param spaceDim Spatial dimension for locations.
    */
-  void _calcElasticConsts(const double* parameters,
+  void _calcElasticConsts(double* const elasticConsts,
+			  const int size,
+			  const double* parameters,
 			  const int numParameters,
 			  const double* totalStrain,
-			  const int numLocs,
 			  const int spaceDim);
 
   // NOT IMPLEMENTED ////////////////////////////////////////////////////

Added: short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc	2007-04-08 02:00:47 UTC (rev 6521)
@@ -0,0 +1,219 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "ElasticStrain1D.hh" // implementation of object methods
+
+#include <assert.h> // USES assert()
+
+// ----------------------------------------------------------------------
+namespace pylith {
+  namespace materials {
+    class _ElasticStrain1D;
+  } // materials
+} // pylith
+
+/** _ElasticStrain1D is a helper class for ElasticStrain1D. We define
+ * it in this implementation file to insulate other objects from these
+ * details.
+ */
+class pylith::materials::_ElasticStrain1D {
+public:
+  /// Number of entries in stress tensor.
+  static const int stressSize;
+
+  /// Number of entries in derivative of elasticity matrix.
+  static const int numElasticConsts;
+
+  /// Values expected in spatial database
+  static const int numDBValues;
+  static const char* namesDBValues[];
+
+  /// Indices (order) of database values
+  static const int didDensity;
+  static const int didVp;
+
+  /// Parameters
+  static const int numParameters;
+  static const char* namesParameters[];
+
+  /// Indices (order) of parameters
+  static const int pidDensity;
+  static const int pidLambda2mu;
+}; // _ElasticStrain1D
+
+const int pylith::materials::_ElasticStrain1D::stressSize = 1;
+const int pylith::materials::_ElasticStrain1D::numElasticConsts = 1;
+const int pylith::materials::_ElasticStrain1D::numDBValues = 2;
+const char* pylith::materials::_ElasticStrain1D::namesDBValues[] =
+  {"density", "vp" };
+const int pylith::materials::_ElasticStrain1D::numParameters = 2;
+const char* pylith::materials::_ElasticStrain1D::namesParameters[] =
+  {"density", "lambda2mu" };
+const int pylith::materials::_ElasticStrain1D::didDensity = 0;
+const int pylith::materials::_ElasticStrain1D::didVp = 1;
+const int pylith::materials::_ElasticStrain1D::pidDensity = 0;
+const int pylith::materials::_ElasticStrain1D::pidLambda2mu = 1;
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::materials::ElasticStrain1D::ElasticStrain1D(void)
+{ // constructor
+  _dimension = 1;
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::materials::ElasticStrain1D::~ElasticStrain1D(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Copy constructor.
+pylith::materials::ElasticStrain1D::ElasticStrain1D(
+						const ElasticStrain1D& m) :
+  ElasticMaterial(m)
+{ // copy constructor
+} // copy constructor
+
+// ----------------------------------------------------------------------
+// Get number of entries in stress tensor.
+int
+pylith::materials::ElasticStrain1D::stressSize(void) const
+{ // stressSize
+  return _ElasticStrain1D::stressSize;
+} // stressSize
+
+// ----------------------------------------------------------------------
+// Get number of elastic constants for material.
+int
+pylith::materials::ElasticStrain1D::numElasticConsts(void) const
+{ // numElasticConsts
+  return _ElasticStrain1D::numElasticConsts;
+} // numElasticConsts
+
+// ----------------------------------------------------------------------
+// Get names of values expected to be in database of parameters for
+const char**
+pylith::materials::ElasticStrain1D::_dbValues(void) const
+{ // _dbValues
+  return _ElasticStrain1D::namesDBValues;
+} // _dbValues
+
+// ----------------------------------------------------------------------
+// Get number of values expected to be in database of parameters for
+int
+pylith::materials::ElasticStrain1D::_numDBValues(void) const
+{ // _numDBValues
+  return _ElasticStrain1D::numDBValues;
+} // _numDBValues
+
+// ----------------------------------------------------------------------
+// Get names of parameters for physical properties.
+const char**
+pylith::materials::ElasticStrain1D::_parameterNames(void) const
+{ // _parameterNames
+  return _ElasticStrain1D::namesParameters;
+} // _parameterNames
+
+// ----------------------------------------------------------------------
+// Get number of parameters for physical properties.
+int
+pylith::materials::ElasticStrain1D::_numParameters(void) const
+{ // _numParameters
+  return _ElasticStrain1D::numParameters;
+} // _numParameters
+
+// ----------------------------------------------------------------------
+// Compute parameters from values in spatial database.
+void
+pylith::materials::ElasticStrain1D::_dbToParameters(double* paramVals,
+						    const int numParams,
+						    const double* dbValues,
+						    const int numValues) const
+{ // computeParameters
+  assert(0 != paramVals);
+  assert(_ElasticStrain1D::numParameters == numParams);
+  assert(0 != dbValues);
+  assert(_ElasticStrain1D::numDBValues == numValues);
+
+  const double density = dbValues[_ElasticStrain1D::didDensity];
+  const double vp = dbValues[_ElasticStrain1D::didVp];
+  const double lambda2mu = density * vp*vp;
+
+  paramVals[_ElasticStrain1D::pidDensity] = density;
+  paramVals[_ElasticStrain1D::pidLambda2mu] = lambda2mu;
+} // computeParameters
+
+// ----------------------------------------------------------------------
+// Compute density at location from parameters.
+void
+pylith::materials::ElasticStrain1D::_calcDensity(double* const density,
+						 const int size,
+						 const double* parameters,
+						 const int numParameters)
+{ // _calcDensity
+  assert(0 != density);
+  assert(1 == size);
+  assert(0 != parameters);
+
+  *density = parameters[_ElasticStrain1D::pidDensity];
+} // _calcDensity
+
+// ----------------------------------------------------------------------
+// Compute stress tensor at location from parameters.
+void
+pylith::materials::ElasticStrain1D::_calcStress(double* const stress,
+						const int size,
+						const double* parameters,
+						const int numParameters,
+						const double* totalStrain,
+						const int spaceDim)
+{ // _calcStress
+  assert(0 != stress);
+  assert(_ElasticStrain1D::stressSize == size);
+  assert(0 != parameters);
+  assert(_ElasticStrain1D::numParameters == numParameters);
+  assert(spaceDim == _dimension);
+
+  const double density = parameters[_ElasticStrain1D::pidDensity];
+  const double lambda2mu = parameters[_ElasticStrain1D::pidLambda2mu];
+
+  const double e11 = totalStrain[0];
+  stress[0] = lambda2mu * e11;
+} // _calcStress
+
+// ----------------------------------------------------------------------
+// Compute derivative of elasticity matrix at location from parameters.
+void
+pylith::materials::ElasticStrain1D::_calcElasticConsts(double* const elasticConsts,
+						       const int size,
+						       const double* parameters,
+						       const int numParameters,
+						       const double* totalStrain,
+						       const int spaceDim)
+{ // _calcElasticConsts
+  assert(0 != elasticConsts);
+  assert(_ElasticStrain1D::numElasticConsts == size);
+  assert(0 != parameters);
+  assert(_ElasticStrain1D::numParameters == numParameters);
+  assert(spaceDim == _dimension);
+ 
+  const double density = parameters[_ElasticStrain1D::pidDensity];
+  const double lambda2mu = parameters[_ElasticStrain1D::pidLambda2mu];
+
+  elasticConsts[0] = lambda2mu;
+} // _calcElasticConsts
+
+
+// End of file 

Added: short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh	2007-04-08 02:00:47 UTC (rev 6521)
@@ -0,0 +1,184 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/materials/ElasticStrain1D.h
+ *
+ * @brief C++ ElasticStrain1D object
+ *
+ * 1-D, linear elastic material with axial strain. The physical
+ * properties are specified using density and compressional-wave
+ * speed. The physical properties are stored internally using density
+ * and lambda + 2 mu, which are directly related to the elasticity
+ * constants used in the finite-element integration.
+ */
+
+#if !defined(pylith_materials_elasticstrain1d_hh)
+#define pylith_materials_elasticstrain1d_hh
+
+#include "ElasticMaterial.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+  namespace materials {
+    class ElasticStrain1D;
+    class TestElasticStrain1D; // unit testing
+  } // materials
+} // pylith
+
+/// 3-D, linear elastic material with axial strain.
+class pylith::materials::ElasticStrain1D : public ElasticMaterial
+{ // class ElasticStrain1D
+  friend class TestElasticStrain1D; // unit testing
+
+  // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+  /// Default constructor
+  ElasticStrain1D(void);
+
+  /// Destructor
+  ~ElasticStrain1D(void);
+
+  /** Create a pointer to a copy of this.
+   *
+   * @returns Pointer to copy
+   */
+  ElasticMaterial* clone(void) const;
+
+  /** Get number of entries in stress tensor.
+   *
+   * 1-D = 1
+   * 2-D = 3
+   * 3-D = 6
+   *
+   * @returns Number of entries in stress tensor.
+   */
+  int stressSize(void) const;
+
+  /** Get number of elastic constants for material.
+   *
+   * 1-D = 1
+   * 2-D = 6
+   * 3-D = 21
+   *
+   * @returns Number of elastic constants
+   */
+  int numElasticConsts(void) const;
+
+  // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+  /** Copy constructor
+   *
+   * @param m Material to copy
+   */
+  ElasticStrain1D(const ElasticStrain1D& m);
+
+  /** Get names of values expected to be in database of parameters for
+   *  physical properties.
+   *
+   * @returns Names of values
+   */
+  const char** _dbValues(void) const;
+
+  /** Get number of values expected to be in database of parameters for
+   *  physical properties.
+   *
+   * @returns Number of values
+   */
+  int _numDBValues(void) const;
+
+  /** Get names of parameters for physical properties.
+   *
+   * @returns Names of parameters
+   */
+  const char** _parameterNames(void) const;
+
+  /** Get number of parameters for physical properties.
+   *
+   * @returns Number of parameters
+   */
+  int _numParameters(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
+   * @param dbValues Array of database values
+   * @param numValues Number of database values
+   */
+  void _dbToParameters(double* paramVals,
+		       const int numParams,
+		       const double* dbValues,
+		       const int numValues) const;
+
+  /** Compute density from parameters.
+   *
+   * @param density Array for density
+   * @param size Size of array for density
+   * @param parameters Parameters at location
+   * @param numParameters Number of parameters
+   */
+  void _calcDensity(double* const density,
+		    const int size,
+		    const double* parameters,
+		    const int numParameters);
+
+  /** Compute stress tensor from parameters.
+   *
+   * @param stress Array for stress tensor
+   * @param size Size of array for stress tensor
+   * @param parameters Parameters at locations.
+   * @param numParameters Number of parameters.
+   * @param totalStrain Total strain at locations.
+   * @param spaceDim Spatial dimension for locations.
+   */
+  void _calcStress(double* const stress,
+		   const int size,
+		   const double* parameters,
+		   const int numParameters,
+		   const double* totalStrain,
+		   const int spaceDim);
+
+  /** Compute derivatives of elasticity matrix from parameters.
+   *
+   * @param elasticConsts Array for elastic constants
+   * @param size Size of array
+   * @param parameters Parameters at locations.
+   * @param numParameters Number of parameters.
+   * @param totalStrain Total strain at locations.
+   * @param spaceDim Spatial dimension for locations.
+   */
+  void _calcElasticConsts(double* const elasticConsts,
+			  const int size,
+			  const double* parameters,
+			  const int numParameters,
+			  const double* totalStrain,
+			  const int spaceDim);
+
+  // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+  /// Not implemented
+  const ElasticStrain1D& operator=(const ElasticStrain1D& m);
+
+}; // class ElasticStrain1D
+
+#include "ElasticStrain1D.icc" // inline methods
+
+#endif // pylith_materials_elasticstrain1d_hh
+
+
+// End of file 

Added: short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.icc	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.icc	2007-04-08 02:00:47 UTC (rev 6521)
@@ -0,0 +1,25 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_materials_elasticstrain1d_hh)
+#error "ElasticStrain1D.icc can only be included from ElasticStrain1D.hh"
+#endif
+
+// Create a pointer to a copy of this.
+inline
+pylith::materials::ElasticMaterial*
+pylith::materials::ElasticStrain1D::clone(void) const {
+  return new ElasticStrain1D(*this);
+}
+
+
+// End of file 

Copied: short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc (from rev 6518, short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.cc)
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.cc	2007-04-07 19:05:26 UTC (rev 6518)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc	2007-04-08 02:00:47 UTC (rev 6521)
@@ -0,0 +1,229 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "ElasticStress1D.hh" // implementation of object methods
+
+#include <assert.h> // USES assert()
+
+// ----------------------------------------------------------------------
+namespace pylith {
+  namespace materials {
+    class _ElasticStress1D;
+  } // materials
+} // pylith
+
+/** _ElasticStress1D is a helper class for ElasticStress1D. We define
+ * it in this implementation file to insulate other objects from these
+ * details.
+ */
+class pylith::materials::_ElasticStress1D {
+public:
+  /// Number of entries in stress tensor.
+  static const int stressSize;
+
+  /// Number of entries in derivative of elasticity matrix.
+  static const int numElasticConsts;
+
+  /// Values expected in spatial database
+  static const int numDBValues;
+  static const char* namesDBValues[];
+
+  /// Indices (order) of database values
+  static const int didDensity;
+  static const int didVs;
+  static const int didVp;
+
+  /// Parameters
+  static const int numParameters;
+  static const char* namesParameters[];
+
+  /// Indices (order) of parameters
+  static const int pidDensity;
+  static const int pidMu;
+  static const int pidLambda;
+}; // _ElasticStress1D
+
+const int pylith::materials::_ElasticStress1D::stressSize = 1;
+const int pylith::materials::_ElasticStress1D::numElasticConsts = 1;
+const int pylith::materials::_ElasticStress1D::numDBValues = 3;
+const char* pylith::materials::_ElasticStress1D::namesDBValues[] =
+  {"density", "vs", "vp" };
+const int pylith::materials::_ElasticStress1D::numParameters = 3;
+const char* pylith::materials::_ElasticStress1D::namesParameters[] =
+  {"density", "mu", "lambda" };
+const int pylith::materials::_ElasticStress1D::didDensity = 0;
+const int pylith::materials::_ElasticStress1D::didVs = 1;
+const int pylith::materials::_ElasticStress1D::didVp = 2;
+const int pylith::materials::_ElasticStress1D::pidDensity = 0;
+const int pylith::materials::_ElasticStress1D::pidMu = 1;
+const int pylith::materials::_ElasticStress1D::pidLambda = 2;
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::materials::ElasticStress1D::ElasticStress1D(void)
+{ // constructor
+  _dimension = 1;
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::materials::ElasticStress1D::~ElasticStress1D(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Copy constructor.
+pylith::materials::ElasticStress1D::ElasticStress1D(
+						const ElasticStress1D& m) :
+  ElasticMaterial(m)
+{ // copy constructor
+} // copy constructor
+
+// ----------------------------------------------------------------------
+// Get number of entries in stress tensor.
+int
+pylith::materials::ElasticStress1D::stressSize(void) const
+{ // stressSize
+  return _ElasticStress1D::stressSize;
+} // stressSize
+
+// ----------------------------------------------------------------------
+// Get number of elastic constants for material.
+int
+pylith::materials::ElasticStress1D::numElasticConsts(void) const
+{ // numElasticConsts
+  return _ElasticStress1D::numElasticConsts;
+} // numElasticConsts
+
+// ----------------------------------------------------------------------
+// Get names of values expected to be in database of parameters for
+const char**
+pylith::materials::ElasticStress1D::_dbValues(void) const
+{ // _dbValues
+  return _ElasticStress1D::namesDBValues;
+} // _dbValues
+
+// ----------------------------------------------------------------------
+// Get number of values expected to be in database of parameters for
+int
+pylith::materials::ElasticStress1D::_numDBValues(void) const
+{ // _numDBValues
+  return _ElasticStress1D::numDBValues;
+} // _numDBValues
+
+// ----------------------------------------------------------------------
+// Get names of parameters for physical properties.
+const char**
+pylith::materials::ElasticStress1D::_parameterNames(void) const
+{ // _parameterNames
+  return _ElasticStress1D::namesParameters;
+} // _parameterNames
+
+// ----------------------------------------------------------------------
+// Get number of parameters for physical properties.
+int
+pylith::materials::ElasticStress1D::_numParameters(void) const
+{ // _numParameters
+  return _ElasticStress1D::numParameters;
+} // _numParameters
+
+// ----------------------------------------------------------------------
+// Compute parameters from values in spatial database.
+void
+pylith::materials::ElasticStress1D::_dbToParameters(double* paramVals,
+						       const int numParams,
+						       const double* dbValues,
+						       const int numValues) const
+{ // computeParameters
+  assert(0 != paramVals);
+  assert(_ElasticStress1D::numParameters == numParams);
+  assert(0 != dbValues);
+  assert(_ElasticStress1D::numDBValues == numValues);
+
+  const double density = dbValues[_ElasticStress1D::didDensity];
+  const double vs = dbValues[_ElasticStress1D::didVs];
+  const double vp = dbValues[_ElasticStress1D::didVp];
+ 
+  const double mu = density * vs*vs;
+  const double lambda = density * vp*vp - 2.0*mu;
+
+  paramVals[_ElasticStress1D::pidDensity] = density;
+  paramVals[_ElasticStress1D::pidMu] = mu;
+  paramVals[_ElasticStress1D::pidLambda] = lambda;
+} // computeParameters
+
+// ----------------------------------------------------------------------
+// Compute density at location from parameters.
+void
+pylith::materials::ElasticStress1D::_calcDensity(double* const density,
+						    const int size,
+						    const double* parameters,
+						    const int numParameters)
+{ // _calcDensity
+  assert(0 != density);
+  assert(1 == size);
+  assert(0 != parameters);
+
+  *density = parameters[_ElasticStress1D::pidDensity];
+} // _calcDensity
+
+// ----------------------------------------------------------------------
+// Compute stress tensor at location from parameters.
+void
+pylith::materials::ElasticStress1D::_calcStress(double* const stress,
+						   const int size,
+						   const double* parameters,
+						   const int numParameters,
+						   const double* totalStrain,
+						   const int spaceDim)
+{ // _calcStress
+  assert(0 != stress);
+  assert(_ElasticStress1D::stressSize == size);
+  assert(0 != parameters);
+  assert(_ElasticStress1D::numParameters == numParameters);
+  assert(spaceDim == _dimension);
+
+  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;
+} // _calcStress
+
+// ----------------------------------------------------------------------
+// Compute derivative of elasticity matrix at location from parameters.
+void
+pylith::materials::ElasticStress1D::_calcElasticConsts(double* const elasticConsts,
+						       const int size,
+						       const double* parameters,
+						       const int numParameters,
+						       const double* totalStrain,
+						       const int spaceDim)
+{ // _calcElasticConsts
+  assert(0 != elasticConsts);
+  assert(_ElasticStress1D::numElasticConsts == size);
+  assert(0 != parameters);
+  assert(_ElasticStress1D::numParameters == numParameters);
+  assert(spaceDim == _dimension);
+ 
+  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);
+} // _calcElasticConsts
+
+
+// End of file 

Copied: short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh (from rev 6518, short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.hh)
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.hh	2007-04-07 19:05:26 UTC (rev 6518)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh	2007-04-08 02:00:47 UTC (rev 6521)
@@ -0,0 +1,185 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/materials/ElasticStress1D.h
+ *
+ * @brief C++ ElasticStress1D object
+ *
+ * 1-D, linear elastic material with axial stress. The physical
+ * properties are specified using density, shear-wave speed, and
+ * compressional-wave speed. The physical properties are stored
+ * internally using density, mu, and lambda, which are directly
+ * related to the elasticity constants used in the finite-element
+ * integration.
+ */
+
+#if !defined(pylith_materials_elasticstress1d_hh)
+#define pylith_materials_elasticstress1d_hh
+
+#include "ElasticMaterial.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+  namespace materials {
+    class ElasticStress1D;
+    class TestElasticStress1D; // unit testing
+  } // materials
+} // pylith
+
+/// 1-D, linear elastic material with axial stres.
+class pylith::materials::ElasticStress1D : public ElasticMaterial
+{ // class ElasticStress1D
+  friend class TestElasticStress1D; // unit testing
+
+  // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+  /// Default constructor
+  ElasticStress1D(void);
+
+  /// Destructor
+  ~ElasticStress1D(void);
+
+  /** Create a pointer to a copy of this.
+   *
+   * @returns Pointer to copy
+   */
+  ElasticMaterial* clone(void) const;
+
+  /** Get number of entries in stress tensor.
+   *
+   * 1-D = 1
+   * 2-D = 3
+   * 3-D = 6
+   *
+   * @returns Number of entries in stress tensor.
+   */
+  int stressSize(void) const;
+
+  /** Get number of elastic constants for material.
+   *
+   * 1-D = 1
+   * 2-D = 6
+   * 3-D = 21
+   *
+   * @returns Number of elastic constants
+   */
+  int numElasticConsts(void) const;
+
+  // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+  /** Copy constructor
+   *
+   * @param m Material to copy
+   */
+  ElasticStress1D(const ElasticStress1D& m);
+
+  /** Get names of values expected to be in database of parameters for
+   *  physical properties.
+   *
+   * @returns Names of values
+   */
+  const char** _dbValues(void) const;
+
+  /** Get number of values expected to be in database of parameters for
+   *  physical properties.
+   *
+   * @returns Number of values
+   */
+  int _numDBValues(void) const;
+
+  /** Get names of parameters for physical properties.
+   *
+   * @returns Names of parameters
+   */
+  const char** _parameterNames(void) const;
+
+  /** Get number of parameters for physical properties.
+   *
+   * @returns Number of parameters
+   */
+  int _numParameters(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
+   * @param dbValues Array of database values
+   * @param numValues Number of database values
+   */
+  void _dbToParameters(double* paramVals,
+		       const int numParams,
+		       const double* dbValues,
+		       const int numValues) const;
+
+  /** Compute density from parameters.
+   *
+   * @param density Array for density
+   * @param size Size of array for density
+   * @param parameters Parameters at location
+   * @param numParameters Number of parameters
+   */
+  void _calcDensity(double* const density,
+		    const int size,
+		    const double* parameters,
+		    const int numParameters);
+
+  /** Compute stress tensor from parameters.
+   *
+   * @param stress Array for stress tensor
+   * @param size Size of array for stress tensor
+   * @param parameters Parameters at locations.
+   * @param numParameters Number of parameters.
+   * @param totalStrain Total strain at locations.
+   * @param spaceDim Spatial dimension for locations.
+   */
+  void _calcStress(double* const stress,
+		   const int size,
+		   const double* parameters,
+		   const int numParameters,
+		   const double* totalStrain,
+		   const int spaceDim);
+
+  /** Compute derivatives of elasticity matrix from parameters.
+   *
+   * @param elasticConsts Array for elastic constants
+   * @param size Size of array
+   * @param parameters Parameters at locations.
+   * @param numParameters Number of parameters.
+   * @param totalStrain Total strain at locations.
+   * @param spaceDim Spatial dimension for locations.
+   */
+  void _calcElasticConsts(double* const elasticConsts,
+			  const int size,
+			  const double* parameters,
+			  const int numParameters,
+			  const double* totalStrain,
+			  const int spaceDim);
+
+  // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+  /// Not implemented
+  const ElasticStress1D& operator=(const ElasticStress1D& m);
+
+}; // class ElasticStress1D
+
+#include "ElasticStress1D.icc" // inline methods
+
+#endif // pylith_materials_elasticstress1d_hh
+
+
+// End of file 

Copied: short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.icc (from rev 6518, short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.icc)
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic1D.icc	2007-04-07 19:05:26 UTC (rev 6518)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.icc	2007-04-08 02:00:47 UTC (rev 6521)
@@ -0,0 +1,25 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_materials_elasticstress1d_hh)
+#error "ElasticStress1D.icc can only be included from ElasticStress1D.hh"
+#endif
+
+// Create a pointer to a copy of this.
+inline
+pylith::materials::ElasticMaterial*
+pylith::materials::ElasticStress1D::clone(void) const {
+  return new ElasticStress1D(*this);
+}
+
+
+// End of file 

Modified: short/3D/PyLith/trunk/libsrc/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Makefile.am	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/libsrc/materials/Makefile.am	2007-04-08 02:00:47 UTC (rev 6521)
@@ -14,8 +14,10 @@
 include $(top_srcdir)/subpackage.am
 
 subpkginclude_HEADERS = \
-	ElasticIsotropic1D.hh \
-	ElasticIsotropic1D.icc \
+	ElasticStrain1D.hh \
+	ElasticStrain1D.icc \
+	ElasticStress1D.hh \
+	ElasticStress1D.icc \
 	ElasticIsotropic3D.hh \
 	ElasticIsotropic3D.icc \
 	ElasticMaterial.hh \

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc	2007-04-07 22:16:54 UTC (rev 6520)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc	2007-04-08 02:00:47 UTC (rev 6521)
@@ -335,19 +335,14 @@
 					ElasticMaterial* material,
 					const ElasticMaterialData& data) const
 { // _testCalcDensity
-  const int numQuadPts = data.numLocs;
-  material->initCellData(numQuadPts);
-  material->_calcDensity(data.parameterData, data.numParameters, data.numLocs);
+  double density = 0;
+  material->_calcDensity(&density, 1, data.parameterData, data.numParameters);
   const double* densityE = data.density;
-  const double* density = material->_density;
   
-  CPPUNIT_ASSERT(0 != density);
   CPPUNIT_ASSERT(0 != densityE);
 
   const double tolerance = 1.0e-06;
-  const double size = numQuadPts;
-  for (int i=0; i < size; ++i)
-    CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, density[i]/densityE[i], tolerance);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, density/densityE[0], tolerance);
 } // _testCalcDensity
 
 // ----------------------------------------------------------------------
@@ -357,18 +352,17 @@
 				       ElasticMaterial* material,
 				       const ElasticMaterialData& data) const
 { // _testCalcElasticConsts
-  const int numQuadPts = data.numLocs;
-  material->initCellData(numQuadPts);
-  material->_calcStress(data.parameterData, data.numParameters, data.strain, 
-			data.numLocs, data.dimension);
+  const int size = material->stressSize();
+  double* stress = (size > 0) ? new double[size] : 0;
+  material->_calcStress(stress, size, 
+			data.parameterData, data.numParameters, 
+			data.strain, data.dimension);
   const double* stressE = data.stress;
-  const double* stress = material->_stress;
   
   CPPUNIT_ASSERT(0 != stress);
   CPPUNIT_ASSERT(0 != stressE);
 
   const double tolerance = 1.0e-06;
-  const double size = numQuadPts * material->stressSize();
   for (int i=0; i < size; ++i)
     if (fabs(stressE[i]) > tolerance)
       CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, stress[i]/stressE[i], 
@@ -376,6 +370,7 @@
     else
       CPPUNIT_ASSERT_DOUBLES_EQUAL(stressE[i], stress[i],
 				   tolerance);
+  delete[] stress; stress = 0;
 } // _testCalcStress
 
 // ----------------------------------------------------------------------
@@ -385,18 +380,17 @@
 				       ElasticMaterial* material,
 				       const ElasticMaterialData& data) const
 { // _testCalcElasticConsts
-  const int numQuadPts = data.numLocs;
-  material->initCellData(numQuadPts);
-  material->_calcElasticConsts(data.parameterData, data.numParameters, 
-			       data.strain, data.numLocs, data.dimension);
+  const int size = material->numElasticConsts();
+  double* elasticConsts = (size > 0) ? new double[size] : 0;
+  material->_calcElasticConsts(elasticConsts, size,
+			       data.parameterData, data.numParameters, 
+			       data.strain, data.dimension);
   const double* elasticConstsE = data.elasticConsts;
-  const double* elasticConsts = material->_elasticConsts;
   
   CPPUNIT_ASSERT(0 != elasticConsts);
   CPPUNIT_ASSERT(0 != elasticConstsE);
 
   const double tolerance = 1.0e-06;
-  const double size = numQuadPts * material->numElasticConsts();
   for (int i=0; i < size; ++i)
     if (fabs(elasticConstsE[i]) > tolerance)
       CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, elasticConsts[i]/elasticConstsE[i], 
@@ -404,6 +398,7 @@
     else
       CPPUNIT_ASSERT_DOUBLES_EQUAL(elasticConstsE[i], elasticConsts[i],
 				   tolerance);
+  delete[] elasticConsts; elasticConsts = 0;
 } // _testCalcElasticConsts