[cig-commits] r15282 - in short/3D/PyLith/trunk: . modulesrc modulesrc/materials templates templates/materials
brad at geodynamics.org
brad at geodynamics.org
Mon Jun 15 15:54:23 PDT 2009
Author: brad
Date: 2009-06-15 15:54:22 -0700 (Mon, 15 Jun 2009)
New Revision: 15282
Added:
short/3D/PyLith/trunk/templates/
short/3D/PyLith/trunk/templates/README
short/3D/PyLith/trunk/templates/materials/
short/3D/PyLith/trunk/templates/materials/Makefile.am
short/3D/PyLith/trunk/templates/materials/PlaneStrainState.cc
short/3D/PyLith/trunk/templates/materials/PlaneStrainState.hh
short/3D/PyLith/trunk/templates/materials/PlaneStrainState.i
short/3D/PyLith/trunk/templates/materials/PlaneStrainState.py
short/3D/PyLith/trunk/templates/materials/README
short/3D/PyLith/trunk/templates/materials/__init__.py
short/3D/PyLith/trunk/templates/materials/configure.ac
short/3D/PyLith/trunk/templates/materials/materialscontrib.i
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/configure.ac
short/3D/PyLith/trunk/modulesrc/Makefile.am
short/3D/PyLith/trunk/modulesrc/materials/Makefile.am
Log:
Setting up constitutive model template.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-06-15 22:25:53 UTC (rev 15281)
+++ short/3D/PyLith/trunk/TODO 2009-06-15 22:54:22 UTC (rev 15282)
@@ -16,8 +16,6 @@
Extensions
- Physical properties
-
Constitutive model
Time dependent Neumann BC
@@ -90,8 +88,6 @@
Extending PyLith [BRAD]
- Constitutive models
-
Material properties
* Reduce memory use with ordering elements by material??
Modified: short/3D/PyLith/trunk/configure.ac
===================================================================
--- short/3D/PyLith/trunk/configure.ac 2009-06-15 22:25:53 UTC (rev 15281)
+++ short/3D/PyLith/trunk/configure.ac 2009-06-15 22:54:22 UTC (rev 15282)
@@ -221,6 +221,7 @@
libsrc/topology/Makefile
libsrc/utils/Makefile
modulesrc/Makefile
+ modulesrc/include/Makefile
modulesrc/bc/Makefile
modulesrc/faults/Makefile
modulesrc/feassemble/Makefile
Modified: short/3D/PyLith/trunk/modulesrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/Makefile.am 2009-06-15 22:25:53 UTC (rev 15281)
+++ short/3D/PyLith/trunk/modulesrc/Makefile.am 2009-06-15 22:54:22 UTC (rev 15282)
@@ -11,6 +11,7 @@
#
SUBDIRS = \
+ include \
bc \
faults \
feassemble \
Modified: short/3D/PyLith/trunk/modulesrc/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/Makefile.am 2009-06-15 22:25:53 UTC (rev 15281)
+++ short/3D/PyLith/trunk/modulesrc/materials/Makefile.am 2009-06-15 22:54:22 UTC (rev 15282)
@@ -30,6 +30,11 @@
GenMaxwellIsotropic3D.i \
PowerLaw3D.i
+swigincludedir = $(pkgdatadir)/swig/$(subpackage)
+swiginclude_HEADERS = \
+ Material.i \
+ ElasticMaterial.i
+
swig_generated = \
materials_wrap.cxx \
materials.py
Added: short/3D/PyLith/trunk/templates/README
===================================================================
--- short/3D/PyLith/trunk/templates/README (rev 0)
+++ short/3D/PyLith/trunk/templates/README 2009-06-15 22:54:22 UTC (rev 15282)
@@ -0,0 +1,3 @@
+This directory contains templates for extending the functionality of
+PyLith. This includes creating user-defined bulk constitutive models
+(see the materials directory).
Property changes on: short/3D/PyLith/trunk/templates/materials
___________________________________________________________________
Name: svn:externals
+ m4 http://geodynamics.org/svn/cig/cs/autoconf/trunk
Added: short/3D/PyLith/trunk/templates/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/templates/materials/Makefile.am (rev 0)
+++ short/3D/PyLith/trunk/templates/materials/Makefile.am 2009-06-15 22:54:22 UTC (rev 15282)
@@ -0,0 +1,88 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+ACLOCAL_AMFLAGS = -I ./m4
+
+SUBDIRS = tests
+
+subpkgpythondir = $(pythondir)/pylith/materials/contrib
+subpkgpyexecdir = $(pyexecdir)/pylith/materials/contrib
+
+
+# LIBRARY --------------------------------------------------------------
+lib_LTLIBRARIES = libmaterialscontrib.la
+
+libmaterialscontrib_la_SOURCES = \
+ PlaneStrainState.cc
+
+noinst_HEADERS = \
+ PlaneStrainState.hh
+
+libmaterialscontrib_la_LDFLAGS = $(AM_LDFLAGS)
+
+libmaterialscontrib_la_LIBADD = \
+ -lpylith \
+ -lproj
+if NO_UNDEFINED
+libmaterialscontrib_la_LIBADD += \
+ $(PYTHON_BLDLIBRARY) $(PYTHON_LIBS) $(PYTHON_SYSLIBS)
+endif
+
+AM_CPPFLAGS = $(PYTHON_EGG_CPPFLAGS) -I$(PYTHON_INCDIR)
+
+INCLUDES =
+
+# MODULE ---------------------------------------------------------------
+
+subpkgpyexec_LTLIBRARIES = _materialscontribmodule.la
+
+subpkgpyexec_PYTHON = \
+ materialscontrib.py \
+ __init__.py
+
+swig_sources = \
+ materialscontrib.i \
+ PlaneStrainState.i
+
+swig_generated = \
+ materialscontrib_wrap.cxx \
+ materialscontrib.py
+
+_materialscontribmodule_la_LDFLAGS = -module -avoid-version \
+ $(AM_LDFLAGS) $(PYTHON_LA_LDFLAGS)
+
+dist__materialscontribmodule_la_SOURCES = $(swig_sources) $(swig_generated)
+
+_materialscontribmodule_la_LIBADD = libmaterialscontrib.la
+if NO_UNDEFINED
+_materialscontribmodule_la_LIBADD += \
+ $(PYTHON_BLDLIBRARY) $(PYTHON_LIBS) $(PYTHON_SYSLIBS)
+endif
+
+INCLUDES += -I$(NUMPY_INCDIR) -I$(PYTHON_INCDIR)
+
+$(srcdir)/materialscontrib_wrap.cxx $(srcdir)/materialscontrib.py: $(swig_sources)
+ $(SWIG) -Wall -c++ -python $(SWIG_FLAGS) $<
+
+
+MAINTAINERCLEANFILES = $(swig_generated)
+
+
+
+# PYTHON ---------------------------------------------------------------
+
+nobase_subpkgpyexec_PYTHON = \
+ __init__.py \
+ PlaneStrainState.py
+
+
+# End of file
Added: short/3D/PyLith/trunk/templates/materials/PlaneStrainState.cc
===================================================================
--- short/3D/PyLith/trunk/templates/materials/PlaneStrainState.cc (rev 0)
+++ short/3D/PyLith/trunk/templates/materials/PlaneStrainState.cc 2009-06-15 22:54:22 UTC (rev 15282)
@@ -0,0 +1,449 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "PlaneStrainState.hh" // implementation of object methods
+
+#include "pylith/materials/Metadata.hh" // USES Metadata
+
+#include "pylith/utils/array.hh" // USES double_array
+#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
+//#include "petsc.h" // USES PetscLogFlops
+
+#include <cassert> // USES assert()
+#include <sstream> // USES std::ostringstream
+#include <stdexcept> // USES std::runtime_error
+
+// ----------------------------------------------------------------------
+// Create a local namespace to use for local constants and other
+// information. This insulates all other classes from this information
+// while preventing clashes with other local constants and data (as
+// long as no other object use the same _PlaneStrainState namespace in
+// the contrib::materials namespace.
+namespace contrib {
+ namespace materials {
+ namespace _PlaneStrainState {
+
+ // Dimension of material.
+ const int dimension = 2;
+
+ // Number of entries in stress tensor.
+ const int tensorSize = 3;
+
+ // Number of elastic constants (for general 3-D elastic material)
+ const int numElasticConsts = 6;
+
+ // These are the physical properties stored during the
+ // simulation and need not coincide with the physical properties
+ // provided by the user. We store Lame's constants and density,
+ // but the user provides Vs, Vp, and density.
+
+ // Number of physical properties.
+ const int numProperties = 3;
+
+ // Physical properties.
+ const Metadata::ParamDescription properties[] = {
+ { "density", 1, pylith::topology::FieldBase::SCALAR },
+ { "mu", 1, pylith::topology::FieldBase::SCALAR },
+ { "lambda", 1, pylith::topology::FieldBase::SCALAR },
+ };
+
+ // Values expected in spatial database
+ const int numDBProperties = 3;
+ const char* dbProperties[] = { "density", "vs", "vp" };
+
+ // These are the state variables stored during the
+ // simulation. Usually, we store only the time-dependent values
+ // needed to compute the behavior at a given point in time. In
+ // this example, however, for illustration purposes we store the
+ // elastic strain (total strain) tensor and the stress tensor.
+
+ /// Number of state variables. we store
+ const int numStateVars = 2;
+
+ /// State variables.
+ const Metadata::ParamDescription stateVars[] = {
+ { "total_strain", tensorSize, pylith::topology::FieldBase::TENSOR },
+ { "stress", tensorSize, pylith::topology::FieldBase::TENSOR },
+ };
+
+ // We do not include a list of values expected to be in the
+ // initial state variables spatial database because the initial
+ // state variables are the initial stress and strain, which are
+ // already provided as special spatial databases.
+
+ } // _PlaneStrainState
+ } // materials
+} // contrib
+
+// Indices of physical properties in the properties array.
+const int contrib::materials::PlaneStrainState::p_density = 0;
+
+const int contrib::materials::PlaneStrainState::p_mu =
+ contrib::materials::PlaneStrainState::p_density + 1;
+
+const int contrib::materials::PlaneStrainState::p_lambda =
+ contrib::materials::PlaneStrainState::p_mu + 1;
+
+// Indices of database values (order must match dbProperties)
+const int contrib::materials::PlaneStrainState::db_density = 0;
+
+const int contrib::materials::PlaneStrainState::db_vs =
+ contrib::materials::PlaneStrainState::db_density + 1;
+
+const int contrib::materials::PlaneStrainState::db_vp =
+ contrib::materials::PlaneStrainState::db_vs + 1;
+
+// Indices of state variables in the state variables array.
+const int contrib::materials::MaxwellIsotropic3D::s_totalStrain = 0;
+
+const int contrib::materials::MaxwellIsotropic3D::s_stress =
+ contrib::materials::PlaneStrainState::s_totalStrain +
+ contrib::materials::_PlaneStrainState::tensorSize;
+
+// ----------------------------------------------------------------------
+// Default constructor.
+contrib::materials::PlaneStrainState::PlaneStrainState(void) :
+ pylith::materials::ElasticMaterial(_PlaneStrainState::dimension,
+ _PlaneStrainState::tensorSize,
+ _PlaneStrainState::numElasticConsts,
+ pylith::materialsMetadata(_PlaneStrainState::properties,
+ _PlaneStrainState::numProperties,
+ _PlaneStrainState::dbProperties,
+ _PlaneStrainState::numDBProperties,
+ _PlaneStrainState::stateVars,
+ _PlaneStrainState::numStateVars,
+ 0, 0))
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+contrib::materials::PlaneStrainState::~PlaneStrainState(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Compute parameters from values in spatial database.
+void
+contrib::materials::PlaneStrainState::_dbToProperties(
+ double* const propValues,
+ const pylith::double_array& dbValues) const
+{ // _dbToProperties
+ // Check consistency of arguments
+ assert(0 != propValues);
+ const int numDBValues = dbValues.size();
+ assert(_PlaneStrainState::numDBProperties == numDBValues);
+
+ // Extract values from array using our defined indices.
+ const double density = dbValues[db_density];
+ const double vs = dbValues[db_vs];
+ const double vp = dbValues[db_vp];
+
+ // Check for reasonable values. If user supplied unreasonable values
+ // throw an exception.
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0) {
+ std::ostringstream msg;
+ msg << "Spatial database returned nonpositive value for physical "
+ << "properties.\n"
+ << "density: " << density << "\n"
+ << "vp: " << vp << "\n"
+ << "vs: " << vs << "\n";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ // Compute physical properties that we store from the user-supplied
+ // physical properties.
+ const double mu = density * vs*vs;
+ const double lambda = density * vp*vp - 2.0*mu;
+
+ // Check for reasonable values. If values are unreasonable throw an
+ // exception.
+ if (lambda <= 0.0) {
+ std::ostringstream msg;
+ msg << "Attempted to set Lame's constant lambda to nonpositive value.\n"
+ << "density: " << density << "\n"
+ << "vp: " << vp << "\n"
+ << "vs: " << vs << "\n";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ // Store computed physical properties in the properties array.
+ propValues[p_density] = density;
+ propValues[p_mu] = mu;
+ propValues[p_lambda] = lambda;
+
+ // PetscLogFlops(6); // Add to logging of number of floating point operations
+} // _dbToProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+contrib::materials::PlaneStrainState::_nondimProperties(double* const values,
+ const int nvalues) const
+{ // _nondimProperties
+ // Check consistency of arguments.
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _PlaneStrainState::numProperties);
+
+ // Get scales needed to nondimensional parameters from the
+ // Nondimensional object.
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+
+ // Use the Nondimensional::nondimensionalize() function to
+ // nondimensionalize the quantities using the appropriate scale.
+ values[p_density] =
+ _normalizer->nondimensionalize(values[p_density], densityScale);
+ values[p_mu] =
+ _normalizer->nondimensionalize(values[p_mu], pressureScale);
+ values[p_lambda] =
+ _normalizer->nondimensionalize(values[p_lambda], pressureScale);
+
+ // PetscLogFlops(3); // Add to logging of number of floating point operations
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+contrib::materials::PlaneStrainState::_dimProperties(double* const values,
+ const int nvalues) const
+{ // _dimProperties
+ // Check consistency of arguments
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _PlaneStrainState::numProperties);
+
+ // Get scales needed to dimensional parameters from the
+ // Nondimensional object.
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+
+ // Use the Nondimensional::dimensionalize() function to
+ // dimensionalize the quantities using the appropriate scale.
+ values[p_density] =
+ _normalizer->dimensionalize(values[p_density], densityScale);
+ values[p_mu] =
+ _normalizer->dimensionalize(values[p_mu], pressureScale);
+ values[p_lambda] =
+ _normalizer->dimensionalize(values[p_lambda], pressureScale);
+
+ // PetscLogFlops(3); // Add to logging of number of floating point operations
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Compute density at location from properties.
+void
+contrib::materials::PlaneStrainState::_calcDensity(double* const density,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars)
+{ // calcDensity
+ // Check consistency of arguments.
+ assert(0 != density);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 == numStateVars);
+
+ // Set density using physical properties (trivial since one our
+ // physical properties is density).
+ density[0] = properties[p_density];
+} // calcDensity
+
+// ----------------------------------------------------------------------
+// Compute stress tensor at location from properties.
+void
+contrib::materials::PlaneStrainState::_calcStress(double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
+ const bool computeStateVars)
+{ // _calcStress
+ // Check consistency of arguments.
+ assert(0 != stress);
+ assert(_PlaneStrainState::tensorSize == stressSize);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 == numStateVars);
+ assert(0 != totalStrain);
+ assert(_PlaneStrainState::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_PlaneStrainState::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_PlaneStrainState::tensorSize == initialStrainSize);
+
+ // Extract the material properties from the properties array.
+ const double density = properties[p_density];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
+
+ // Compute a convenient constant.
+ const double mu2 = 2.0*mu;
+
+ // Compute the current strains accounting for the initial strain.
+ const double e11 = totalStrain[0] - initialStrain[0];
+ const double e22 = totalStrain[1] - initialStrain[1];
+ const double e12 = totalStrain[2] - initialStrain[2];
+
+ // Compute another convenient constant.
+ const double s12 = lambda * (e11 + e22);
+
+ // Compute the stresses and store them in the stress array.
+ stress[0] = s12 + mu2*e11 + initialStress[0];
+ stress[1] = s12 + mu2*e22 + initialStress[1];
+ stress[2] = mu2 * e12 + initialStress[2];
+
+ // PetscLogFlops(14); // Add to logging of number of floating point operations
+} // _calcStress
+
+// ----------------------------------------------------------------------
+// Compute elastic constants at location from properties.
+void
+contrib::materials::PlaneStrainState::_calcElasticConsts(
+ double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
+{ // calcElasticConsts
+ // Check consistency of arguments.
+ assert(0 != elasticConsts);
+ assert(_PlaneStrainState::numElasticConsts == numElasticConsts);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 == numStateVars);
+ assert(0 != totalStrain);
+ assert(_PlaneStrainState::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_PlaneStrainState::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_PlaneStrainState::tensorSize == initialStrainSize);
+
+ // Extract the material properties from the properties array.
+ const double density = properties[p_density];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
+
+ // Compute a couple convenient constants.
+ const double mu2 = 2.0 * mu;
+ const double lambda2mu = lambda + mu2;
+
+ // Compute the elastic constants and store them in the elastic
+ // constants array.
+ elasticConsts[0] = lambda2mu; // C1111
+ elasticConsts[1] = lambda; // C1122
+ elasticConsts[2] = 0; // C1112
+ elasticConsts[3] = lambda2mu; // C2222
+ elasticConsts[4] = 0; // C2212
+ elasticConsts[5] = mu2; // C1212
+
+ // PetscLogFlops(2); // Add to logging of number of floating point operations
+} // calcElasticConsts
+
+// ----------------------------------------------------------------------
+// Get stable time step for implicit time integration.
+double
+contrib::materials::PlaneStrainState::stableTimeStepImplicit(
+ const topology::Mesh& mesh) {
+ // Override the ElasticMaterial::stableTimeStepImplicit() function
+ // (which calls _stableTimeStepImplicit() for each quadrature point
+ // ) with an optimized calculation of the stable time step. This is
+ // possible because the stable time step for an elastic material is
+ // infinite, so we can simply return a very large number.
+ return pylith::PYLITH_MAXDOUBLE;
+}
+
+// ----------------------------------------------------------------------
+// Get stable time step for implicit time integration.
+double
+contrib::materials::PlaneStrainState::_stableTimeStepImplicit(
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars) const
+{ // _stableTimeStepImplicit Return the stable time step for this
+ // material given its current state. This function will never be
+ // called because we provide the stableTimeStepImplicit() function,
+ // but we implement this function to satisfy the requirements of an
+ // interface for an elastic material (which is defined by the
+ // Material and ElasticMaterial objects).
+ return pylith::PYLITH_MAXDOUBLE;
+} // _stableTimeStepImplicit
+
+// ----------------------------------------------------------------------
+// Update state variables.
+void
+contrib::materials::PlaneStrainState::_updateStateVars(
+ double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
+{ // _updateStateVars
+ // Check consistency of arguments.
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != totalStrain);
+ assert(_PlaneStrainState::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_PlaneStrainState::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_PlaneStrainState::tensorSize == initialStrainSize);
+
+ // Store the tensor size as a local value.
+ const int tensorSize = _tensorSize;
+
+ // Store the total strain in the state variable array.
+ for (int iComp=0; iComp < tensorSize; ++iComp)
+ stateVars[s_totalStrain+iComp] = totalStrain[iComp];
+
+ // Use the _calcStress function to compute the stresses.
+ const bool computeStateVars = false; // We are computing the state vars here!
+ _calcStress(&stateVars[s_stress], tensorSize,
+ properties, numProperties,
+ stateVars, numStateVars,
+ totalStrain, strainSize,
+ initialStress, initialStressSize,
+ initialStrain, initialStrainSize,
+ computeStateVars);
+} // _updateStateVars
+
+
+// End of file
Added: short/3D/PyLith/trunk/templates/materials/PlaneStrainState.hh
===================================================================
--- short/3D/PyLith/trunk/templates/materials/PlaneStrainState.hh (rev 0)
+++ short/3D/PyLith/trunk/templates/materials/PlaneStrainState.hh 2009-06-15 22:54:22 UTC (rev 15282)
@@ -0,0 +1,262 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/* @brief C++ PlaneStrainState object that stores stress and strain as
+ * state variables.
+ *
+ * This objects demonstrates how to extend PyLith by adding bulk
+ * constitutive models. This PlaneStrainState object is identical to
+ * the ElasticPlaneStrain object in PyLith except that it stores the
+ * current stress and strain as state variables.
+ *
+ * 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
+ * related to the elasticity constants used in the finite-element
+ * integration.
+ *
+ * $\sigma - \sigma_0 = C (\epsilon - \epsilon_0)
+ *
+ * This implies that when $\epsilon = \epsilon_0$, $\sigma =
+ * \sigma_0$.
+ */
+
+#if !defined(pylith_materials_planestrainstate_hh)
+#define pylith_materials_planestrainstate_hh
+
+#include "pylith/materials/ElasticMaterial.hh" // ISA ElasticMaterial
+
+namespace contrib {
+ namespace materials {
+ class PlaneStrainState;
+ } // materials
+} // contrib
+
+// PlaneStrainState -----------------------------------------------------
+class contrib::materials::PlaneStrainState :
+ public pylith::materials::ElasticMaterial
+{ // class PlaneStrainState
+ friend class TestPlaneStrainState; // unit testing
+
+ // --------------------------------------------------------------------
+ // All of these functions are required to satisfy the PyLith
+ // interface for a bulk constitutive model.
+ // --------------------------------------------------------------------
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor
+ PlaneStrainState(void);
+
+ /// Destructor
+ ~PlaneStrainState(void);
+
+ /** Get stable time step for implicit time integration.
+ *
+ * Default is MAXDOUBLE (or 1.0e+30 if MAXFLOAT is not defined in math.h).
+ *
+ * This function is optional but provides an optimized
+ * implementation of the more general
+ * ElasticMaterial::stableTimeStepImplicit().
+ *
+ * @param mesh Finite-element mesh.
+ * @returns Time step
+ */
+ double stableTimeStepImplicit(const topology::Mesh& mesh);
+
+ // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+ /** Compute properties from values in spatial database.
+ *
+ * Order of values in arrays matches order used in dbValues() and
+ * parameterNames().
+ *
+ * @param propValues Array of property values.
+ * @param dbValues Array of database values.
+ */
+ void _dbToProperties(double* const propValues,
+ const pylith::double_array& dbValues) const;
+
+ /** Nondimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _nondimProperties(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _dimProperties(double* const values,
+ const int nvalues) const;
+
+ /** Compute density from properties.
+ *
+ * @param density Array for density.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ */
+ void _calcDensity(double* const density,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars);
+
+ /** Compute stress tensor from properties and state variables. If
+ * the state variables are from the previous time step, then the
+ * computeStateVars flag should be set to true so that the state
+ * variables are updated (but not stored) when computing the
+ * stresses.
+ *
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ * @param computeStateVars Flag indicating to compute updated state variables.
+ */
+ void _calcStress(double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
+ const bool computeStateVars);
+
+ /** Compute derivatives of elasticity matrix from properties.
+ *
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ */
+ void _calcElasticConsts(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ /** Get stable time step for implicit time integration.
+ *
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ *
+ * @returns Time step
+ */
+ double _stableTimeStepImplicit(const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars) const;
+
+ // --------------------------------------------------------------------
+ // Optional function in the PyLith interface for a bulk constitutive
+ // model. Even though this function is optional, for it to be used
+ // it the interface must exactly matched the one specified in
+ // ElasticMaterial.
+ // --------------------------------------------------------------------
+
+ // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+ /** Update state variables (for next time step).
+ *
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ */
+ void _updateStateVars(double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+ // --------------------------------------------------------------------
+ // We use these constants for consistent access into the arrays of
+ // physical properties and state variables.
+ // --------------------------------------------------------------------
+
+ static const int p_density;
+ static const int p_mu;
+ static const int p_lambda;
+ static const int db_density;
+ static const int db_vs;
+ static const int db_vp;
+
+ static const int s_totalStrain;
+ static const int s_stress;
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ PlaneStrainState(const PlaneStrainState& m);
+
+ /// Not implemented
+ const PlaneStrainState& operator=(const PlaneStrainState& m);
+
+}; // class PlaneStrainState
+
+#endif // pylith_materials_planestrainstate_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/templates/materials/PlaneStrainState.i
===================================================================
--- short/3D/PyLith/trunk/templates/materials/PlaneStrainState.i (rev 0)
+++ short/3D/PyLith/trunk/templates/materials/PlaneStrainState.i 2009-06-15 22:54:22 UTC (rev 15282)
@@ -0,0 +1,206 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+// SWIG interface to C++ UniformVelModel object.
+
+/* This is nearly identical to the C++ PlaneStrainState header
+ * file. There are a few important differences required by SWIG:
+ *
+ * (1) Instead of forward declaring the PlaneStrainState class, we
+ * embed the class definition within the namespace declarations.
+ *
+ * (2) We only include public members and methods and implementations
+ * of abstract methods, because this is an interface file.
+ */
+
+namespace contrib {
+ namespace materials {
+
+ class PlaneStrainState : public pylith::materials::ElasticMaterial
+ { // class PlaneStrainState
+
+ // PUBLIC METHODS /////////////////////////////////////////////////
+ public :
+
+ /// Default constructor
+ PlaneStrainState(void);
+
+ /// Destructor
+ ~PlaneStrainState(void);
+
+ /** Get stable time step for implicit time integration.
+ *
+ * Default is MAXDOUBLE (or 1.0e+30 if MAXFLOAT is not defined
+ * in math.h).
+ *
+ * This function is optional but provides an optimized
+ * implementation of the more general
+ * ElasticMaterial::stableTimeStepImplicit().
+ *
+ * @param mesh Finite-element mesh.
+ * @returns Time step
+ */
+ double stableTimeStepImplicit(const topology::Mesh& mesh);
+
+ // PROTECTED METHODS //////////////////////////////////////////////
+ protected :
+
+ /** Compute properties from values in spatial database.
+ *
+ * Order of values in arrays matches order used in dbValues() and
+ * parameterNames().
+ *
+ * @param propValues Array of property values.
+ * @param dbValues Array of database values.
+ */
+ void _dbToProperties(double* const propValues,
+ const pylith::double_array& dbValues) const;
+
+ /** Nondimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _nondimProperties(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ void _dimProperties(double* const values,
+ const int nvalues) const;
+
+ /** Compute density from properties.
+ *
+ * @param density Array for density.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ */
+ void _calcDensity(double* const density,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars);
+
+ /** Compute stress tensor from properties and state variables. If
+ * the state variables are from the previous time step, then the
+ * computeStateVars flag should be set to true so that the state
+ * variables are updated (but not stored) when computing the
+ * stresses.
+ *
+ * @param stress Array for stress tensor.
+ * @param stressSize Size of stress tensor.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ * @param computeStateVars Flag indicating to compute updated
+ * state variables.
+ */
+ void _calcStress(double* const stress,
+ const int stressSize,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
+ const bool computeStateVars);
+
+ /** Compute derivatives of elasticity matrix from properties.
+ *
+ * @param elasticConsts Array for elastic constants.
+ * @param numElasticConsts Number of elastic constants.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ */
+ void _calcElasticConsts(double* const elasticConsts,
+ const int numElasticConsts,
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ /** Get stable time step for implicit time integration.
+ *
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ *
+ * @returns Time step
+ */
+ double _stableTimeStepImplicit(const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars) const;
+
+ // PROTECTED METHODS //////////////////////////////////////////////
+ protected :
+
+ /** Update state variables (for next time step).
+ *
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param totalStrain Total strain at location.
+ * @param strainSize Size of strain tensor.
+ * @param initialStress Initial stress tensor at location.
+ * @param initialStressSize Size of initial stress array.
+ * @param initialStrain Initial strain tensor at location.
+ * @param initialStrainSize Size of initial strain array.
+ */
+ void _updateStateVars(double* const stateVars,
+ const int numStateVars,
+ const double* properties,
+ const int numProperties,
+ const double* totalStrain,
+ const int strainSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize);
+
+ }; // class PlaneStrainState
+
+ } // materials
+} // contrib
+
+
+// End of file
Added: short/3D/PyLith/trunk/templates/materials/PlaneStrainState.py
===================================================================
--- short/3D/PyLith/trunk/templates/materials/PlaneStrainState.py (rev 0)
+++ short/3D/PyLith/trunk/templates/materials/PlaneStrainState.py 2009-06-15 22:54:22 UTC (rev 15282)
@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @brief Python object implementing 1-D isotropic linear elastic
+## material for plane strain.
+##
+## Factory: material.
+
+from ElasticMaterial import ElasticMaterial # ISA ElasticMaterial
+
+# Import the SWIG module PlanseStrainState object and rename it
+# ModulePlaneStrainState so that it doesn't clash with the local
+# Python class with the same name.
+from materials import PlaneStrainState as ModulePlaneStrainState
+
+# PlaneStrainState class
+class PlaneStrainState(ElasticMaterial, ModulePlaneStrainState):
+ """
+ Python object implementing 2-D isotropic linear elastic material for
+ plane strain.
+
+ Factory: material.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="planestrainstate"):
+ """
+ Constructor.
+ """
+ ElasticMaterial.__init__(self, name)
+ # Set the fields that are available for output. These are the
+ # stored physical properties, state variables, and the total
+ # strain tensor and the stress tensor. For bulk elasticity
+ # materials we can compute the stresses and strains in a general
+ # fashion, so they need not be stored as they are in this example.
+ #
+ # There are no vertex fields because the constitutive model
+ # operations on quantities evaluated at the quadrature points.
+ self.availableFields = \
+ {'vertex': \
+ {'info': [],
+ 'data': []},
+ 'cell': \
+ {'info': ["mu", "lambda", "density"],
+ 'data': ["total_strain", "stress"]}}
+ self._loggingPrefix = "MaPlSn " # Prefix that appears in PETSc logging
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _createModuleObj(self):
+ """
+ Call constructor for module object for access to C++ object.
+ """
+ ModulePlaneStrainState.__init__(self)
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def material(): # The name of this function MUST be 'material'.
+ """
+ Factory associated with PlaneStrainState.
+ """
+ return PlaneStrainState() # Return our object
+
+
+# End of file
Added: short/3D/PyLith/trunk/templates/materials/README
===================================================================
--- short/3D/PyLith/trunk/templates/materials/README (rev 0)
+++ short/3D/PyLith/trunk/templates/materials/README 2009-06-15 22:54:22 UTC (rev 15282)
@@ -0,0 +1,12 @@
+PyLith interface files are installed to $PREFIX/shared/pylith/swig.
+
+Functions in bulk material routines are called for every quadrature
+point of every cell of the material, so efficient code is essential.
+
+The consistency checks of user input should throw exceptions with a
+meaningful error message so that the user can fix the
+mistake. Consistency checks in the other routines should use the
+assert() macro so that they can be removed during compilation. To
+remove the assert() calls during compilation configure with -DNDEBUG
+added to the CFLAGS and CXXFLAGS environment variables.
+
Added: short/3D/PyLith/trunk/templates/materials/__init__.py
===================================================================
--- short/3D/PyLith/trunk/templates/materials/__init__.py (rev 0)
+++ short/3D/PyLith/trunk/templates/materials/__init__.py 2009-06-15 22:54:22 UTC (rev 15282)
@@ -0,0 +1,21 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file __init__.py
+##
+## @brief Python module initialization
+
+__all__ = ['PlaneStrainState',
+ ]
+
+
+# End of file
Added: short/3D/PyLith/trunk/templates/materials/configure.ac
===================================================================
--- short/3D/PyLith/trunk/templates/materials/configure.ac (rev 0)
+++ short/3D/PyLith/trunk/templates/materials/configure.ac 2009-06-15 22:54:22 UTC (rev 15282)
@@ -0,0 +1,79 @@
+# -*- autoconf -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+AC_PREREQ(2.59)
+AC_INIT([materialscontrib], [0.0.1], [cig-short at geodynamics.org])
+AC_CONFIG_HEADER([portinfo])
+AC_CONFIG_AUX_DIR([./aux-config])
+#AC_CONFIG_MACRO_DIR([./m4])
+AM_INIT_AUTOMAKE([foreign])
+
+# ----------------------------------------------------------------------
+# C/C++/libtool/install
+AC_PROG_CXX
+AC_PROG_CC
+AC_DISABLE_STATIC
+
+AC_PROG_LIBTOOL
+if test "$allow_undefined_flag" = unsupported; then
+ # See issue119.
+ AM_LDFLAGS="-no-undefined $AM_LDFLAGS"
+fi
+AM_CONDITIONAL([NO_UNDEFINED], [test "$allow_undefined_flag" = unsupported])
+AC_SUBST(AM_LDFLAGS)
+
+AC_PROG_INSTALL
+
+# PYTHON
+AM_PATH_PYTHON([2.3])
+CIT_PYTHON_SYSCONFIG
+
+# SWIG
+CIT_NUMPY_PYTHON_MODULE
+CIT_NUMPY_INCDIR
+AC_PROG_SWIG(1.3.33)
+
+# PYLITH
+AC_LANG(C++)
+AC_CHECK_HEADER([pylith/materials/ElasticMaterial.hh], [], [
+ AC_MSG_ERROR([PyLith ElasticMaterial header not found; try CPPFLAGS="-I<PyLith include dir>"])
+])
+AC_MSG_CHECKING([for materials::ElasticIsotropic3D in -lpylith])
+AC_REQUIRE_CPP
+AC_COMPILE_IFELSE(
+ [AC_LANG_PROGRAM([[#include <pylith/materials/ElasticIsotropic3D.hh>]],
+ [[pylith::materials::ElasticIsotropic3D material;]])],
+ [AC_MSG_RESULT(yes)],
+ [AC_MSG_RESULT(no)
+ AC_MSG_ERROR([PyLith library not found; try LDFLAGS="-L<PyLith lib dir>"])
+])
+if test "x$PYLITH_SWIG_DIR" != "x" ; then
+ AC_CHECK_FILE([$PYLITH_SWIG_DIR/materials/ElasticMaterial.i], [], [
+ AC_MSG_ERROR([PyLith ElasticMaterial.i SWIG interface file not found])])
+ AC_SUBST([SWIG_FLAGS], ["-I$PYLITH_SWIG_DIR $SWIG_FLAGS"])
+else
+ AC_CHECK_FILE([materials/ElasticMaterial.i], [], [
+ AC_MSG_ERROR([PyLith ElasticMaterial.i SWIG interface file not found; Try setting PYLITH_SWIG_DIR=<directory containing materials/ElasticMaterial.i>])])
+fi
+
+
+# ENDIANNESS
+AC_C_BIGENDIAN
+
+# ----------------------------------------------------------------------
+AC_CONFIG_FILES([Makefile
+ tests/Makefile])
+
+AC_OUTPUT
+
+
+# End of file
Added: short/3D/PyLith/trunk/templates/materials/materialscontrib.i
===================================================================
--- short/3D/PyLith/trunk/templates/materials/materialscontrib.i (rev 0)
+++ short/3D/PyLith/trunk/templates/materials/materialscontrib.i 2009-06-15 22:54:22 UTC (rev 15282)
@@ -0,0 +1,52 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// Define the spatialdatacontrib SWIG interface module.
+
+// Set module name
+%module materialscontrib
+
+// Header files for module C++ code.
+%{
+#include "PlaneStrainState.hh"
+%}
+
+// Convert standard C++ exceptions to Python exceptions.
+%include "exception.i"
+%exception {
+ try {
+ $action
+ } catch (const std::exception& err) {
+ SWIG_exception(SWIG_RuntimeError, err.what());
+ } // try/catch
+} // exception
+
+%include "typemaps.i"
+%include "include/doublearray.i"
+
+// Numpy interface stuff
+%{
+#define SWIG_FILE_WITH_INIT
+%}
+%include "include/numpy.i"
+%init %{
+import_array();
+%}
+
+
+// Interface files.
+%include "materials/Material.i"
+%include "materials/ElasticMaterial.i"
+%include "PlaneStrainState.i"
+
+
+// End of file
More information about the CIG-COMMITS
mailing list