[cig-commits] r13421 - in short/3D/PyLith/trunk: libsrc/bc libsrc/faults libsrc/feassemble libsrc/materials modulesrc/materials pylith/bc pylith/faults pylith/feassemble pylith/materials pylith/problems unittests/libtests/faults unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Sat Nov 29 15:58:23 PST 2008
Author: brad
Date: 2008-11-29 15:58:23 -0800 (Sat, 29 Nov 2008)
New Revision: 13421
Modified:
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.cc
short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.hh
short/3D/PyLith/trunk/libsrc/bc/DirichletPoints.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc
short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh
short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.cc
short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.hh
short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc
short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc
short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh
short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh
short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.cc
short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.hh
short/3D/PyLith/trunk/libsrc/feassemble/Constraint.cc
short/3D/PyLith/trunk/libsrc/feassemble/Constraint.hh
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/Material.cc
short/3D/PyLith/trunk/libsrc/materials/Material.hh
short/3D/PyLith/trunk/modulesrc/materials/materials.pyxe.src
short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py
short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py
short/3D/PyLith/trunk/pylith/bc/DirichletPoints.py
short/3D/PyLith/trunk/pylith/bc/Neumann.py
short/3D/PyLith/trunk/pylith/faults/Fault.py
short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
short/3D/PyLith/trunk/pylith/feassemble/Constraint.py
short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py
short/3D/PyLith/trunk/pylith/materials/Material.py
short/3D/PyLith/trunk/pylith/problems/Explicit.py
short/3D/PyLith/trunk/pylith/problems/Formulation.py
short/3D/PyLith/trunk/pylith/problems/Implicit.py
short/3D/PyLith/trunk/pylith/problems/Problem.py
short/3D/PyLith/trunk/pylith/problems/TimeDependent.py
short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
Log:
Worked on nondimensionalization.
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -20,6 +20,7 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include <Selection.hh> // USES submesh algorithms
@@ -162,6 +163,11 @@
mesh->getRealSection("coordinates");
assert(!coordinates.isNull());
+ assert(0 != _normalizer);
+ const double densityScale = _normalizer->densityScale();
+ const double velocityScale =
+ _normalizer->lengthScale() / _normalizer->timeScale();
+
for(Mesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
++c_iter) {
@@ -184,9 +190,17 @@
<< "using spatial database " << _db->label() << ".";
throw std::runtime_error(msg.str());
} // if
- const double constTangential =
- (3 == numValues) ? queryData[0]*queryData[2] : 0.0;
- const double constNormal = queryData[0]*queryData[1];
+ // Nondimensionalize damping constants
+ const double densityN =
+ _normalizer->nondimensionalize(queryData[0], densityScale);
+ const double vpN =
+ _normalizer->nondimensionalize(queryData[1], velocityScale);
+ const double vsN = (3 == numValues) ?
+ _normalizer->nondimensionalize(queryData[2], velocityScale) :
+ 0.0;
+
+ const double constTangential = densityN * vsN;
+ const double constNormal = densityN * vpN;
const int numTangential = spaceDim-1;
for (int iDim=0; iDim < numTangential; ++iDim)
dampingConstsLocal[iDim] = constTangential;
Modified: short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -16,6 +16,7 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include <Selection.hh> // USES submesh algorithms
@@ -117,6 +118,11 @@
double_array queryValues(2*numFixedDOF);
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->lengthScale();
+ const double velocityScale =
+ _normalizer->lengthScale() / _normalizer->timeScale();
+
for (Mesh::label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter) {
@@ -133,6 +139,8 @@
msg << ") using spatial database " << _db->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ for (int i=0; i < numFixedDOF; ++i)
+ _normalizer->nondimensionalize(queryValues[i], lengthScale);
err = _dbRate->query(&queryValues[numFixedDOF], numFixedDOF, vCoords,
spaceDim, cs);
@@ -144,6 +152,8 @@
msg << ") using spatial database " << _dbRate->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ for (int i=0; i < numFixedDOF; ++i)
+ _normalizer->nondimensionalize(queryValues[numFixedDOF+i], velocityScale);
_values->updatePoint(*v_iter, &queryValues[0]);
} // for
@@ -310,14 +320,31 @@
const ALE::Obj<Mesh>& mesh,
topology::FieldsManager* const fields)
{ // getVertexField
-
+ assert(0 != fieldType);
+ assert(0 != name);
+ assert(!_boundaryMesh.isNull());
+ assert(!_values.isNull());
+ assert(0 != _normalizer);
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ _boundaryMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ ALE::Obj<real_section_type> field = 0;
+ int fiberDim = 0;
+ double scale = 0.0;
if (0 == strcasecmp(name, "initial")) {
*fieldType = VECTOR_FIELD;
- _buffer = _values->getFibration(0);
+ field = _values->getFibration(0);
+ fiberDim =
+ (vertices->size() > 0) ? field->getFiberDimension(*vertices->begin()) : 0;
+ scale = _normalizer->lengthScale();
} else if (0 == strcasecmp(name, "rate-of-change")) {
*fieldType = VECTOR_FIELD;
- _buffer = _values->getFibration(1);
+ field = _values->getFibration(0);
+ fiberDim =
+ (vertices->size() > 0) ? field->getFiberDimension(*vertices->begin()) : 0;
+ scale = _normalizer->lengthScale() / _normalizer->timeScale();
} else {
std::ostringstream msg;
msg << "Unknown field '" << name << "' requested for Dirichlet BC '"
@@ -325,6 +352,33 @@
throw std::runtime_error(msg.str());
} // else
+ // Allocate buffer if necessary
+ if (_buffer.isNull()) {
+ _buffer = new real_section_type(_boundaryMesh->comm(),
+ _boundaryMesh->debug());
+ _buffer->setChart(real_section_type::chart_type(
+ *std::min_element(vertices->begin(),
+ vertices->end()),
+ *std::max_element(vertices->begin(),
+ vertices->end())+1));
+ _buffer->setFiberDimension(vertices, fiberDim);
+ _boundaryMesh->allocate(_buffer);
+ } // if
+
+ // dimensionalize values
+ double_array pointValues(fiberDim);
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter) {
+ assert(fiberDim == field->getFiberDimension(*v_iter));
+ assert(fiberDim == _buffer->getFiberDimension(*v_iter));
+ const real_section_type::value_type* values =
+ field->restrictPoint(*v_iter);
+ for (int i=0; i < fiberDim; ++i)
+ pointValues[i] = _normalizer->dimensionalize(values[i], scale);
+ _buffer->updatePointAll(*v_iter, &pointValues[0]);
+ } // for
+
return _buffer;
} // getVertexField
Modified: short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.hh 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.hh 2008-11-29 23:58:23 UTC (rev 13421)
@@ -148,7 +148,7 @@
/// Initial values and rate of change of values at DOF.
ALE::Obj<real_section_type> _values;
- ALE::Obj<real_section_type> _buffer; ///< Buffer for values.
+ ALE::Obj<real_section_type> _buffer; ///< Buffer for output.
ALE::Obj<Mesh> _boundaryMesh; ///< Boundary mesh.
int_array _fixedDOF; ///< Indices of fixed degrees of freedom
Modified: short/3D/PyLith/trunk/libsrc/bc/DirichletPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/DirichletPoints.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/bc/DirichletPoints.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -16,6 +16,7 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include <string.h> // USES strcpy()
#include <assert.h> // USES assert()
@@ -95,6 +96,11 @@
assert(!coordinates.isNull());
const int spaceDim = cs->spaceDim();
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->lengthScale();
+ const double velocityScale =
+ _normalizer->lengthScale() / _normalizer->timeScale();
+
_valuesInitial.resize(numPoints*numFixedDOF);
_valuesRate.resize(numPoints*numFixedDOF);
double_array queryValues(numFixedDOF);
@@ -113,7 +119,8 @@
throw std::runtime_error(msg.str());
} // if
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- _valuesInitial[numFixedDOF*iPoint+iDOF] = queryValues[iDOF];
+ _valuesInitial[numFixedDOF*iPoint+iDOF] =
+ _normalizer->nondimensionalize(queryValues[iDOF], lengthScale);
err = _dbRate->query(&queryValues[0], numFixedDOF, vCoords,
spaceDim, cs);
@@ -126,7 +133,8 @@
throw std::runtime_error(msg.str());
} // if
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- _valuesRate[numFixedDOF*iPoint+iDOF] = queryValues[iDOF];
+ _valuesRate[numFixedDOF*iPoint+iDOF] =
+ _normalizer->nondimensionalize(queryValues[iDOF], velocityScale);
} // for
_db->close();
_dbRate->close();
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -20,6 +20,7 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include <Selection.hh> // USES submesh algorithms
@@ -100,8 +101,6 @@
const Mesh::label_sequence::iterator cellsBegin = cells->begin();
const Mesh::label_sequence::iterator cellsEnd = cells->end();
- // std::cout << "cellsBegin: " << *cellsBegin << std::endl;
- // std::cout << "cellsEnd: " << *cellsEnd << std::endl;
const int boundaryDepth = _boundaryMesh->depth()-1; //depth of boundary cells
// Make sure surface cells are compatible with quadrature.
@@ -185,6 +184,9 @@
assert(!coordinates.isNull());
// coordinates->view("Mesh coordinates from Neumann::initialize");
+ assert(0 != _normalizer);
+ const double pressureScale = _normalizer->pressureScale();
+
// Loop over cells in boundary mesh, compute orientations, and then
// compute corresponding traction vector in global coordinates
// (store values in _tractionGlobal).
@@ -196,18 +198,9 @@
const double_array& quadPts = _quadrature->quadPts();
_boundaryMesh->restrictClosure(coordinates, *c_iter,
&cellVertices[0], cellVertices.size());
- /* Debugging stuff
- std::cout << "cellVertices: " << std::endl;
- for(int iTest = 0; iTest < numBasis; ++iTest) {
- for(int iDim = 0; iDim < spaceDim; ++iDim) {
- std::cout << " " << cellVertices[iDim+spaceDim*iTest];
- } // for
- std::cout << std::endl;
- } // for
- */
cellTractionsGlobal = 0.0;
- for(int iQuad = 0, iRef=0, iSpace=0; iQuad < numQuadPts;
+ for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts;
++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
// Get traction vector in local coordinate system at quadrature point
const int err = _db->query(&tractionDataLocal[0], spaceDim,
@@ -222,6 +215,8 @@
<< "using spatial database " << _db->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ _normalizer->nondimensionalize(&tractionDataLocal[0], spaceDim,
+ pressureScale);
// Compute Jacobian and determinant at quadrature point, then get
// orientation.
@@ -362,11 +357,27 @@
const ALE::Obj<Mesh>& mesh,
topology::FieldsManager* const fields)
{ // cellField
- _tractions->view("TRACTIONS");
+ assert(0 != fieldType);
+ assert(0 != name);
+ assert(!_boundaryMesh.isNull());
+ assert(!_tractions.isNull());
+ assert(0 != _normalizer);
+ const ALE::Obj<Mesh::label_sequence>& cells = _boundaryMesh->heightStratum(1);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int spaceDim = _quadrature->spaceDim();
+
+ ALE::Obj<real_section_type> field = 0;
+ int fiberDim = 0;
+ double scale = 0.0;
if (0 == strcasecmp(name, "tractions")) {
*fieldType = OTHER_FIELD;
- return _tractions;
+ field = _tractions;
+ fiberDim = spaceDim * numQuadPts;
+ scale = _normalizer->pressureScale();
} else {
std::ostringstream msg;
msg << "Unknown field '" << name << "' requested for Neumann BC '"
@@ -374,7 +385,32 @@
throw std::runtime_error(msg.str());
} // else
- return _tractions;
+ // Allocate buffer if necessary
+ if (_buffer.isNull()) {
+ _buffer = new real_section_type(_boundaryMesh->comm(), _boundaryMesh->debug());
+ assert(!_buffer.isNull());
+ _buffer->setChart(real_section_type::chart_type(
+ *std::min_element(cells->begin(), cells->end()),
+ *std::max_element(cells->begin(), cells->end())+1));
+ _buffer->setFiberDimension(cells, fiberDim);
+ _boundaryMesh->allocate(_buffer);
+ } // if
+
+ // dimensionalize values
+ double_array cellValues(fiberDim);
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ assert(fiberDim == field->getFiberDimension(*c_iter));
+ assert(fiberDim == _buffer->getFiberDimension(*c_iter));
+ const real_section_type::value_type* values =
+ field->restrictPoint(*c_iter);
+ for (int i=0; i < fiberDim; ++i)
+ cellValues[i] = _normalizer->dimensionalize(values[i], scale);
+ _buffer->updatePointAll(*c_iter, &cellValues[0]);
+ } // for
+
+ return _buffer;
} // cellField
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2008-11-29 23:58:23 UTC (rev 13421)
@@ -141,6 +141,8 @@
/// Traction vector in global coordinates at integration points.
ALE::Obj<real_section_type> _tractions;
+ ALE::Obj<real_section_type> _buffer; ///< Buffer for output.
+
}; // class Neumann
// #include "Neumann.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -19,6 +19,7 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include <assert.h> // USES assert()
#include <sstream> // USES std::ostringstream
@@ -56,9 +57,10 @@
// Initialize slip time function.
void
pylith::faults::BruneSlipFn::initialize(
- const ALE::Obj<Mesh>& faultMesh,
- const spatialdata::geocoords::CoordSys* cs,
- const double originTime)
+ const ALE::Obj<Mesh>& faultMesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ const spatialdata::units::Nondimensional& normalizer,
+ const double originTime)
{ // initialize
assert(!faultMesh.isNull());
assert(0 != cs);
@@ -127,8 +129,12 @@
faultMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
+ const double lengthScale = normalizer.lengthScale();
+ const double timeScale = normalizer.timeScale();
+ const double velocityScale =
+ normalizer.lengthScale() / normalizer.timeScale();
+
double_array paramsVertex(fiberDim);
-
for (Mesh::label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter) {
@@ -148,6 +154,8 @@
msg << ") using spatial database " << _dbFinalSlip->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ normalizer.nondimensionalize(¶msVertex[indexFinalSlip], spaceDim,
+ lengthScale);
err = _dbPeakRate->query(¶msVertex[indexPeakRate], 1,
coordsVertex, spaceDim, cs);
if (err) {
@@ -158,6 +166,8 @@
msg << ") using spatial database " << _dbPeakRate->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ normalizer.nondimensionalize(¶msVertex[indexPeakRate], 1,
+ velocityScale);
err = _dbSlipTime->query(¶msVertex[indexSlipTime], 1,
coordsVertex, spaceDim, cs);
@@ -169,6 +179,9 @@
msg << ") using spatial database " << _dbSlipTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ normalizer.nondimensionalize(¶msVertex[indexSlipTime], 1,
+ timeScale);
+
// add origin time to rupture time
paramsVertex[indexSlipTime] += originTime;
Modified: short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh 2008-11-29 23:58:23 UTC (rev 13421)
@@ -39,6 +39,9 @@
namespace spatialdb {
class SpatialDB;
} // spatialdb
+ namespace units {
+ class Nondimensional;
+ } // units
} // spatialdata
/// C++ implementation of Brune slip time function.
@@ -78,10 +81,12 @@
*
* @param faultMesh Finite-element mesh of fault.
* @param cs Coordinate system for mesh.
+ * @param normalizer Nondimensionalization of scales.
* @param originTime Origin time for earthquake source.
*/
void initialize(const ALE::Obj<Mesh>& faultMesh,
const spatialdata::geocoords::CoordSys* cs,
+ const spatialdata::units::Nondimensional& normalizer,
const double originTime =0.0);
/** Get slip on fault surface at time t.
Modified: short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -19,6 +19,7 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include <assert.h> // USES assert()
#include <sstream> // USES std::ostringstream
@@ -53,9 +54,10 @@
// Initialize slip time function.
void
pylith::faults::ConstRateSlipFn::initialize(
- const ALE::Obj<Mesh>& faultMesh,
- const spatialdata::geocoords::CoordSys* cs,
- const double originTime)
+ const ALE::Obj<Mesh>& faultMesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ const spatialdata::units::Nondimensional& normalizer,
+ const double originTime)
{ // initialize
assert(!faultMesh.isNull());
assert(0 != cs);
@@ -116,8 +118,11 @@
faultMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
+ const double timeScale = normalizer.timeScale();
+ const double velocityScale =
+ normalizer.lengthScale() / normalizer.timeScale();
+
double_array paramsVertex(fiberDim);
-
for (Mesh::label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter) {
@@ -137,6 +142,8 @@
msg << ") using spatial database " << _dbSlipRate->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ normalizer.nondimensionalize(¶msVertex[indexSlipRate], spaceDim,
+ velocityScale);
err = _dbSlipTime->query(¶msVertex[indexSlipTime], 1,
coordsVertex, spaceDim, cs);
@@ -148,6 +155,8 @@
msg << ") using spatial database " << _dbSlipTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ normalizer.nondimensionalize(¶msVertex[indexSlipTime], spaceDim,
+ timeScale);
// add origin time to rupture time
paramsVertex[indexSlipTime] += originTime;
Modified: short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.hh 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.hh 2008-11-29 23:58:23 UTC (rev 13421)
@@ -38,6 +38,9 @@
namespace spatialdb {
class SpatialDB;
} // spatialdb
+ namespace units {
+ class Nondimensional;
+ } // units
} // spatialdata
/// C++ implementation of ConstRate slip time function.
@@ -71,10 +74,12 @@
*
* @param faultMesh Finite-element mesh of fault.
* @param cs Coordinate system for mesh.
+ * @param normalizer Nondimensionalization of scales.
* @param originTime Origin time for earthquake source.
*/
void initialize(const ALE::Obj<Mesh>& faultMesh,
const spatialdata::geocoords::CoordSys* cs,
+ const spatialdata::units::Nondimensional& normalizer,
const double originTime =0.0);
/** Get slip on fault surface at time t.
Modified: short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -16,6 +16,8 @@
#include "SlipTimeFn.hh" // USES SlipTimeFn
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
#include <assert.h> // USES assert()
// ----------------------------------------------------------------------
@@ -61,11 +63,13 @@
// Initialize slip time function.
void
pylith::faults::EqKinSrc::initialize(
- const ALE::Obj<Mesh>& faultMesh,
- const spatialdata::geocoords::CoordSys* cs)
+ const ALE::Obj<Mesh>& faultMesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ const spatialdata::units::Nondimensional& normalizer)
{ // initialize
+ normalizer.nondimensionalize(&_originTime, 1, normalizer.timeScale());
assert(0 != _slipfn);
- _slipfn->initialize(faultMesh, cs, _originTime);
+ _slipfn->initialize(faultMesh, cs, normalizer, _originTime);
} // initialize
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh 2008-11-29 23:58:23 UTC (rev 13421)
@@ -39,10 +39,12 @@
namespace spatialdb {
class SpatialDB;
} // spatialdb
-
namespace geocoords {
class CoordSys;
- } // geocoords
+ } // geocoord
+ namespace units {
+ class Nondimensional;
+ } // units
} // spatialdata
/// C++ oject for managing parameters for a kinematic earthquake source.
@@ -81,9 +83,11 @@
*
* @param faultMesh Finite-element mesh of fault.
* @param cs Coordinate system for mesh
+ * @param normalizer Nondimensionalization of scales.
*/
void initialize(const ALE::Obj<Mesh>& faultMesh,
- const spatialdata::geocoords::CoordSys* cs);
+ const spatialdata::geocoords::CoordSys* cs,
+ const spatialdata::units::Nondimensional& normalizer);
/** Get slip on fault surface at time t.
*
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -25,6 +25,7 @@
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/spatialdb/SpatialDB.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include <Completion.hh> // USES completeSection
#include <Selection.hh> // Algorithms for submeshes
@@ -86,7 +87,6 @@
CohesiveTopology::createParallel(&_faultMesh, &_cohesiveToFault, mesh, id(),
_useLagrangeConstraints());
//_faultMesh->getLabel("height")->view("Fault mesh height");
-
//_faultMesh->view("FAULT MESH");
// Setup pseudo-stiffness of cohesive cells to improve conditioning
@@ -105,25 +105,24 @@
++s_iter) {
EqKinSrc* src = s_iter->second;
assert(0 != src);
- src->initialize(_faultMesh, cs);
+ src->initialize(_faultMesh, cs, *_normalizer);
} // for
// Allocate slip field
const ALE::Obj<Mesh::label_sequence>& vertices = _faultMesh->depthStratum(0);
_slip = new real_section_type(_faultMesh->comm(), _faultMesh->debug());
- _slip->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(),
- vertices->end()),
- *std::max_element(vertices->begin(), vertices->end())+1));
+ _slip->setChart(real_section_type::chart_type(
+ *std::min_element(vertices->begin(), vertices->end()),
+ *std::max_element(vertices->begin(), vertices->end())+1));
_slip->setFiberDimension(vertices, cs->spaceDim());
_faultMesh->allocate(_slip);
assert(!_slip.isNull());
// Allocate cumulative slip field
_cumSlip = new real_section_type(_faultMesh->comm(), _faultMesh->debug());
- _cumSlip->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(),
- vertices->end()),
- *std::max_element(vertices->begin(),
- vertices->end())+1));
+ _cumSlip->setChart(real_section_type::chart_type(
+ *std::min_element(vertices->begin(), vertices->end()),
+ *std::max_element(vertices->begin(), vertices->end())+1));
_cumSlip->setFiberDimension(vertices, cs->spaceDim());
_faultMesh->allocate(_cumSlip);
assert(!_cumSlip.isNull());
@@ -840,6 +839,9 @@
_faultMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
+ assert(0 != _normalizer);
+ const double pressureScale = _normalizer->pressureScale();
+
double_array matprops(numStiffnessVals);
int count = 0;
for (Mesh::label_sequence::iterator v_iter=vertices->begin();
@@ -862,8 +864,8 @@
const double density = matprops[0];
const double vs = matprops[1];
const double mu = density * vs*vs;
- //const double mu = 1.0;
- _pseudoStiffness->updatePoint(*v_iter, &mu);
+ const double muN = _normalizer->nondimensionalize(mu, pressureScale);
+ _pseudoStiffness->updatePoint(*v_iter, &muN);
} // for
PetscLogFlops(count * 2);
@@ -1005,10 +1007,9 @@
if (tractions->isNull() ||
fiberDim != (*tractions)->getFiberDimension(*fvertices->begin())) {
*tractions = new real_section_type(_faultMesh->comm(), _faultMesh->debug());
- (*tractions)->setChart(real_section_type::chart_type(*std::min_element(fvertices->begin(),
- fvertices->end()),
- *std::max_element(fvertices->begin(),
- fvertices->end())+1));
+ (*tractions)->setChart(real_section_type::chart_type(
+ *std::min_element(fvertices->begin(), fvertices->end()),
+ *std::max_element(fvertices->begin(), fvertices->end())+1));
(*tractions)->setFiberDimension(fvertices, fiberDim);
_faultMesh->allocate(*tractions);
assert(!tractions->isNull());
@@ -1063,7 +1064,9 @@
_faultMesh->debug());
const ALE::Obj<Mesh::label_sequence>& vertices =
_faultMesh->depthStratum(0);
- _bufferVertexScalar->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(), vertices->end()), *std::max_element(vertices->begin(), vertices->end())+1));
+ _bufferVertexScalar->setChart(real_section_type::chart_type(
+ *std::min_element(vertices->begin(), vertices->end()),
+ *std::max_element(vertices->begin(), vertices->end())+1));
_bufferVertexScalar->setFiberDimension(vertices, fiberDim);
_faultMesh->allocate(_bufferVertexScalar);
} // if
@@ -1081,7 +1084,9 @@
_faultMesh->debug());
const ALE::Obj<Mesh::label_sequence>& vertices =
_faultMesh->depthStratum(0);
- _bufferVertexVector->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(), vertices->end()), *std::max_element(vertices->begin(), vertices->end())+1));
+ _bufferVertexVector->setChart(real_section_type::chart_type(
+ *std::min_element(vertices->begin(), vertices->end()),
+ *std::max_element(vertices->begin(), vertices->end())+1));
_bufferVertexVector->setFiberDimension(vertices, fiberDim);
_faultMesh->allocate(_bufferVertexVector);
} // if
Modified: short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -19,6 +19,7 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include <assert.h> // USES assert()
#include <sstream> // USES std::ostringstream
@@ -56,9 +57,10 @@
// Initialize slip time function.
void
pylith::faults::LiuCosSlipFn::initialize(
- const ALE::Obj<Mesh>& faultMesh,
- const spatialdata::geocoords::CoordSys* cs,
- const double originTime)
+ const ALE::Obj<Mesh>& faultMesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ const spatialdata::units::Nondimensional& normalizer,
+ const double originTime)
{ // initialize
assert(!faultMesh.isNull());
assert(0 != cs);
@@ -127,8 +129,10 @@
faultMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
+ const double lengthScale = normalizer.lengthScale();
+ const double timeScale = normalizer.timeScale();
+
double_array paramsVertex(fiberDim);
-
for (Mesh::label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter) {
@@ -148,6 +152,9 @@
msg << ") using spatial database " << _dbFinalSlip->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ normalizer.nondimensionalize(¶msVertex[indexFinalSlip], spaceDim,
+ lengthScale);
+
err = _dbRiseTime->query(¶msVertex[indexRiseTime], 1,
coordsVertex, spaceDim, cs);
if (err) {
@@ -158,6 +165,8 @@
msg << ") using spatial database " << _dbRiseTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ normalizer.nondimensionalize(¶msVertex[indexRiseTime], spaceDim,
+ timeScale);
err = _dbSlipTime->query(¶msVertex[indexSlipTime], 1,
coordsVertex, spaceDim, cs);
@@ -169,6 +178,8 @@
msg << ") using spatial database " << _dbSlipTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ normalizer.nondimensionalize(¶msVertex[indexSlipTime], spaceDim,
+ timeScale);
// add origin time to rupture time
paramsVertex[indexSlipTime] += originTime;
Modified: short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh 2008-11-29 23:58:23 UTC (rev 13421)
@@ -37,6 +37,9 @@
namespace spatialdb {
class SpatialDB;
} // spatialdb
+ namespace units {
+ class Nondimensional;
+ } // units
} // spatialdata
/// C++ implementation of LiuCos slip time function.
@@ -79,9 +82,11 @@
* @param faultMesh Finite-element mesh of fault.
* @param cs Coordinate system for mesh.
* @param originTime Origin time for earthquake source.
+ * @param normalizer Nondimensionalization of scales.
*/
void initialize(const ALE::Obj<Mesh>& faultMesh,
const spatialdata::geocoords::CoordSys* cs,
+ const spatialdata::units::Nondimensional& normalizer,
const double originTime =0.0);
/** Get slip on fault surface at time t.
Modified: short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh 2008-11-29 23:58:23 UTC (rev 13421)
@@ -35,6 +35,9 @@
namespace geocoords {
class CoordSys;
} // geocoords
+ namespace units {
+ class Nondimensional;
+ } // units
} // spatialdata
/// C++ abstract base class for Fault object.
@@ -56,11 +59,13 @@
*
* @param faultMesh Finite-element mesh of fault.
* @param cs Coordinate system for mesh
+ * @param normalizer Nondimensionalization of scales.
* @param originTime Origin time for earthquake source.
*/
virtual
void initialize(const ALE::Obj<Mesh>& faultMesh,
const spatialdata::geocoords::CoordSys* cs,
+ const spatialdata::units::Nondimensional& normalizer,
const double originTime =0.0) = 0;
/** Get slip on fault surface at time t.
Modified: short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -19,6 +19,7 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include <assert.h> // USES assert()
#include <sstream> // USES std::ostringstream
@@ -53,9 +54,10 @@
// Initialize slip time function.
void
pylith::faults::StepSlipFn::initialize(
- const ALE::Obj<Mesh>& faultMesh,
- const spatialdata::geocoords::CoordSys* cs,
- const double originTime)
+ const ALE::Obj<Mesh>& faultMesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ const spatialdata::units::Nondimensional& normalizer,
+ const double originTime)
{ // initialize
assert(!faultMesh.isNull());
assert(0 != cs);
@@ -116,8 +118,10 @@
faultMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
+ const double lengthScale = normalizer.lengthScale();
+ const double timeScale = normalizer.timeScale();
+
double_array paramsVertex(fiberDim);
-
for (Mesh::label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter) {
@@ -137,6 +141,8 @@
msg << ") using spatial database " << _dbFinalSlip->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ normalizer.nondimensionalize(¶msVertex[indexFinalSlip], spaceDim,
+ lengthScale);
err = _dbSlipTime->query(¶msVertex[indexSlipTime], 1,
coordsVertex, spaceDim, cs);
@@ -148,6 +154,8 @@
msg << ") using spatial database " << _dbSlipTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
+ normalizer.nondimensionalize(¶msVertex[indexSlipTime], spaceDim,
+ timeScale);
// add origin time to rupture time
paramsVertex[indexSlipTime] += originTime;
Modified: short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.hh 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.hh 2008-11-29 23:58:23 UTC (rev 13421)
@@ -37,6 +37,9 @@
namespace spatialdb {
class SpatialDB;
} // spatialdb
+ namespace units {
+ class Nondimensional;
+ } // units
} // spatialdata
/// C++ implementation of Step slip time function.
@@ -69,10 +72,12 @@
*
* @param faultMesh Finite-element mesh of fault.
* @param cs Coordinate system for mesh.
+ * @param normalizer Nondimensionalization of scales.
* @param originTime Origin time for earthquake source.
*/
void initialize(const ALE::Obj<Mesh>& faultMesh,
const spatialdata::geocoords::CoordSys* cs,
+ const spatialdata::units::Nondimensional& normalizer,
const double originTime =0.0);
/** Get slip on fault surface at time t.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Constraint.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Constraint.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Constraint.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -14,9 +14,12 @@
#include "Constraint.hh" // implementation of object methods
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::feassemble::Constraint::Constraint(void) :
+ _normalizer(new spatialdata::units::Nondimensional),
_useSolnIncr(false)
{ // constructor
} // constructor
@@ -25,7 +28,19 @@
// Destructor.
pylith::feassemble::Constraint::~Constraint(void)
{ // destructor
+ delete _normalizer; _normalizer = 0;
} // destructor
+// ----------------------------------------------------------------------
+// Set manager of scales used to nondimensionalize problem.
+void
+pylith::feassemble::Constraint::normalizer(const spatialdata::units::Nondimensional& dim)
+{ // normalizer
+ if (0 == _normalizer)
+ _normalizer = new spatialdata::units::Nondimensional(dim);
+ else
+ *_normalizer = dim;
+} // normalizer
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Constraint.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Constraint.hh 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Constraint.hh 2008-11-29 23:58:23 UTC (rev 13421)
@@ -29,6 +29,11 @@
} // feassemble
} // pylith
+namespace spatialdata {
+ namespace units {
+ class Nondimensional; // USES Nondimensional
+ } // units
+} // spatialdata
/// C++ abstract base class defining interface for constraints applied
/// to finite-elements.
@@ -46,6 +51,12 @@
virtual
~Constraint(void);
+ /** Set manager of scales used to nondimensionalize problem.
+ *
+ * @param dim Nondimensionalizer.
+ */
+ void normalizer(const spatialdata::units::Nondimensional& dim);
+
/** Set number of degrees of freedom that are constrained at points in field.
*
* @param field Solution field
@@ -86,6 +97,8 @@
// PROTECTED MEMBERS //////////////////////////////////////////////////
protected :
+ spatialdata::units::Nondimensional* _normalizer; ///< Nondimensionalizer.
+
/// Flag indicating whether to set constraints for a total field
/// solution or an incremental field solution
bool _useSolnIncr;
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -25,6 +25,7 @@
#include "petscmat.h" // USES PetscMat
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimendional
#include "spatialdata/spatialdb/GravityField.hh" // USES GravityField
#include <assert.h> // USES assert()
@@ -240,6 +241,11 @@
// Get density at quadrature points for this cell
const double_array& density = _material->calcDensity();
+ assert(0 != _normalizer);
+ const double gravityScale =
+ _normalizer->pressureScale() / (_normalizer->lengthScale() *
+ _normalizer->densityScale());
+
// Compute action for element body forces
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
double gravVec[spaceDim];
@@ -250,6 +256,7 @@
coords, spaceDim, cs);
if (err)
throw std::runtime_error("Unable to get gravity vector for point.");
+ _normalizer->nondimensionalize(gravVec, spaceDim, gravityScale);
const double wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
for (int iBasis=0, iQ=iQuad*numBasis;
iBasis < numBasis; ++iBasis) {
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -17,6 +17,7 @@
#include "Quadrature.hh" // USES Quadrature
#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include "spatialdata/spatialdb/GravityField.hh" // USES GravityField
#include <assert.h> // USES assert()
@@ -26,6 +27,7 @@
pylith::feassemble::Integrator::Integrator(void) :
_dt(-1.0),
_quadrature(0),
+ _normalizer(new spatialdata::units::Nondimensional),
_gravityField(0),
_cellVector(0),
_cellMatrix(0),
@@ -39,6 +41,7 @@
pylith::feassemble::Integrator::~Integrator(void)
{ // destructor
delete _quadrature; _quadrature = 0;
+ delete _normalizer; _normalizer = 0;
delete[] _cellVector; _cellVector = 0;
delete[] _cellMatrix; _cellMatrix = 0;
} // destructor
@@ -57,6 +60,17 @@
} // quadrature
// ----------------------------------------------------------------------
+// Set manager of scales used to nondimensionalize problem.
+void
+pylith::feassemble::Integrator::normalizer(const spatialdata::units::Nondimensional& dim)
+{ // normalizer
+ if (0 == _normalizer)
+ _normalizer = new spatialdata::units::Nondimensional(dim);
+ else
+ *_normalizer = dim;
+} // normalizer
+
+// ----------------------------------------------------------------------
// Set gravity field.
void
pylith::feassemble::Integrator::gravityField(spatialdata::spatialdb::GravityField* const gravityField)
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2008-11-29 23:58:23 UTC (rev 13421)
@@ -47,6 +47,9 @@
namespace geocoords {
class CoordSys; // USES CoordSys
} // geocoords
+ namespace units {
+ class Nondimensional; // USES Nondimensional
+ } // units
} // spatialdata
class pylith::feassemble::Integrator
@@ -70,6 +73,12 @@
*/
void quadrature(const Quadrature* q);
+ /** Set manager of scales used to nondimensionalize problem.
+ *
+ * @param dim Nondimensionalizer.
+ */
+ void normalizer(const spatialdata::units::Nondimensional& dim);
+
/** Set gravity field. Gravity Field should already be initialized.
*
* @param g Gravity field.
@@ -216,6 +225,7 @@
Quadrature* _quadrature; ///< Quadrature for integrating finite-element
+ spatialdata::units::Nondimensional* _normalizer; ///< Nondimensionalizer.
spatialdata::spatialdb::GravityField* _gravityField; ///< Gravity field.
/// Vector local to cell containing result of integration action
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -17,6 +17,8 @@
#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 <assert.h> // USES assert()
@@ -92,8 +94,8 @@
// Compute properties from values in spatial database.
void
pylith::materials::ElasticIsotropic3D::_dbToProperties(
- double* const propValues,
- const double_array& dbValues) const
+ double* const propValues,
+ const double_array& dbValues) const
{ // _dbToProperties
assert(0 != propValues);
const int numDBValues = dbValues.size();
@@ -103,7 +105,7 @@
const double vs = dbValues[_ElasticIsotropic3D::didVs];
const double vp = dbValues[_ElasticIsotropic3D::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"
<< "density: " << density << "\n"
@@ -115,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";
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -68,9 +68,11 @@
// ----------------------------------------------------------------------
// 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)
+pylith::materials::Material::initialize(
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ pylith::feassemble::Quadrature* quadrature,
+ const spatialdata::units::Nondimensional& normalizer)
{ // initialize
assert(0 != _db);
assert(0 != cs);
@@ -152,8 +154,8 @@
for (int iQuadPt=0, index=0;
iQuadPt < numQuadPts;
++iQuadPt, index+=spaceDim) {
- const int err = _db->query(&queryData[0], numValues, &quadPts[index],
- spaceDim, cs);
+ int err = _db->query(&queryData[0], numValues, &quadPts[index],
+ spaceDim, cs);
if (err) {
std::ostringstream msg;
@@ -166,26 +168,32 @@
throw std::runtime_error(msg.str());
} // if
_dbToProperties(&cellData[totalPropsQuadPt*iQuadPt], queryData);
+#if 0
+ _nondimProperties(&cellData[totalPropsQuadPt*iQuadPt], totalPropsQuadPt);
+#endif
if (0 != _initialStateDB) {
- const int err2 = _initialStateDB->query(&initialStateQueryData[0],
- initialStateSize,
- &quadPts[index],
- spaceDim, cs);
+ err = _initialStateDB->query(&initialStateQueryData[0],
+ initialStateSize,
+ &quadPts[index], spaceDim, cs);
- if (err2) {
+ if (err) {
std::ostringstream msg;
- msg << "Could not find initial state values at \n"
- << "(";
+ msg << "Could not find initial state values at \n" << "(";
for (int i=0; i < spaceDim; ++i)
msg << " " << quadPts[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],
&initialStateQueryData[0],
initialStateSize * sizeof(double));
+#endif
} // if
} // for
@@ -277,16 +285,21 @@
double_array fieldCell(fiberDim * 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);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+#if 0
+ _dimProperties(&propVals[totalPropsQuadPt*iQuad], totalPropsQuadPt);
+#endif
memcpy(&fieldCell[iQuad*fiberDim],
- &propVals[iQuad*_totalPropsQuadPt+propOffset],
+ &propVals[iQuad*totalPropsQuadPt+propOffset],
fiberDim*sizeof(double));
+ } // for
(*field)->updatePoint(*c_iter, &fieldCell[0]);
} // for
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.hh 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.hh 2008-11-29 23:58:23 UTC (rev 13421)
@@ -48,6 +48,9 @@
namespace geocoords {
class CoordSys; // forward declaration
} // geocoords
+ namespace units {
+ class Nondimensional; // forward declaration
+ } // units
} // spatialdata
/// C++ abstract base class for Material object.
@@ -148,10 +151,12 @@
* @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);
+ pylith::feassemble::Quadrature* quadrature,
+ const spatialdata::units::Nondimensional& normalizer);
/** Get flag indicating whether Jacobian matrix must be reformed for
* current state.
Modified: short/3D/PyLith/trunk/modulesrc/materials/materials.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/materials.pyxe.src 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/modulesrc/materials/materials.pyxe.src 2008-11-29 23:58:23 UTC (rev 13421)
@@ -71,24 +71,27 @@
return
- def initialize(self, mesh, cs, quadrature):
+ def initialize(self, mesh, cs, quadrature, normalizer):
"""
Initialize material by getting physical property parameters from
database.
"""
# create shim for method 'initialize'
- #embed{ void Material_initialize(void* objVptr, void* meshVptr, void* csVptr, void* qVptr)
+ #embed{ void Material_initialize(void* objVptr, void* meshVptr, void* csVptr, void* qVptr, void* nVptr)
try {
assert(0 != objVptr);
assert(0 != meshVptr);
assert(0 != csVptr);
assert(0 != qVptr);
+ assert(0 != nVptr);
ALE::Obj<pylith::Mesh>* mesh = (ALE::Obj<pylith::Mesh>*) meshVptr;
spatialdata::geocoords::CoordSys* cs =
(spatialdata::geocoords::CoordSys*) csVptr;
pylith::feassemble::Quadrature* quadrature =
(pylith::feassemble::Quadrature*) qVptr;
- ((pylith::materials::Material*) objVptr)->initialize(*mesh, cs, quadrature);
+ spatialdata::units::Nondimensional* normalizer =
+ (spatialdata::units::Nondimensional*) nVptr;
+ ((pylith::materials::Material*) objVptr)->initialize(*mesh, cs, quadrature, *normalizer);
} catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError,
const_cast<char*>(err.what()));
@@ -113,8 +116,12 @@
raise TypeError, \
"Argument must be extension module type " \
"'pylith::feassemble::Quadrature'."
+ if not normalizer.name == "spatialdata_units_Nondimensional":
+ raise TypeError, \
+ "Argument must be extension module type " \
+ "'spatialdata::units::Nondimensional'."
Material_initialize(self.thisptr, ptrFromHandle(mesh), ptrFromHandle(cs),
- ptrFromHandle(quadrature))
+ ptrFromHandle(quadrature), ptrFromHandle(normalizer))
return
Modified: short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -95,17 +95,17 @@
return
- def initialize(self, totalTime, numTimeSteps):
+ def initialize(self, totalTime, numTimeSteps, normalizer):
"""
Initialize AbsorbingDampers boundary condition.
"""
logEvent = "%sinit" % self._loggingPrefix
self._logger.eventBegin(logEvent)
- Integrator.initialize(self, totalTime, numTimeSteps)
+ Integrator.initialize(self, totalTime, numTimeSteps, normalizer)
self.cppHandle.quadrature = self.quadrature.cppHandle
- BoundaryCondition.initialize(self, totalTime, numTimeSteps)
+ BoundaryCondition.initialize(self, totalTime, numTimeSteps, normalizer)
self._logger.eventEnd(logEvent)
return
Modified: short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/bc/BoundaryCondition.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -116,7 +116,7 @@
return
- def initialize(self, totalTime, numTimeSteps):
+ def initialize(self, totalTime, numTimeSteps, normalizer):
"""
Initialize boundary condition.
"""
Modified: short/3D/PyLith/trunk/pylith/bc/DirichletPoints.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/DirichletPoints.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/bc/DirichletPoints.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -121,7 +121,7 @@
return
- def initialize(self, totalTime, numTimeSteps):
+ def initialize(self, totalTime, numTimeSteps, normalizer):
"""
Initialize DirichletPoints boundary condition.
"""
@@ -133,7 +133,7 @@
self.dbRate.initialize()
self.cppHandle.dbRate = self.dbRate.cppHandle
- BoundaryCondition.initialize(self, totalTime, numTimeSteps)
+ BoundaryCondition.initialize(self, totalTime, numTimeSteps, normalizer)
self._logger.eventEnd(logEvent)
return
Modified: short/3D/PyLith/trunk/pylith/bc/Neumann.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/Neumann.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/bc/Neumann.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -108,17 +108,17 @@
return
- def initialize(self, totalTime, numTimeSteps):
+ def initialize(self, totalTime, numTimeSteps, normalizer):
"""
Initialize Neumann boundary condition.
"""
logEvent = "%sinit" % self._loggingPrefix
self._logger.eventBegin(logEvent)
- Integrator.initialize(self, totalTime, numTimeSteps)
+ Integrator.initialize(self, totalTime, numTimeSteps, normalizer)
self.cppHandle.quadrature = self.quadrature.cppHandle
- BoundaryCondition.initialize(self, totalTime, numTimeSteps)
+ BoundaryCondition.initialize(self, totalTime, numTimeSteps, normalizer)
from pylith.topology.Mesh import Mesh
self.boundaryMesh = Mesh()
Modified: short/3D/PyLith/trunk/pylith/faults/Fault.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/Fault.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/faults/Fault.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -176,7 +176,7 @@
return
- def initialize(self, totalTime, numTimeSteps):
+ def initialize(self, totalTime, numTimeSteps, normalizer):
"""
Initialize fault.
"""
Modified: short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -136,7 +136,7 @@
return
- def initialize(self, totalTime, numTimeSteps):
+ def initialize(self, totalTime, numTimeSteps, normalizer):
"""
Initialize cohesive elements.
"""
@@ -144,11 +144,11 @@
self._logger.eventBegin(logEvent)
self._info.log("Initializing fault '%s'." % self.label)
- Integrator.initialize(self, totalTime, numTimeSteps)
+ Integrator.initialize(self, totalTime, numTimeSteps, normalizer)
for eqsrc in self.eqsrcs.components():
eqsrc.initialize()
- FaultCohesive.initialize(self, totalTime, numTimeSteps)
+ FaultCohesive.initialize(self, totalTime, numTimeSteps, normalizer)
self._logger.eventEnd(logEvent)
return
Modified: short/3D/PyLith/trunk/pylith/feassemble/Constraint.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/Constraint.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/feassemble/Constraint.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -80,13 +80,16 @@
return
- def initialize(self, totalTime, numTimeSteps):
+ def initialize(self, totalTime, numTimeSteps, normalizer):
"""
Do initialization.
"""
logEvent = "%sinit" % self._loggingPrefix
self._logger.eventBegin(logEvent)
+ assert(None != self.cppHandle)
+ #self.cppHandle.normalizer = normalizer.cppHandle
+
self._logger.eventEnd(logEvent)
return
Modified: short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/Integrator.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/feassemble/Integrator.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -86,7 +86,7 @@
return
- def initialize(self, totalTime, numTimeSteps):
+ def initialize(self, totalTime, numTimeSteps, normalizer):
"""
Do initialization.
"""
@@ -94,6 +94,7 @@
self._logger.eventBegin(logEvent)
assert(None != self.cppHandle)
+ #self.cppHandle.normalizer = normalizer.cppHandle
if None != self.gravityField:
self.cppHandle.gravityField = self.gravityField.cppHandle
Modified: short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/feassemble/IntegratorElasticity.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -79,7 +79,7 @@
return
- def initialize(self, totalTime, numTimeSteps):
+ def initialize(self, totalTime, numTimeSteps, normalizer):
"""
Initialize material properties.
"""
@@ -89,9 +89,9 @@
self._info.log("Initializing integrator for material '%s'." % \
self.material.label)
- Integrator.initialize(self, totalTime, numTimeSteps)
+ Integrator.initialize(self, totalTime, numTimeSteps, normalizer)
- self.material.initialize(self.mesh, totalTime, numTimeSteps)
+ self.material.initialize(self.mesh, totalTime, numTimeSteps, normalizer)
self.output.initialize(self.quadrature)
self.output.writeInfo()
self.output.open(totalTime, numTimeSteps)
Modified: short/3D/PyLith/trunk/pylith/materials/Material.py
===================================================================
--- short/3D/PyLith/trunk/pylith/materials/Material.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/materials/Material.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -127,7 +127,7 @@
return
- def initialize(self, mesh, totalTime, numTimeSteps):
+ def initialize(self, mesh, totalTime, numTimeSteps, normalizer):
"""
Initialize material property manager.
"""
@@ -144,7 +144,8 @@
self.initialStateDB.initialize()
self.cppHandle.initialStateDB = self.initialStateDB.cppHandle
self.cppHandle.initialize(mesh.cppHandle, mesh.coordsys.cppHandle,
- self.quadrature.cppHandle)
+ self.quadrature.cppHandle,
+ normalizer.cppHandle)
self._logger.eventEnd(logEvent)
return
Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -59,14 +59,14 @@
return ElasticityExplicit()
- def initialize(self, dimension):
+ def initialize(self, dimension, normalizer):
"""
Initialize problem for explicit time integration.
"""
logEvent = "%sinit" % self._loggingPrefix
self._logger.eventBegin(logEvent)
- Formulation.initialize(self, dimension)
+ Formulation.initialize(self, dimension, normalizer)
self._info.log("Creating other fields and matrices.")
self.fields.addReal("dispTpdt")
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -160,7 +160,7 @@
return
- def initialize(self, dimension):
+ def initialize(self, dimension, normalizer):
"""
Create integrators for each element family.
"""
@@ -183,12 +183,12 @@
self._info.log("Initializing integrators.")
for integrator in self.integrators:
integrator.gravityField = self.gravityField
- integrator.initialize(totalTime, numTimeSteps)
+ integrator.initialize(totalTime, numTimeSteps, normalizer)
self._debug.log(resourceUsageString())
self._info.log("Initializing constraints.")
for constraint in self.constraints:
- constraint.initialize(totalTime, numTimeSteps)
+ constraint.initialize(totalTime, numTimeSteps, normalizer)
self._debug.log(resourceUsageString())
self._info.log("Setting up solution output.")
Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -84,14 +84,14 @@
return ElasticityImplicit()
- def initialize(self, dimension):
+ def initialize(self, dimension, normalizer):
"""
Initialize problem for implicit time integration.
"""
logEvent = "%sinit" % self._loggingPrefix
self._logger.eventBegin(logEvent)
- Formulation.initialize(self, dimension)
+ Formulation.initialize(self, dimension, normalizer)
self._info.log("Creating other fields.")
self._debug.log(resourceUsageString())
Modified: short/3D/PyLith/trunk/pylith/problems/Problem.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Problem.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/problems/Problem.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -70,6 +70,7 @@
## @li \b useGravity Gravity on (true) or off (false).
##
## \b Facilities
+ ## @li \b normalizer Nondimensionalizer for problem.
## @li \b materials Materials in problem.
## @li \b bc Boundary conditions.
## @li \b interfaces Interior surfaces with constraints or
@@ -86,6 +87,12 @@
validator=pyre.inventory.choice([1,2,3]))
dimension.meta['tip'] = "Spatial dimension of problem space."
+ from spatialdata.units.Nondimensional import Nondimensional
+ normalizer = pyre.inventory.facility("normalizer",
+ family="nondimensional",
+ factory=Nondimensional)
+ normalizer.meta['tip'] = "Nondimensionalizer for problem."
+
from pylith.materials.Homogeneous import Homogeneous
materials = pyre.inventory.facilityArray("materials",
itemFactory=materialFactory,
@@ -215,6 +222,7 @@
Set members based using inventory.
"""
Component._configure(self)
+ self.normalizer = self.inventory.normalizer
self.dimension = self.inventory.dimension
self.materials = self.inventory.materials
self.bc = self.inventory.bc
@@ -258,4 +266,5 @@
"""
return Problem()
+
# End of file
Modified: short/3D/PyLith/trunk/pylith/problems/TimeDependent.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/TimeDependent.py 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/pylith/problems/TimeDependent.py 2008-11-29 23:58:23 UTC (rev 13421)
@@ -111,7 +111,8 @@
self._logger.eventBegin(logEvent)
self._info.log("Initializing problem.")
- self.formulation.initialize(self.dimension)
+ self.normalizer.initialize()
+ self.formulation.initialize(self.dimension, self.normalizer)
self._logger.eventEnd(logEvent)
return
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -23,6 +23,7 @@
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestBruneSlipFn );
@@ -410,12 +411,14 @@
ioPeakRate.filename(peakRateFilename);
dbPeakRate.ioHandler(&ioPeakRate);
+ spatialdata::units::Nondimensional normalizer;
+
// setup BruneSlipFn
slipfn->dbFinalSlip(&dbFinalSlip);
slipfn->dbSlipTime(&dbSlipTime);
slipfn->dbPeakRate(&dbPeakRate);
- slipfn->initialize(*faultMesh, &cs, originTime);
+ slipfn->initialize(*faultMesh, &cs, normalizer, originTime);
} // _initialize
// ----------------------------------------------------------------------
@@ -438,7 +441,8 @@
cs.setSpaceDim(spaceDim);
// Create fault mesh
- ALE::Obj<Mesh> faultMesh = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
+ ALE::Obj<Mesh> faultMesh =
+ new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
ALE::Obj<ALE::Mesh> faultBd = NULL;
const bool useLagrangeConstraints = true;
CohesiveTopology::createFault(faultMesh, faultBd,
@@ -449,7 +453,7 @@
data.faultId,
useLagrangeConstraints);
CPPUNIT_ASSERT(!faultMesh.isNull());
- // Need to copy coordinates from mesh to fault mesh since we are not
+ // Need to copy coordinates from mesh to fault mesh since we are
// using create() instead of createParallel().
faultMesh->setRealSection("coordinates",
mesh->getRealSection("coordinates"));
@@ -470,6 +474,8 @@
ioPeakRate.filename(data.peakRateFilename);
dbPeakRate.ioHandler(&ioPeakRate);
+ spatialdata::units::Nondimensional normalizer;
+
// setup BruneSlipFn
BruneSlipFn slipfn;
slipfn.dbFinalSlip(&dbFinalSlip);
@@ -478,7 +484,7 @@
const double originTime = 5.353;
- slipfn.initialize(faultMesh, &cs, originTime);
+ slipfn.initialize(faultMesh, &cs, normalizer, originTime);
const double tolerance = 1.0e-06;
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -23,6 +23,7 @@
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestConstRateSlipFn );
@@ -326,11 +327,13 @@
ioSlipTime.filename(slipTimeFilename);
dbSlipTime.ioHandler(&ioSlipTime);
+ spatialdata::units::Nondimensional normalizer;
+
// setup ConstRateSlipFn
slipfn->dbSlipRate(&dbSlipRate);
slipfn->dbSlipTime(&dbSlipTime);
- slipfn->initialize(*faultMesh, &cs, originTime);
+ slipfn->initialize(*faultMesh, &cs, normalizer, originTime);
} // _initialize
// ----------------------------------------------------------------------
@@ -385,9 +388,11 @@
slipfn.dbSlipRate(&dbSlipRate);
slipfn.dbSlipTime(&dbSlipTime);
+ spatialdata::units::Nondimensional normalizer;
+
const double originTime = 5.353;
- slipfn.initialize(faultMesh, &cs, originTime);
+ slipfn.initialize(faultMesh, &cs, normalizer, originTime);
const double tolerance = 1.0e-06;
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -24,6 +24,7 @@
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestEqKinSrc );
@@ -243,6 +244,8 @@
ioPeakRate.filename(peakRateFilename);
dbPeakRate.ioHandler(&ioPeakRate);
+ spatialdata::units::Nondimensional normalizer;
+
// setup EqKinSrc
slipfn->dbFinalSlip(&dbFinalSlip);
slipfn->dbSlipTime(&dbSlipTime);
@@ -250,7 +253,7 @@
eqsrc->originTime(originTime);
eqsrc->slipfn(slipfn);
- eqsrc->initialize(*faultMesh, &cs);
+ eqsrc->initialize(*faultMesh, &cs, normalizer);
} // _initialize
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -23,6 +23,7 @@
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestLiuCosSlipFn );
@@ -407,12 +408,14 @@
ioRiseTime.filename(riseTimeFilename);
dbRiseTime.ioHandler(&ioRiseTime);
+ spatialdata::units::Nondimensional normalizer;
+
// setup LiuCosSlipFn
slipfn->dbFinalSlip(&dbFinalSlip);
slipfn->dbSlipTime(&dbSlipTime);
slipfn->dbRiseTime(&dbRiseTime);
- slipfn->initialize(*faultMesh, &cs, originTime);
+ slipfn->initialize(*faultMesh, &cs, normalizer, originTime);
} // _initialize
// ----------------------------------------------------------------------
@@ -473,9 +476,11 @@
slipfn.dbSlipTime(&dbSlipTime);
slipfn.dbRiseTime(&dbRiseTime);
+ spatialdata::units::Nondimensional normalizer;
+
const double originTime = 5.353;
- slipfn.initialize(faultMesh, &cs, originTime);
+ slipfn.initialize(faultMesh, &cs, normalizer, originTime);
const double tolerance = 1.0e-06;
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -23,6 +23,7 @@
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestStepSlipFn );
@@ -317,11 +318,13 @@
ioSlipTime.filename(slipTimeFilename);
dbSlipTime.ioHandler(&ioSlipTime);
+ spatialdata::units::Nondimensional normalizer;
+
// setup StepSlipFn
slipfn->dbFinalSlip(&dbFinalSlip);
slipfn->dbSlipTime(&dbSlipTime);
- slipfn->initialize(*faultMesh, &cs, originTime);
+ slipfn->initialize(*faultMesh, &cs, normalizer, originTime);
} // _initialize
// ----------------------------------------------------------------------
@@ -376,9 +379,10 @@
slipfn.dbFinalSlip(&dbFinalSlip);
slipfn.dbSlipTime(&dbSlipTime);
+ spatialdata::units::Nondimensional normalizer;
const double originTime = 5.353;
- slipfn.initialize(faultMesh, &cs, originTime);
+ slipfn.initialize(faultMesh, &cs, normalizer, originTime);
const double tolerance = 1.0e-06;
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -25,6 +25,7 @@
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
#include "spatialdata/spatialdb/GravityField.hh" // USES GravityField
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include <math.h> // USES fabs()
@@ -325,10 +326,12 @@
spatialdata::spatialdb::SimpleDB db;
db.ioHandler(&iohandler);
+ spatialdata::units::Nondimensional normalizer;
+
_material->id(_data->matId);
_material->label(_data->matLabel);
_material->db(&db);
- _material->initialize(*mesh, &cs, _quadrature);
+ _material->initialize(*mesh, &cs, _quadrature, normalizer);
integrator->quadrature(_quadrature);
integrator->gravityField(_gravityField);
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc 2008-11-28 21:22:34 UTC (rev 13420)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc 2008-11-29 23:58:23 UTC (rev 13421)
@@ -26,6 +26,7 @@
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
#include "spatialdata/spatialdb/GravityField.hh" // USES GravityField
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include <math.h> // USES fabs()
@@ -326,10 +327,12 @@
spatialdata::spatialdb::SimpleDB db;
db.ioHandler(&iohandler);
+ spatialdata::units::Nondimensional normalizer;
+
_material->id(_data->matId);
_material->label(_data->matLabel);
_material->db(&db);
- _material->initialize(*mesh, &cs, _quadrature);
+ _material->initialize(*mesh, &cs, _quadrature, normalizer);
integrator->quadrature(_quadrature);
integrator->gravityField(_gravityField);
More information about the CIG-COMMITS
mailing list