[cig-commits] r13439 - short/3D/PyLith/trunk/libsrc/materials
brad at geodynamics.org
brad at geodynamics.org
Tue Dec 2 12:06:36 PST 2008
Author: brad
Date: 2008-12-02 12:06:35 -0800 (Tue, 02 Dec 2008)
New Revision: 13439
Modified:
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/ElasticPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh
short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh
short/3D/PyLith/trunk/libsrc/materials/Material.cc
short/3D/PyLith/trunk/libsrc/materials/Material.hh
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
Log:
Nondimensionzlized physical properties.
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2008-12-02 20:06:35 UTC (rev 13439)
@@ -21,7 +21,7 @@
#include "petsc.h" // USES PetscLogFlops
-#include <assert.h> // USES assert()
+#include <cassert> // USES assert()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
@@ -107,7 +107,8 @@
if (density <= 0.0 || vs <= 0.0 || vp <= 0.0) {
std::ostringstream msg;
- msg << "Spatial database returned negative value for physical properties.\n"
+ msg << "Spatial database returned nonpositive value for physical "
+ << "properties.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -134,6 +135,88 @@
} // _dbToProperties
// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+pylith::materials::ElasticIsotropic3D::_nondimProperties(double* const values,
+ const int nvalues) const
+{ // _nondimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticIsotropic3D::numProperties);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ values[_ElasticIsotropic3D::pidDensity] =
+ _normalizer->nondimensionalize(values[_ElasticIsotropic3D::pidDensity],
+ densityScale);
+ values[_ElasticIsotropic3D::pidMu] =
+ _normalizer->nondimensionalize(values[_ElasticIsotropic3D::pidMu],
+ pressureScale);
+ values[_ElasticIsotropic3D::pidLambda] =
+ _normalizer->nondimensionalize(values[_ElasticIsotropic3D::pidLambda],
+ pressureScale);
+
+ PetscLogFlops(3);
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+pylith::materials::ElasticIsotropic3D::_dimProperties(double* const values,
+ const int nvalues) const
+{ // _dimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticIsotropic3D::numProperties);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ values[_ElasticIsotropic3D::pidDensity] =
+ _normalizer->dimensionalize(values[_ElasticIsotropic3D::pidDensity],
+ densityScale);
+ values[_ElasticIsotropic3D::pidMu] =
+ _normalizer->dimensionalize(values[_ElasticIsotropic3D::pidMu],
+ pressureScale);
+ values[_ElasticIsotropic3D::pidLambda] =
+ _normalizer->dimensionalize(values[_ElasticIsotropic3D::pidLambda],
+ pressureScale);
+
+ PetscLogFlops(3);
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize initial state.
+void
+pylith::materials::ElasticIsotropic3D::_nondimInitState(double* const values,
+ const int nvalues) const
+{ // _nondimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticIsotropic3D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _nondimInitState
+
+// ----------------------------------------------------------------------
+// Dimensionalize initial state.
+void
+pylith::materials::ElasticIsotropic3D::_dimInitState(double* const values,
+ const int nvalues) const
+{ // _dimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticIsotropic3D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _dimInitState
+
+// ----------------------------------------------------------------------
// Compute density at location from properties.
void
pylith::materials::ElasticIsotropic3D::_calcDensity(
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.hh 2008-12-02 20:06:35 UTC (rev 13439)
@@ -63,6 +63,38 @@
void _dbToProperties(double* const propValues,
const 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;
+
+ /** Nondimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _nondimInitState(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _dimInitState(double* const values,
+ const int nvalues) const;
+
/** Compute density from properties.
*
* @param density Array for density.
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2008-12-02 20:06:35 UTC (rev 13439)
@@ -18,8 +18,8 @@
#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
#include "pylith/utils/sievetypes.hh" // USES Mesh
-#include <string.h> // USES memcpy()
-#include <assert.h> // USES assert()
+#include <cstring> // USES memcpy()
+#include <cassert> // USES assert()
// ----------------------------------------------------------------------
// Default constructor.
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc 2008-12-02 20:06:35 UTC (rev 13439)
@@ -16,9 +16,11 @@
#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
#include "petsc.h" // USES PetscLogFlops
-#include <assert.h> // USES assert()
+#include <cassert> // USES assert()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
@@ -102,9 +104,10 @@
const double vs = dbValues[_ElasticPlaneStrain::didVs];
const double vp = dbValues[_ElasticPlaneStrain::didVp];
- if (density < 0.0 || vs < 0.0 || vp < 0.0) {
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0) {
std::ostringstream msg;
- msg << "Spatial database returned negative value for physical properties.\n"
+ msg << "Spatial database returned nonpositive value for physical "
+ << "properties.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -114,9 +117,9 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- if (lambda < 0.0) {
+ if (lambda <= 0.0) {
std::ostringstream msg;
- msg << "Attempted to set Lame's constant lambda to negative value.\n"
+ msg << "Attempted to set Lame's constant lambda to nonpositive value.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -131,6 +134,80 @@
} // _dbToProperties
// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+pylith::materials::ElasticPlaneStrain::_nondimProperties(double* const values,
+ const int nvalues) const
+{ // _nondimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticPlaneStrain::numProperties);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ values[_ElasticPlaneStrain::pidDensity] =
+ _normalizer->nondimensionalize(values[_ElasticPlaneStrain::pidDensity],
+ densityScale);
+ values[_ElasticPlaneStrain::pidMu] =
+ _normalizer->nondimensionalize(values[_ElasticPlaneStrain::pidMu],
+ pressureScale);
+ values[_ElasticPlaneStrain::pidLambda] =
+ _normalizer->nondimensionalize(values[_ElasticPlaneStrain::pidLambda],
+ pressureScale);
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+pylith::materials::ElasticPlaneStrain::_dimProperties(double* const values,
+ const int nvalues) const
+{ // _dimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticPlaneStrain::numProperties);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ values[_ElasticPlaneStrain::pidDensity] =
+ _normalizer->dimensionalize(values[_ElasticPlaneStrain::pidDensity],
+ densityScale);
+ values[_ElasticPlaneStrain::pidMu] =
+ _normalizer->dimensionalize(values[_ElasticPlaneStrain::pidMu],
+ pressureScale);
+ values[_ElasticPlaneStrain::pidLambda] =
+ _normalizer->dimensionalize(values[_ElasticPlaneStrain::pidLambda],
+ pressureScale);
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize initial state.
+void
+pylith::materials::ElasticPlaneStrain::_nondimInitState(double* const values,
+ const int nvalues) const
+{ // _nondimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticPlaneStrain::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(values, nvalues, pressureScale);
+} // _nondimInitState
+
+// ----------------------------------------------------------------------
+// Dimensionalize initial state.
+void
+pylith::materials::ElasticPlaneStrain::_dimInitState(double* const values,
+ const int nvalues) const
+{ // _dimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticPlaneStrain::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(values, nvalues, pressureScale);
+} // _dimInitState
+
+// ----------------------------------------------------------------------
// Compute density at location from properties.
void
pylith::materials::ElasticPlaneStrain::_calcDensity(
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.hh 2008-12-02 20:06:35 UTC (rev 13439)
@@ -63,6 +63,38 @@
void _dbToProperties(double* const propValues,
const 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;
+
+ /** Nondimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _nondimInitState(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _dimInitState(double* const values,
+ const int nvalues) const;
+
/** Compute density from properties.
*
* @param density Array for density.
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc 2008-12-02 20:06:35 UTC (rev 13439)
@@ -16,9 +16,11 @@
#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
#include "petsc.h" // USES PetscLogFlops
-#include <assert.h> // USES assert()
+#include <cassert> // USES assert()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
@@ -100,9 +102,10 @@
const double vs = dbValues[_ElasticPlaneStress::didVs];
const double vp = dbValues[_ElasticPlaneStress::didVp];
- if (density < 0.0 || vs < 0.0 || vp < 0.0) {
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0) {
std::ostringstream msg;
- msg << "Spatial database returned negative value for physical properties.\n"
+ msg << "Spatial database returned nonpositive value for physical "
+ << "properties.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -112,9 +115,9 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- if (lambda < 0.0) {
+ if (lambda <= 0.0) {
std::ostringstream msg;
- msg << "Attempted to set Lame's constant lambda to negative value.\n"
+ msg << "Attempted to set Lame's constant lambda to nonpositive value.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -129,6 +132,88 @@
} // _dbToProperties
// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+pylith::materials::ElasticPlaneStress::_nondimProperties(double* const values,
+ const int nvalues) const
+{ // _nondimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticPlaneStress::numProperties);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ values[_ElasticPlaneStress::pidDensity] =
+ _normalizer->nondimensionalize(values[_ElasticPlaneStress::pidDensity],
+ densityScale);
+ values[_ElasticPlaneStress::pidMu] =
+ _normalizer->nondimensionalize(values[_ElasticPlaneStress::pidMu],
+ pressureScale);
+ values[_ElasticPlaneStress::pidLambda] =
+ _normalizer->nondimensionalize(values[_ElasticPlaneStress::pidLambda],
+ pressureScale);
+
+ PetscLogFlops(3);
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+pylith::materials::ElasticPlaneStress::_dimProperties(double* const values,
+ const int nvalues) const
+{ // _dimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticPlaneStress::numProperties);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ values[_ElasticPlaneStress::pidDensity] =
+ _normalizer->dimensionalize(values[_ElasticPlaneStress::pidDensity],
+ densityScale);
+ values[_ElasticPlaneStress::pidMu] =
+ _normalizer->dimensionalize(values[_ElasticPlaneStress::pidMu],
+ pressureScale);
+ values[_ElasticPlaneStress::pidLambda] =
+ _normalizer->dimensionalize(values[_ElasticPlaneStress::pidLambda],
+ pressureScale);
+
+ PetscLogFlops(3);
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize initial state.
+void
+pylith::materials::ElasticPlaneStress::_nondimInitState(double* const values,
+ const int nvalues) const
+{ // _nondimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticPlaneStress::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _nondimInitState
+
+// ----------------------------------------------------------------------
+// Dimensionalize initial state.
+void
+pylith::materials::ElasticPlaneStress::_dimInitState(double* const values,
+ const int nvalues) const
+{ // _dimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticPlaneStress::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _dimInitState
+
+// ----------------------------------------------------------------------
// Compute density at location from properties.
void
pylith::materials::ElasticPlaneStress::_calcDensity(double* const density,
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.hh 2008-12-02 20:06:35 UTC (rev 13439)
@@ -63,6 +63,38 @@
void _dbToProperties(double* const propertyValues,
const 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;
+
+ /** Nondimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _nondimInitState(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _dimInitState(double* const values,
+ const int nvalues) const;
+
/** Compute density from properties.
*
* @param density Array for density.
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.cc 2008-12-02 20:06:35 UTC (rev 13439)
@@ -16,9 +16,11 @@
#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
#include "petsc.h" // USES PetscLogFlops
-#include <assert.h> // USES assert()
+#include <cassert> // USES assert()
// ----------------------------------------------------------------------
namespace pylith {
@@ -98,9 +100,10 @@
const double vp = dbValues[_ElasticStrain1D::didVp];
const double vs = dbValues[_ElasticStrain1D::didVs];
- if (density < 0.0 || vs < 0.0 || vp < 0.0) {
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0) {
std::ostringstream msg;
- msg << "Spatial database returned negative value for physical properties.\n"
+ msg << "Spatial database returned nonpositive value for physical "
+ << "properties.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -110,9 +113,9 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- if (lambda < 0.0) {
+ if (lambda <= 0.0) {
std::ostringstream msg;
- msg << "Attempted to set Lame's constant lambda to negative value.\n"
+ msg << "Attempted to set Lame's constant lambda to nonpositive value.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -127,6 +130,88 @@
} // _dbToProperties
// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+pylith::materials::ElasticStrain1D::_nondimProperties(double* const values,
+ const int nvalues) const
+{ // _nondimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticStrain1D::numProperties);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ values[_ElasticStrain1D::pidDensity] =
+ _normalizer->nondimensionalize(values[_ElasticStrain1D::pidDensity],
+ densityScale);
+ values[_ElasticStrain1D::pidMu] =
+ _normalizer->nondimensionalize(values[_ElasticStrain1D::pidMu],
+ pressureScale);
+ values[_ElasticStrain1D::pidLambda] =
+ _normalizer->nondimensionalize(values[_ElasticStrain1D::pidLambda],
+ pressureScale);
+
+ PetscLogFlops(3);
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+pylith::materials::ElasticStrain1D::_dimProperties(double* const values,
+ const int nvalues) const
+{ // _dimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticStrain1D::numProperties);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ values[_ElasticStrain1D::pidDensity] =
+ _normalizer->dimensionalize(values[_ElasticStrain1D::pidDensity],
+ densityScale);
+ values[_ElasticStrain1D::pidMu] =
+ _normalizer->dimensionalize(values[_ElasticStrain1D::pidMu],
+ pressureScale);
+ values[_ElasticStrain1D::pidLambda] =
+ _normalizer->dimensionalize(values[_ElasticStrain1D::pidLambda],
+ pressureScale);
+
+ PetscLogFlops(3);
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize initial state.
+void
+pylith::materials::ElasticStrain1D::_nondimInitState(double* const values,
+ const int nvalues) const
+{ // _nondimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticStrain1D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _nondimInitState
+
+// ----------------------------------------------------------------------
+// Dimensionalize initial state.
+void
+pylith::materials::ElasticStrain1D::_dimInitState(double* const values,
+ const int nvalues) const
+{ // _dimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticStrain1D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _dimInitState
+
+// ----------------------------------------------------------------------
// Compute density at location from properties.
void
pylith::materials::ElasticStrain1D::_calcDensity(double* const density,
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStrain1D.hh 2008-12-02 20:06:35 UTC (rev 13439)
@@ -62,6 +62,38 @@
void _dbToProperties(double* const propValues,
const 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;
+
+ /** Nondimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _nondimInitState(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _dimInitState(double* const values,
+ const int nvalues) const;
+
/** Compute density from properties.
*
* @param density Array for density.
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.cc 2008-12-02 20:06:35 UTC (rev 13439)
@@ -16,9 +16,11 @@
#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
#include "petsc.h" // USES PetscLogFlops
-#include <assert.h> // USES assert()
+#include <cassert> // USES assert()
// ----------------------------------------------------------------------
namespace pylith {
@@ -97,9 +99,10 @@
const double vs = dbValues[_ElasticStress1D::didVs];
const double vp = dbValues[_ElasticStress1D::didVp];
- if (density < 0.0 || vs < 0.0 || vp < 0.0) {
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0) {
std::ostringstream msg;
- msg << "Spatial database returned negative value for physical properties.\n"
+ msg << "Spatial database returned nonpositive value for physical "
+ << "properties.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -109,9 +112,9 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- if (lambda < 0.0) {
+ if (lambda <= 0.0) {
std::ostringstream msg;
- msg << "Attempted to set Lame's constant lambda to negative value.\n"
+ msg << "Attempted to set Lame's constant lambda to nonpositive value.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -126,6 +129,88 @@
} // _dbToProperties
// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+pylith::materials::ElasticStress1D::_nondimProperties(double* const values,
+ const int nvalues) const
+{ // _nondimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticStress1D::numProperties);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ values[_ElasticStress1D::pidDensity] =
+ _normalizer->nondimensionalize(values[_ElasticStress1D::pidDensity],
+ densityScale);
+ values[_ElasticStress1D::pidMu] =
+ _normalizer->nondimensionalize(values[_ElasticStress1D::pidMu],
+ pressureScale);
+ values[_ElasticStress1D::pidLambda] =
+ _normalizer->nondimensionalize(values[_ElasticStress1D::pidLambda],
+ pressureScale);
+
+ PetscLogFlops(3);
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+pylith::materials::ElasticStress1D::_dimProperties(double* const values,
+ const int nvalues) const
+{ // _dimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticStress1D::numProperties);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ values[_ElasticStress1D::pidDensity] =
+ _normalizer->dimensionalize(values[_ElasticStress1D::pidDensity],
+ densityScale);
+ values[_ElasticStress1D::pidMu] =
+ _normalizer->dimensionalize(values[_ElasticStress1D::pidMu],
+ pressureScale);
+ values[_ElasticStress1D::pidLambda] =
+ _normalizer->dimensionalize(values[_ElasticStress1D::pidLambda],
+ pressureScale);
+
+ PetscLogFlops(3);
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize initial state.
+void
+pylith::materials::ElasticStress1D::_nondimInitState(double* const values,
+ const int nvalues) const
+{ // _nondimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticStress1D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _nondimInitState
+
+// ----------------------------------------------------------------------
+// Dimensionalize initial state.
+void
+pylith::materials::ElasticStress1D::_dimInitState(double* const values,
+ const int nvalues) const
+{ // _dimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _ElasticStress1D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _dimInitState
+
+// ----------------------------------------------------------------------
// Compute density at location from properties.
void
pylith::materials::ElasticStress1D::_calcDensity(double* const density,
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticStress1D.hh 2008-12-02 20:06:35 UTC (rev 13439)
@@ -63,6 +63,38 @@
void _dbToProperties(double* const propValues,
const 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;
+
+ /** Nondimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _nondimInitState(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _dimInitState(double* const values,
+ const int nvalues) const;
+
/** Compute density from properties.
*
* @param density Array for density.
Modified: short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc 2008-12-02 20:06:35 UTC (rev 13439)
@@ -17,12 +17,14 @@
#include "ViscoelasticMaxwell.hh" // USES computeVisStrain
#include "pylith/utils/array.hh" // USES double_array
-#include "pylith/utils/constdefs.h" // USES double_array
+#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
#include "petsc.h" // USES PetscLogFlops
-#include <assert.h> // USES assert()
-#include <string.h> // USES memcpy()
+#include <cassert> // USES assert()
+#include <cstring> // USES memcpy()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
@@ -137,9 +139,10 @@
const double vp = dbValues[_GenMaxwellIsotropic3D::didVp];
const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
- if (density < 0.0 || vs < 0.0 || vp < 0.0) {
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0) {
std::ostringstream msg;
- msg << "Spatial database returned negative value for physical properties.\n"
+ msg << "Spatial database returned nonpositive value for physical "
+ << "properties.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -149,9 +152,9 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- if (lambda < 0.0) {
+ if (lambda <= 0.0) {
std::ostringstream msg;
- msg << "Attempted to set Lame's constant lambda to negative value.\n"
+ msg << "Attempted to set Lame's constant lambda to nonpositive value.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -190,6 +193,96 @@
} // _dbToProperties
// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+pylith::materials::GenMaxwellIsotropic3D::_nondimProperties(double* const values,
+ const int nvalues) const
+{ // _nondimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _totalPropsQuadPt);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ const double timeScale = _normalizer->timeScale();
+ values[_GenMaxwellIsotropic3D::pidDensity] =
+ _normalizer->nondimensionalize(values[_GenMaxwellIsotropic3D::pidDensity],
+ densityScale);
+ values[_GenMaxwellIsotropic3D::pidMuTot] =
+ _normalizer->nondimensionalize(values[_GenMaxwellIsotropic3D::pidMuTot],
+ pressureScale);
+ values[_GenMaxwellIsotropic3D::pidLambdaTot] =
+ _normalizer->nondimensionalize(values[_GenMaxwellIsotropic3D::pidLambdaTot],
+ pressureScale);
+ const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
+ _normalizer->nondimensionalize(&values[_GenMaxwellIsotropic3D::pidMaxwellTime],
+ numMaxwellModels, timeScale);
+
+ PetscLogFlops(3+1*numMaxwellModels);
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+pylith::materials::GenMaxwellIsotropic3D::_dimProperties(double* const values,
+ const int nvalues) const
+{ // _dimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _totalPropsQuadPt);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ const double timeScale = _normalizer->timeScale();
+ values[_GenMaxwellIsotropic3D::pidDensity] =
+ _normalizer->dimensionalize(values[_GenMaxwellIsotropic3D::pidDensity],
+ densityScale);
+ values[_GenMaxwellIsotropic3D::pidMuTot] =
+ _normalizer->dimensionalize(values[_GenMaxwellIsotropic3D::pidMuTot],
+ pressureScale);
+ values[_GenMaxwellIsotropic3D::pidLambdaTot] =
+ _normalizer->dimensionalize(values[_GenMaxwellIsotropic3D::pidLambdaTot],
+ pressureScale);
+ const int numMaxwellModels = _GenMaxwellIsotropic3D::numMaxwellModels;
+ _normalizer->dimensionalize(&values[_GenMaxwellIsotropic3D::pidMaxwellTime],
+ numMaxwellModels, timeScale);
+
+ PetscLogFlops(3+1*numMaxwellModels);
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize initial state.
+void
+pylith::materials::GenMaxwellIsotropic3D::_nondimInitState(double* const values,
+ const int nvalues) const
+{ // _nondimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _GenMaxwellIsotropic3D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _nondimInitState
+
+// ----------------------------------------------------------------------
+// Dimensionalize initial state.
+void
+pylith::materials::GenMaxwellIsotropic3D::_dimInitState(double* const values,
+ const int nvalues) const
+{ // _dimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _GenMaxwellIsotropic3D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _dimInitState
+
+// ----------------------------------------------------------------------
// Compute density at location from properties.
void
pylith::materials::GenMaxwellIsotropic3D::_calcDensity(
Modified: short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.hh 2008-12-02 20:06:35 UTC (rev 13439)
@@ -92,6 +92,38 @@
void _dbToProperties(double* const propValues,
const 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;
+
+ /** Nondimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _nondimInitState(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _dimInitState(double* const values,
+ const int nvalues) const;
+
/** Compute density from properties.
*
* @param density Array for density.
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.cc 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.cc 2008-12-02 20:06:35 UTC (rev 13439)
@@ -18,12 +18,13 @@
#include "pylith/utils/array.hh" // USES double_array, std::vector
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include "pylith/utils/sievetypes.hh" // USES Mesh
-#include <string.h> // USES memcpy()
+#include <cstring> // USES memcpy()
#include <strings.h> // USES strcasecmp()
-#include <assert.h> // USES assert()
+#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
@@ -36,6 +37,7 @@
const PropMetaData* properties,
const int numProperties) :
_dt(0.0),
+ _normalizer(new spatialdata::units::Nondimensional),
_totalPropsQuadPt(0),
_dimension(0),
_tensorSize(tensorSize),
@@ -60,19 +62,31 @@
// Destructor.
pylith::materials::Material::~Material(void)
{ // destructor
+ delete _normalizer; _normalizer = 0;
+
// Python db object owns database, so just set pointer to null
_db = 0;
_initialStateDB = 0;
} // destructor
// ----------------------------------------------------------------------
+// Set scales used to nondimensionalize physical properties.
+void
+pylith::materials::Material::normalizer(const spatialdata::units::Nondimensional& dim)
+{ // normalizer
+ if (0 == _normalizer)
+ _normalizer = new spatialdata::units::Nondimensional(dim);
+ else
+ *_normalizer = dim;
+} // normalizer
+
+// ----------------------------------------------------------------------
// Get physical property parameters and initial state (if used) from database.
void
pylith::materials::Material::initialize(
const ALE::Obj<Mesh>& mesh,
const spatialdata::geocoords::CoordSys* cs,
- pylith::feassemble::Quadrature* quadrature,
- const spatialdata::units::Nondimensional& normalizer)
+ pylith::feassemble::Quadrature* quadrature)
{ // initialize
assert(0 != _db);
assert(0 != cs);
@@ -95,7 +109,8 @@
const int numQuadPts = quadrature->numQuadPts();
const int spaceDim = quadrature->spaceDim();
-
+ double_array quadPtsGlobal(numQuadPts*spaceDim);
+
// Fiber dimension is number of quadrature points times number of
// values per parameter
const int totalPropsQuadPt = _totalPropsQuadPt;
@@ -103,15 +118,9 @@
_properties->setFiberDimension(cells, fiberDim);
mesh->allocate(_properties);
- // Fiber dimension for initial stresses is number of quadrature points times
- // initial state size.
- const int initialStateFiberDim = _initialStateSize * numQuadPts;
-
- // Container for data returned in query of initial state database
const int initialStateSize = _initialStateSize;
+ const int initialStateFiberDim = initialStateSize * numQuadPts;
double_array initialStateQueryData(initialStateSize);
-
- // Container of initial state values at cell's quadrature points
double_array initialStateCellData(initialStateFiberDim);
// If initial state is being used, create a section to hold it.
@@ -120,9 +129,9 @@
else {
_initialState = new real_section_type(mesh->comm(), mesh->debug());
assert(!_initialState.isNull());
- _initialState->setChart(real_section_type::chart_type(*std::min_element(cells->begin(), cells->end()),
- *std::max_element(cells->begin(), cells->end())+1));
-
+ _initialState->setChart(real_section_type::chart_type(
+ *std::min_element(cells->begin(), cells->end()),
+ *std::max_element(cells->begin(), cells->end())+1));
_initialState->setFiberDimension(cells, initialStateFiberDim);
mesh->allocate(_initialState);
@@ -135,12 +144,12 @@
const int numValues = _numDBValues;
_db->open();
_db->queryVals(_dbValues, numValues);
+
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->lengthScale();
- // Container for data returned in query of database
- double_array queryData(numValues);
-
- // Container of parameters at cell's quadrature points.
- double_array cellData(fiberDim);
+ double_array queryData(numValues); // data returned in query
+ double_array cellData(fiberDim); // Parameters at cell's quad pts
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
@@ -148,52 +157,48 @@
// Compute geometry information for current cell
quadrature->computeGeometry(mesh, coordinates, *c_iter);
- const double_array& quadPts = quadrature->quadPts();
+ const double_array& quadPtsNonDim = quadrature->quadPts();
+ quadPtsGlobal = quadPtsNonDim;
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
+ lengthScale);
// Loop over quadrature points in cell and query database
for (int iQuadPt=0, index=0;
iQuadPt < numQuadPts;
++iQuadPt, index+=spaceDim) {
- int err = _db->query(&queryData[0], numValues, &quadPts[index],
+ int err = _db->query(&queryData[0], numValues, &quadPtsGlobal[index],
spaceDim, cs);
-
if (err) {
std::ostringstream msg;
msg << "Could not find parameters for physical properties at \n"
<< "(";
for (int i=0; i < spaceDim; ++i)
- msg << " " << quadPts[index+i];
+ msg << " " << quadPtsGlobal[index+i];
msg << ") in material " << _label << "\n"
<< "using spatial database '" << _db->label() << "'.";
throw std::runtime_error(msg.str());
} // if
_dbToProperties(&cellData[totalPropsQuadPt*iQuadPt], queryData);
-#if 0
_nondimProperties(&cellData[totalPropsQuadPt*iQuadPt], totalPropsQuadPt);
-#endif
if (0 != _initialStateDB) {
err = _initialStateDB->query(&initialStateQueryData[0],
initialStateSize,
- &quadPts[index], spaceDim, cs);
-
+ &quadPtsGlobal[index], spaceDim, cs);
if (err) {
std::ostringstream msg;
msg << "Could not find initial state values at \n" << "(";
for (int i=0; i < spaceDim; ++i)
- msg << " " << quadPts[index+i];
+ msg << " " << quadPtsGlobal[index+i];
msg << ") in material " << _label << "\n"
<< "using spatial database '" << _initialStateDB->label() << "'.";
throw std::runtime_error(msg.str());
} // if
-#if 0 // nondimensionalize initial state
- _nondimInitialState(&initialStateCellData[initialStateSize*iQuadPt],
- &initialStateQueryData[0], initialStateSize);
-#else
- memcpy(&initialStateCellData[initialStateSize * iQuadPt],
+ // nondimensionalize initial state
+ _nondimInitState(&initialStateQueryData[0], initialStateSize);
+ memcpy(&initialStateCellData[iQuadPt*initialStateSize],
&initialStateQueryData[0],
- initialStateSize * sizeof(double));
-#endif
+ initialStateSize*sizeof(double));
} // if
} // for
@@ -275,29 +280,30 @@
if (field->isNull() ||
totalFiberDim != (*field)->getFiberDimension(*cells->begin())) {
*field = new real_section_type(mesh->comm(), mesh->debug());
- (*field)->setChart(real_section_type::chart_type(*std::min_element(cells->begin(), cells->end()),
- *std::max_element(cells->begin(), cells->end())+1));
+ (*field)->setChart(real_section_type::chart_type(
+ *std::min_element(cells->begin(), cells->end()),
+ *std::max_element(cells->begin(), cells->end())+1));
(*field)->setFiberDimension(cells, totalFiberDim);
mesh->allocate(*field);
} // if
// Buffer for property at cell's quadrature points
- double_array fieldCell(fiberDim * numQuadPts);
+ const int totalPropsQuadPt = _totalPropsQuadPt;
+ double_array fieldCell(fiberDim*numQuadPts);
+ double_array propertiesCell(totalPropsQuadPt*numQuadPts);
// Loop over cells
- const int totalPropsQuadPt = _totalPropsQuadPt;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
++c_iter) {
- const real_section_type::value_type* propVals =
- _properties->restrictPoint(*c_iter);
+ _properties->restrictPoint(*c_iter,
+ &propertiesCell[0], propertiesCell.size());
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
-#if 0
- _dimProperties(&propVals[totalPropsQuadPt*iQuad], totalPropsQuadPt);
-#endif
+ _dimProperties(&propertiesCell[iQuad*totalPropsQuadPt],
+ totalPropsQuadPt);
memcpy(&fieldCell[iQuad*fiberDim],
- &propVals[iQuad*totalPropsQuadPt+propOffset],
+ &propertiesCell[iQuad*totalPropsQuadPt+propOffset],
fiberDim*sizeof(double));
} // for
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.hh 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.hh 2008-12-02 20:06:35 UTC (rev 13439)
@@ -145,18 +145,22 @@
*/
double timeStep(void) const;
+ /** Set scales used to nondimensionalize physical properties.
+ *
+ * @param dim Nondimensionalizer
+ */
+ void normalizer(const spatialdata::units::Nondimensional& dim);
+
/** Initialize material by getting physical property parameters from
* database.
*
* @param mesh PETSc mesh
* @param cs Coordinate system associated with mesh
* @param quadrature Quadrature for finite-element integration
- * @param normalizer Nondimensionalization of scales.
*/
void initialize(const ALE::Obj<Mesh>& mesh,
const spatialdata::geocoords::CoordSys* cs,
- pylith::feassemble::Quadrature* quadrature,
- const spatialdata::units::Nondimensional& normalizer);
+ pylith::feassemble::Quadrature* quadrature);
/** Get flag indicating whether Jacobian matrix must be reformed for
* current state.
@@ -202,6 +206,42 @@
void _dbToProperties(double* const propValues,
const double_array& dbValues) const = 0;
+ /** Nondimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ virtual
+ void _nondimProperties(double* const values,
+ const int nvalues) const = 0;
+
+ /** Dimensionalize properties.
+ *
+ * @param values Array of property values.
+ * @param nvalues Number of values.
+ */
+ virtual
+ void _dimProperties(double* const values,
+ const int nvalues) const = 0;
+
+ /** Nondimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ virtual
+ void _nondimInitState(double* const values,
+ const int nvalues) const = 0;
+
+ /** Dimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ virtual
+ void _dimInitState(double* const values,
+ const int nvalues) const = 0;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
@@ -221,6 +261,8 @@
/// Section containing the initial state variables for the material.
ALE::Obj<real_section_type> _initialState;
+
+ spatialdata::units::Nondimensional* _normalizer; ///< Nondimensionalizer
int _totalPropsQuadPt; ///< Total # of property values per quad point.
int _dimension; ///< Spatial dimension associated with material.
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2008-12-02 20:06:35 UTC (rev 13439)
@@ -18,10 +18,12 @@
#include "pylith/utils/array.hh" // USES double_array
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
#include "petsc.h" // USES PetscLogFlops
-#include <assert.h> // USES assert()
-#include <string.h> // USES memcpy()
+#include <cassert> // USES assert()
+#include <cstring> // USES memcpy()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
@@ -116,9 +118,10 @@
const double vp = dbValues[_MaxwellIsotropic3D::didVp];
const double viscosity = dbValues[_MaxwellIsotropic3D::didViscosity];
- if (density < 0.0 || vs < 0.0 || vp < 0.0 || viscosity < 0.0) {
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || viscosity <= 0.0) {
std::ostringstream msg;
- msg << "Spatial database returned negative value for physical properties.\n"
+ msg << "Spatial database returned nonpositive value for physical "
+ << "properties.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n"
@@ -129,9 +132,9 @@
const double mu = density * vs*vs;
const double lambda = density * vp*vp - 2.0*mu;
- if (lambda < 0.0) {
+ if (lambda <= 0.0) {
std::ostringstream msg;
- msg << "Attempted to set Lame's constant lambda to negative value.\n"
+ msg << "Attempted to set Lame's constant lambda to nonpositive value.\n"
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n";
@@ -151,6 +154,96 @@
} // _dbToProperties
// ----------------------------------------------------------------------
+// Nondimensionalize properties.
+void
+pylith::materials::MaxwellIsotropic3D::_nondimProperties(double* const values,
+ const int nvalues) const
+{ // _nondimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _totalPropsQuadPt);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ const double timeScale = _normalizer->timeScale();
+ values[_MaxwellIsotropic3D::pidDensity] =
+ _normalizer->nondimensionalize(values[_MaxwellIsotropic3D::pidDensity],
+ densityScale);
+ values[_MaxwellIsotropic3D::pidMu] =
+ _normalizer->nondimensionalize(values[_MaxwellIsotropic3D::pidMu],
+ pressureScale);
+ values[_MaxwellIsotropic3D::pidLambda] =
+ _normalizer->nondimensionalize(values[_MaxwellIsotropic3D::pidLambda],
+ pressureScale);
+ values[_MaxwellIsotropic3D::pidMaxwellTime] =
+ _normalizer->nondimensionalize(values[_MaxwellIsotropic3D::pidMaxwellTime],
+ timeScale);
+
+ PetscLogFlops(4);
+} // _nondimProperties
+
+// ----------------------------------------------------------------------
+// Dimensionalize properties.
+void
+pylith::materials::MaxwellIsotropic3D::_dimProperties(double* const values,
+ const int nvalues) const
+{ // _dimProperties
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _totalPropsQuadPt);
+
+ const double densityScale = _normalizer->densityScale();
+ const double pressureScale = _normalizer->pressureScale();
+ const double timeScale = _normalizer->timeScale();
+ values[_MaxwellIsotropic3D::pidDensity] =
+ _normalizer->dimensionalize(values[_MaxwellIsotropic3D::pidDensity],
+ densityScale);
+ values[_MaxwellIsotropic3D::pidMu] =
+ _normalizer->dimensionalize(values[_MaxwellIsotropic3D::pidMu],
+ pressureScale);
+ values[_MaxwellIsotropic3D::pidLambda] =
+ _normalizer->dimensionalize(values[_MaxwellIsotropic3D::pidLambda],
+ pressureScale);
+ values[_MaxwellIsotropic3D::pidMaxwellTime] =
+ _normalizer->dimensionalize(values[_MaxwellIsotropic3D::pidMaxwellTime],
+ timeScale);
+
+ PetscLogFlops(4);
+} // _dimProperties
+
+// ----------------------------------------------------------------------
+// Nondimensionalize initial state.
+void
+pylith::materials::MaxwellIsotropic3D::_nondimInitState(double* const values,
+ const int nvalues) const
+{ // _nondimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _MaxwellIsotropic3D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->nondimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _nondimInitState
+
+// ----------------------------------------------------------------------
+// Dimensionalize initial state.
+void
+pylith::materials::MaxwellIsotropic3D::_dimInitState(double* const values,
+ const int nvalues) const
+{ // _dimInitState
+ assert(0 != _normalizer);
+ assert(0 != values);
+ assert(nvalues == _MaxwellIsotropic3D::numInitialStateDBValues);
+
+ const double pressureScale = _normalizer->pressureScale();
+ _normalizer->dimensionalize(values, nvalues, pressureScale);
+
+ PetscLogFlops(nvalues);
+} // _dimInitState
+
+// ----------------------------------------------------------------------
// Compute density at location from properties.
void
pylith::materials::MaxwellIsotropic3D::_calcDensity(double* const density,
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2008-12-02 16:09:37 UTC (rev 13438)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.hh 2008-12-02 20:06:35 UTC (rev 13439)
@@ -82,6 +82,38 @@
void _dbToProperties(double* const propValues,
const 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;
+
+ /** Nondimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _nondimInitState(double* const values,
+ const int nvalues) const;
+
+ /** Dimensionalize initial state.
+ *
+ * @param values Array of initial state values.
+ * @param nvalues Number of values.
+ */
+ void _dimInitState(double* const values,
+ const int nvalues) const;
+
/** Compute density from properties.
*
* @param density Array for density.
More information about the CIG-COMMITS
mailing list