[cig-commits] r11167 - in short/3D/PyLith/trunk: libsrc/faults unittests/libtests/faults unittests/libtests/faults/data unittests/pytests/bc
brad at geodynamics.org
brad at geodynamics.org
Fri Feb 15 13:37:20 PST 2008
Author: brad
Date: 2008-02-15 13:37:20 -0800 (Fri, 15 Feb 2008)
New Revision: 11167
Added:
short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichletBoundary.py
Modified:
short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc
short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh
short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.icc
short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.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/FaultCohesiveKin.hh
short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.cc
short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh
short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.hh
short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.hh
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataHex8.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4e.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4e.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3d.cc
Log:
Switched fault parameters to sections over fault mesh from sections over entire mesh.
Modified: short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -24,12 +24,22 @@
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
+namespace pylith {
+ namespace faults {
+ namespace _BruneSlipFn {
+ const int offsetPeakRate = 0;
+ const int offsetSlipTime = 1;
+ } // _BruneSlipFn
+ } // faults
+} // pylith
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::BruneSlipFn::BruneSlipFn(void) :
_dbFinalSlip(0),
_dbSlipTime(0),
- _dbPeakRate(0)
+ _dbPeakRate(0),
+ _spaceDim(0)
{ // constructor
} // constructor
@@ -45,61 +55,39 @@
// ----------------------------------------------------------------------
// Initialize slip time function.
void
-pylith::faults::BruneSlipFn::initialize(const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<Mesh>& faultMesh,
- const std::set<Mesh::point_type>& vertices,
-
- const spatialdata::geocoords::CoordSys* cs)
+pylith::faults::BruneSlipFn::initialize(
+ const ALE::Obj<Mesh>& faultMesh,
+ const spatialdata::geocoords::CoordSys* cs)
{ // initialize
- typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
-
- assert(!mesh.isNull());
+ assert(!faultMesh.isNull());
assert(0 != cs);
assert(0 != _dbFinalSlip);
assert(0 != _dbSlipTime);
assert(0 != _dbPeakRate);
- const int spaceDim = cs->spaceDim();
+ _spaceDim = cs->spaceDim();
+ const int spaceDim = _spaceDim;
+ const int indexFinalSlip = 0;
+ const int indexPeakRate = spaceDim + _BruneSlipFn::offsetPeakRate;
+ const int indexSlipTime = spaceDim + _BruneSlipFn::offsetSlipTime;
- // Create and allocate sections for parameters
- if (0 != _parameters) // Can't delete NULL pointer that holds reference
- delete _parameters;
- _parameters = new topology::FieldsManager(mesh);
- if (0 == _parameters)
- throw std::runtime_error("Could not create manager for parameters of "
- "Brune slip time function.");
- assert(0 != _parameters);
-
- // Parameter: final slip
- _parameters->addReal("final slip");
- const ALE::Obj<real_section_type>& finalSlip =
- _parameters->getReal("final slip");
- assert(!finalSlip.isNull());
- const vert_iterator vBegin = vertices.begin();
- const vert_iterator vEnd = vertices.end();
- for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter)
- finalSlip->setFiberDimension(*v_iter, spaceDim);
- mesh->allocate(finalSlip);
+ // Get vertices in fault mesh
+ const ALE::Obj<Mesh::label_sequence>& vertices = faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- // Parameter: slip initiation time
- _parameters->addReal("slip time");
- const ALE::Obj<real_section_type>& slipTime =
- _parameters->getReal("slip time");
- assert(!slipTime.isNull());
- // reuse atlas from finalSlip (but use local fiber dimension)
- for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter)
- slipTime->setFiberDimension(*v_iter, 1);
- mesh->allocate(slipTime);
- slipTime->getAtlas()->setAtlas(finalSlip->getAtlas()->getAtlas());
+ const int fiberDim = spaceDim + 2;
+ _parameters = new real_section_type(faultMesh->comm(), faultMesh->debug());
+ _parameters->addSpace(); // final slip
+ _parameters->addSpace(); // peak slip rate
+ _parameters->addSpace(); // slip time
+ assert(3 == _parameters->getNumSpaces());
+ _parameters->setFiberDimension(vertices, fiberDim);
+ _parameters->setFiberDimension(vertices, spaceDim, 0); // final slip
+ _parameters->setFiberDimension(vertices, 1, 1); // peak slip rate
+ _parameters->setFiberDimension(vertices, 1, 2); // slip time
+ faultMesh->allocate(_parameters);
+ assert(!_parameters.isNull());
- // Parameter: peak slip rate
- _parameters->addReal("peak rate");
- const ALE::Obj<real_section_type>& peakRate =
- _parameters->getReal("peak rate");
- assert(!peakRate.isNull());
- peakRate->setAtlas(slipTime->getAtlas()); // reuse atlas from slipTime
- peakRate->allocateStorage();
-
// Open databases and set query values
_dbFinalSlip->open();
switch (spaceDim)
@@ -134,51 +122,53 @@
// Get coordinates of vertices
const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
+ faultMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- // Query databases for parameters
- double_array slipData(spaceDim);
- double slipTimeData;
- double peakRateData;
- for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter) {
+ double_array paramsVertex(fiberDim);
+
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter) {
+
// Get coordinates of vertex
- const real_section_type::value_type* vCoords =
+ const real_section_type::value_type* coordsVertex =
coordinates->restrictPoint(*v_iter);
+ assert(0 != coordsVertex);
- int err = _dbFinalSlip->query(&slipData[0], spaceDim,
- vCoords, spaceDim, cs);
+ int err = _dbFinalSlip->query(¶msVertex[indexFinalSlip], spaceDim,
+ coordsVertex, spaceDim, cs);
if (err) {
std::ostringstream msg;
msg << "Could not find final slip at (";
for (int i=0; i < spaceDim; ++i)
- msg << " " << vCoords[i];
+ msg << " " << coordsVertex[i];
msg << ") using spatial database " << _dbFinalSlip->label() << ".";
throw std::runtime_error(msg.str());
} // if
- finalSlip->updatePoint(*v_iter, &slipData[0]);
-
- err = _dbSlipTime->query(&slipTimeData, 1, vCoords, spaceDim, cs);
+ err = _dbPeakRate->query(¶msVertex[indexPeakRate], 1,
+ coordsVertex, spaceDim, cs);
if (err) {
std::ostringstream msg;
- msg << "Could not find slip initiation time at (";
+ msg << "Could not find peak slip rate at (";
for (int i=0; i < spaceDim; ++i)
- msg << " " << vCoords[i];
- msg << ") using spatial database " << _dbSlipTime->label() << ".";
+ msg << " " << coordsVertex[i];
+ msg << ") using spatial database " << _dbPeakRate->label() << ".";
throw std::runtime_error(msg.str());
} // if
- slipTime->updatePoint(*v_iter, &slipTimeData);
- err = _dbPeakRate->query(&peakRateData, 1, vCoords, spaceDim, cs);
+ err = _dbSlipTime->query(¶msVertex[indexSlipTime], 1,
+ coordsVertex, spaceDim, cs);
if (err) {
std::ostringstream msg;
- msg << "Could not find peak slip rate at (";
+ msg << "Could not find slip initiation time at (";
for (int i=0; i < spaceDim; ++i)
- msg << " " << vCoords[i];
- msg << ") using spatial database " << _dbPeakRate->label() << ".";
+ msg << " " << coordsVertex[i];
+ msg << ") using spatial database " << _dbSlipTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
- peakRate->updatePoint(*v_iter, &peakRateData);
+
+ _parameters->updatePoint(*v_iter, ¶msVertex[0]);
} // for
// Close databases
@@ -187,65 +177,62 @@
_dbPeakRate->close();
// Allocate slip field
- _slipField = new real_section_type(mesh->comm(), mesh->debug());
- for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter)
- _slipField->setFiberDimension(*v_iter, spaceDim);
- mesh->allocate(_slipField);
+ _slip = new real_section_type(faultMesh->comm(), faultMesh->debug());
+ _slip->setFiberDimension(vertices, spaceDim);
+ faultMesh->allocate(_slip);
+ assert(!_slip.isNull());
} // initialize
// ----------------------------------------------------------------------
// Get slip on fault surface at time t.
const ALE::Obj<pylith::real_section_type>&
pylith::faults::BruneSlipFn::slip(const double t,
- const std::set<Mesh::point_type>& vertices)
+ const ALE::Obj<Mesh>& faultMesh)
{ // slip
- typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
+ assert(!_parameters.isNull());
+ assert(!_slip.isNull());
+ assert(!faultMesh.isNull());
- assert(0 != _parameters);
- assert(!_slipField.isNull());
+ const int spaceDim = _spaceDim;
+ const int indexFinalSlip = 0;
+ const int indexPeakRate = spaceDim + _BruneSlipFn::offsetPeakRate;
+ const int indexSlipTime = spaceDim + _BruneSlipFn::offsetSlipTime;
+
+ double_array slipValues(spaceDim);
- // Get parameters
- const ALE::Obj<real_section_type>& finalSlip =
- _parameters->getReal("final slip");
- assert(!finalSlip.isNull());
+ // Get vertices in fault mesh
+ const ALE::Obj<Mesh::label_sequence>& vertices = faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- const ALE::Obj<real_section_type>& slipTime =
- _parameters->getReal("slip time");
- assert(!slipTime.isNull());
+ int count = 0;
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter, ++count) {
+ const real_section_type::value_type* paramsVertex =
+ _parameters->restrictPoint(*v_iter);
+ assert(0 != paramsVertex);
- const ALE::Obj<real_section_type>& peakRate =
- _parameters->getReal("peak rate");
- assert(!peakRate.isNull());
+ const double* finalSlip = ¶msVertex[indexFinalSlip];
+ const double peakRate = paramsVertex[indexPeakRate];
+ const double slipTime = paramsVertex[indexSlipTime];
+
+ double finalSlipMag = 0.0;
+ for (int i=0; i < spaceDim; ++i)
+ finalSlipMag += finalSlip[i]*finalSlip[i];
+ finalSlipMag = sqrt(finalSlipMag);
- double_array slipValues(3);
- const vert_iterator vBegin = vertices.begin();
- const vert_iterator vEnd = vertices.end();
- const int vSize = vertices.size();
- for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter) {
- // Get values of parameters at vertex
- const int numSlipValues = finalSlip->getFiberDimension(*v_iter);
- const real_section_type::value_type* vFinalSlip =
- finalSlip->restrictPoint(*v_iter);
- const real_section_type::value_type* vSlipTime =
- slipTime->restrictPoint(*v_iter);
- const real_section_type::value_type* vPeakRate =
- peakRate->restrictPoint(*v_iter);
+ const double slip = _slipFn(t-slipTime, finalSlipMag, peakRate);
+ const double scale = finalSlipMag > 0.0 ? slip / finalSlipMag : 0.0;
+ for (int i=0; i < spaceDim; ++i)
+ slipValues[i] = scale * finalSlip[i];
- double vFinalSlipMag = 0.0;
- for (int iSlip=0; iSlip < numSlipValues; ++iSlip)
- vFinalSlipMag += vFinalSlip[iSlip]*vFinalSlip[iSlip];
- vFinalSlipMag = sqrt(vFinalSlipMag);
- const double vSlip = _slip(t-vSlipTime[0], vFinalSlipMag, vPeakRate[0]);
- const double scale = vFinalSlipMag > 0.0 ? vSlip / vFinalSlipMag : 0.0;
- for (int iSlip=0; iSlip < numSlipValues; ++iSlip)
- slipValues[iSlip] = scale * vFinalSlip[iSlip];
-
// Update field
- _slipField->updatePoint(*v_iter, &slipValues[0]);
+ _slip->updatePoint(*v_iter, &slipValues[0]);
} // for
- PetscLogFlopsNoCheck(vSize * (10 + 3*finalSlip->getFiberDimension(*vBegin)));
- return _slipField;
+ PetscLogFlopsNoCheck(count * (2+8 + 3*spaceDim));
+
+ return _slip;
} // slip
// ----------------------------------------------------------------------
@@ -253,56 +240,54 @@
const ALE::Obj<pylith::real_section_type>&
pylith::faults::BruneSlipFn::slipIncr(const double t0,
const double t1,
- const std::set<Mesh::point_type>& vertices)
+ const ALE::Obj<Mesh>& faultMesh)
{ // slipIncr
- typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
+ assert(!_parameters.isNull());
+ assert(!_slip.isNull());
+ assert(!faultMesh.isNull());
- assert(0 != _parameters);
- assert(!_slipField.isNull());
+ const int spaceDim = _spaceDim;
+ const int indexFinalSlip = 0;
+ const int indexPeakRate = spaceDim + _BruneSlipFn::offsetPeakRate;
+ const int indexSlipTime = spaceDim + _BruneSlipFn::offsetSlipTime;
+
+ double_array slipValues(spaceDim);
- // Get parameters
- const ALE::Obj<real_section_type>& finalSlip =
- _parameters->getReal("final slip");
- assert(!finalSlip.isNull());
+ // Get vertices in fault mesh
+ const ALE::Obj<Mesh::label_sequence>& vertices = faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- const ALE::Obj<real_section_type>& slipTime =
- _parameters->getReal("slip time");
- assert(!slipTime.isNull());
+ int count = 0;
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter, ++count) {
+ const real_section_type::value_type* paramsVertex =
+ _parameters->restrictPoint(*v_iter);
+ assert(0 != paramsVertex);
- const ALE::Obj<real_section_type>& peakRate =
- _parameters->getReal("peak rate");
- assert(!peakRate.isNull());
+ const double* finalSlip = ¶msVertex[indexFinalSlip];
+ const double peakRate = paramsVertex[indexPeakRate];
+ const double slipTime = paramsVertex[indexSlipTime];
+
+ double finalSlipMag = 0.0;
+ for (int i=0; i < spaceDim; ++i)
+ finalSlipMag += finalSlip[i]*finalSlip[i];
+ finalSlipMag = sqrt(finalSlipMag);
- double_array slipValues(3);
- const vert_iterator vBegin = vertices.begin();
- const vert_iterator vEnd = vertices.end();
- const int vSize = vertices.size();
- for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter) {
- // Get values of parameters at vertex
- const int numSlipValues = finalSlip->getFiberDimension(*v_iter);
- const real_section_type::value_type* vFinalSlip =
- finalSlip->restrictPoint(*v_iter);
- const real_section_type::value_type* vSlipTime =
- slipTime->restrictPoint(*v_iter);
- const real_section_type::value_type* vPeakRate =
- peakRate->restrictPoint(*v_iter);
+ const double slip0 = _slipFn(t0-slipTime, finalSlipMag, peakRate);
+ const double slip1 = _slipFn(t1-slipTime, finalSlipMag, peakRate);
+ const double scale = finalSlipMag > 0.0 ?
+ (slip1 - slip0) / finalSlipMag : 0.0;
+ for (int i=0; i < spaceDim; ++i)
+ slipValues[i] = scale * finalSlip[i];
- double vFinalSlipMag = 0.0;
- for (int iSlip=0; iSlip < numSlipValues; ++iSlip)
- vFinalSlipMag += vFinalSlip[iSlip]*vFinalSlip[iSlip];
- vFinalSlipMag = sqrt(vFinalSlipMag);
- const double vSlip0 = _slip(t0-vSlipTime[0], vFinalSlipMag, vPeakRate[0]);
- const double vSlip1 = _slip(t1-vSlipTime[0], vFinalSlipMag, vPeakRate[0]);
- const double scale = vFinalSlipMag > 0.0 ? (vSlip1 - vSlip0) / vFinalSlipMag : 0.0;
- for (int iSlip=0; iSlip < numSlipValues; ++iSlip)
- slipValues[iSlip] = scale * vFinalSlip[iSlip];
-
// Update field
- _slipField->updatePoint(*v_iter, &slipValues[0]);
+ _slip->updatePoint(*v_iter, &slipValues[0]);
} // for
- PetscLogFlopsNoCheck(vSize * (19 + 3*finalSlip->getFiberDimension(*vBegin)));
- return _slipField;
+ PetscLogFlopsNoCheck(count * (3+2*8 + 3*spaceDim));
+
+ return _slip;
} // slipIncr
// ----------------------------------------------------------------------
@@ -310,7 +295,7 @@
ALE::Obj<pylith::real_section_type>
pylith::faults::BruneSlipFn::finalSlip(void)
{ // finalSlip
- return _parameters->getReal("final slip");
+ return _parameters->getFibration(0);
} // finalSlip
// ----------------------------------------------------------------------
@@ -318,7 +303,7 @@
ALE::Obj<pylith::real_section_type>
pylith::faults::BruneSlipFn::slipTime(void)
{ // slipTime
- return _parameters->getReal("slip time");
+ return _parameters->getFibration(2);
} // slipTime
Modified: short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh 2008-02-15 21:37:20 UTC (rev 11167)
@@ -76,33 +76,29 @@
/** Initialize slip time function.
*
- * @param mesh Finite-element mesh.
* @param faultMesh Finite-element mesh of fault.
- * @param vertices Vertices where slip will be prescribed.
* @param cs Coordinate system for mesh
*/
- void initialize(const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<Mesh>& faultMesh,
- const std::set<Mesh::point_type>& vertices,
+ void initialize(const ALE::Obj<Mesh>& faultMesh,
const spatialdata::geocoords::CoordSys* cs);
/** Get slip on fault surface at time t.
*
* @param t Time t.
- * @param vertices Vertices where slip will be prescribed.
+ * @param faultMesh Mesh over fault surface.
*/
const ALE::Obj<real_section_type>& slip(const double t,
- const std::set<Mesh::point_type>& vertices);
+ const ALE::Obj<Mesh>& faultMesh);
/** Get slip increment on fault surface between time t0 and t1.
*
* @param t0 Time t.
* @param t1 Time t+dt.
- * @param vertices Vertices where slip will be prescribed.
+ * @param faultMesh Mesh over fault surface.
*/
const ALE::Obj<real_section_type>& slipIncr(const double t0,
const double t1,
- const std::set<Mesh::point_type>& vertices);
+ const ALE::Obj<Mesh>& faultMesh);
/** Get final slip.
*
@@ -137,15 +133,19 @@
* @returns Slip at point at time t
*/
static
- double _slip(const double t,
- const double finalSlip,
- const double peakRate);
+ double _slipFn(const double t,
+ const double finalSlip,
+ const double peakRate);
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
- ALE::Obj<real_section_type> _slipField; ///< Slip field on fault surface
+ /// Parameters for Brune slip time function.
+ /// Final slip (vector), peak slip rate (scalar), slip time (scalar).
+ ALE::Obj<real_section_type> _parameters;
+ ALE::Obj<real_section_type> _slip; ///< Slip field on fault surface
+
/// Spatial database for final slip
spatialdata::spatialdb::SpatialDB* _dbFinalSlip;
@@ -155,6 +155,8 @@
/// Spatial database for peak slip rate
spatialdata::spatialdb::SpatialDB* _dbPeakRate;
+ int _spaceDim; ///< Spatial dimension for slip field.
+
}; // class BruneSlipFn
#include "BruneSlipFn.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.icc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.icc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -41,9 +41,9 @@
// Compute slip using slip time function.
inline
double
-pylith::faults::BruneSlipFn::_slip(const double t,
- const double finalSlip,
- const double peakRate) {
+pylith::faults::BruneSlipFn::_slipFn(const double t,
+ const double finalSlip,
+ const double peakRate) {
double slip = 0.0;
if (t > 0) {
assert(peakRate > 0.0);
Modified: short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -643,13 +643,18 @@
// ----------------------------------------------------------------------
// Form a parallel fault mesh using the cohesive cell information
void
-pylith::faults::CohesiveTopology::createParallel(ALE::Obj<Mesh>* fault,
- const ALE::Obj<Mesh>& mesh,
- const int materialId,
- const bool constraintCell)
+pylith::faults::CohesiveTopology::createParallel(
+ ALE::Obj<Mesh>* fault,
+ std::map<Mesh::point_type, Mesh::point_type>* cohesiveToFault,
+ const ALE::Obj<Mesh>& mesh,
+ const int materialId,
+ const bool constraintCell)
{
assert(0 != fault);
+ assert(0 != cohesiveToFault);
+
*fault = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
+ cohesiveToFault->clear();
const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
const ALE::Obj<sieve_type> faultSieve =
@@ -695,11 +700,21 @@
for(int i=0; i < faceSize; ++i, ++v_iter)
faultSieve->addArrow(*v_iter, face, color++);
} // if/else
-
+ (*cohesiveToFault)[*c_iter] = face;
++face;
} // for
(*fault)->setSieve(faultSieve);
(*fault)->stratify();
+
+ const ALE::Obj<Mesh::label_sequence>& faultCells =
+ (*fault)->heightStratum(0);
+ assert(!faultCells.isNull());
+ for (Mesh::label_sequence::iterator c_iter=cBegin,
+ f_iter=faultCells->begin();
+ c_iter != cEnd;
+ ++c_iter, ++f_iter)
+ (*cohesiveToFault)[*c_iter] = *f_iter;
+
#if 1
(*fault)->setRealSection("coordinates", mesh->getRealSection("coordinates"));
#else
Modified: short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.hh 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.hh 2008-02-15 21:37:20 UTC (rev 11167)
@@ -58,14 +58,16 @@
/** Create (distributed) fault mesh from cohesive cells.
*
- * @param fault Finite-element mesh of fault (output)
- * @param mesh Finite-element mesh
+ * @param fault Finite-element mesh of fault (output).
+ * @param cohesiveToFault Mapping of cohesive cell to fault mesh cell.
+ * @param mesh Finite-element mesh.
* @param materialId Material id for cohesive elements.
* @param constraintCell True if creating cells constrained with
- * Lagrange multipliers that require extra vertices, false otherwise
+ * Lagrange multipliers that require extra vertices, false otherwise.
*/
static
void createParallel(ALE::Obj<Mesh>* fault,
+ std::map<Mesh::point_type, Mesh::point_type>* cohesiveToFault,
const ALE::Obj<Mesh>& mesh,
const int materialId,
const bool constraintCell =false);
Modified: short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -44,23 +44,21 @@
// Initialize slip time function.
void
pylith::faults::EqKinSrc::initialize(
- const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<Mesh>& faultMesh,
- const std::set<Mesh::point_type>& vertices,
- const spatialdata::geocoords::CoordSys* cs)
+ const ALE::Obj<Mesh>& faultMesh,
+ const spatialdata::geocoords::CoordSys* cs)
{ // initialize
assert(0 != _slipfn);
- _slipfn->initialize(mesh, faultMesh, vertices, cs);
+ _slipfn->initialize(faultMesh, cs);
} // initialize
// ----------------------------------------------------------------------
// Get slip on fault surface at time t.
const ALE::Obj<pylith::real_section_type>&
pylith::faults::EqKinSrc::slip(const double t,
- const std::set<Mesh::point_type>& vertices)
+ const ALE::Obj<Mesh>& faultMesh)
{ // slip
assert(0 != _slipfn);
- return _slipfn->slip(t, vertices);
+ return _slipfn->slip(t, faultMesh);
} // slip
// ----------------------------------------------------------------------
@@ -68,10 +66,10 @@
const ALE::Obj<pylith::real_section_type>&
pylith::faults::EqKinSrc::slipIncr(const double t0,
const double t1,
- const std::set<Mesh::point_type>& vertices)
+ const ALE::Obj<Mesh>& faultMesh)
{ // slip
assert(0 != _slipfn);
- return _slipfn->slipIncr(t0, t1, vertices);
+ return _slipfn->slipIncr(t0, t1, faultMesh);
} // slip
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh 2008-02-15 21:37:20 UTC (rev 11167)
@@ -68,36 +68,32 @@
/** Initialize slip time function.
*
- * @param mesh Finite-element mesh.
* @param faultMesh Finite-element mesh of fault.
- * @param vertices Vertices where slip will be prescribed.
* @param cs Coordinate system for mesh
*/
virtual
- void initialize(const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<Mesh>& faultMesh,
- const std::set<Mesh::point_type>& vertices,
+ void initialize(const ALE::Obj<Mesh>& faultMesh,
const spatialdata::geocoords::CoordSys* cs);
/** Get slip on fault surface at time t.
*
* @param t Time t.
- * @param vertices Vertices where slip will be prescribed.
+ * @param faultMesh Finite-element mesh of fault.
*/
virtual
const ALE::Obj<real_section_type>& slip(const double t,
- const std::set<Mesh::point_type>& vertices);
+ const ALE::Obj<Mesh>& faultMesh);
/** Get increment of slip on fault surface between time t0 and t1.
*
* @param t Time t.
- * @param vertices Vertices where slip will be prescribed.
+ * @param faultMesh Finite-element mesh of fault.
*/
virtual
const ALE::Obj<real_section_type>& slipIncr(
const double t0,
const double t1,
- const std::set<Mesh::point_type>& vertices);
+ const ALE::Obj<Mesh>& faultMesh);
/** Get final slip.
*
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -27,6 +27,7 @@
#include "spatialdata/spatialdb/SpatialDB.hh" // USES CoordSys
#include <Distribution.hh> // USES completeSection
+#include <Selection.hh> // Algorithms for submeshes
#include <math.h> // USES pow(), sqrt()
#include <assert.h> // USES assert()
@@ -74,286 +75,23 @@
throw std::runtime_error("Normal direction for fault orientation must be "
"a vector with 3 components.");
- CohesiveTopology::createParallel(&_faultMesh, mesh, id(),
- _useLagrangeConstraints());
+ CohesiveTopology::createParallel(&_faultMesh, &_cohesiveToFault, mesh, id(),
+ _useLagrangeConstraints());
- /* First find vertices associated with Lagrange multiplier
- * constraints in cohesive cells and compute the orientation of the
- * fault at these locations.
- */
+ //_faultMesh->view("FAULT MESH");
- const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
- assert(!sieve.isNull());
- const ALE::Obj<ALE::Mesh::label_sequence>& cellsCohesive =
- mesh->getLabelStratum("material-id", id());
- assert(!cellsCohesive.isNull());
- const Mesh::label_sequence::iterator cellsCohesiveBegin =
- cellsCohesive->begin();
- const Mesh::label_sequence::iterator cellsCohesiveEnd =
- cellsCohesive->end();
-
- // Create set of vertices associated with Lagrange multiplier constraints
- _constraintVert.clear();
- for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
- c_iter != cellsCohesiveEnd;
- ++c_iter) {
- // Vertices for each cohesive cell are in three groups of N.
- // 0 to N-1: vertices on negative side of the fault
- // N-1 to 2N-1: vertices on positive side of the fault
- // 2N to 3N-1: vertices associated with constraint forces
- const ALE::Obj<sieve_type::traits::coneSequence>& cone =
- sieve->cone(*c_iter);
- assert(!cone.isNull());
- const sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
- const int coneSize = cone->size();
- assert(coneSize % 3 == 0);
- sieve_type::traits::coneSequence::iterator v_iter = vBegin;
- // Skip over non-constraint vertices
- for (int i=0, numSkip=2*coneSize/3; i < numSkip; ++i)
- ++v_iter;
- // Add constraint vertices to set
- for(int i=0, numConstraintVert=coneSize/3;
- i < numConstraintVert;
- ++i, ++v_iter)
- _constraintVert.insert(*v_iter);
- } // for
-
- // Create orientation section for constraint vertices
- const int cohesiveDim = _faultMesh->getDimension();
- const int spaceDim = cs->spaceDim();
- const int orientationSize = spaceDim*spaceDim;
- _orientation = new real_section_type(mesh->comm(), mesh->debug());
- assert(!_orientation.isNull());
- for (int iDim=0; iDim <= cohesiveDim; ++iDim)
- _orientation->addSpace();
- assert(cohesiveDim+1 == _orientation->getNumSpaces());
- const std::set<Mesh::point_type>::const_iterator vertConstraintBegin =
- _constraintVert.begin();
- const std::set<Mesh::point_type>::const_iterator vertConstraintEnd =
- _constraintVert.end();
- const int vertConstraintSize = _constraintVert.size();
- for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
- v_iter != vertConstraintEnd;
- ++v_iter) {
- _orientation->setFiberDimension(*v_iter, orientationSize);
-
- // Create fibrations, one for each direction
- for (int iDim=0; iDim <= cohesiveDim; ++iDim)
- _orientation->setFiberDimension(*v_iter, spaceDim, iDim);
- } // for
- mesh->allocate(_orientation);
-
- // Compute orientation of fault at constraint vertices
-
- // Get section containing coordinates of vertices
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
-
- // Set orientation function
- assert(cohesiveDim == _quadrature->cellDim());
- assert(spaceDim == _quadrature->spaceDim());
-
- // Loop over cohesive cells, computing orientation at constraint vertices
- const int numBasis = _quadrature->numBasis();
- const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
- const double_array& verticesRef = cellGeometry.vertices();
- const int jacobianSize = (cohesiveDim > 0) ? spaceDim * cohesiveDim : 1;
- double_array jacobian(jacobianSize);
- double jacobianDet = 0;
- double_array vertexOrientation(orientationSize);
- double_array cohesiveVertices(3*numBasis*spaceDim);
- double_array faceVertices(numBasis*spaceDim);
-
- for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
- c_iter != cellsCohesiveEnd;
- ++c_iter) {
- mesh->restrict(coordinates, *c_iter,
- &cohesiveVertices[0], cohesiveVertices.size());
- const int size = numBasis*spaceDim;
- for (int i=0, offset=2*size; i < size; ++i)
- faceVertices[i] = cohesiveVertices[offset+i];
-
- const ALE::Obj<sieve_type::traits::coneSequence>& cone =
- sieve->cone(*c_iter);
- assert(!cone.isNull());
- const sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
- const sieve_type::traits::coneSequence::iterator vEnd = cone->end();
-
- // skip over non-constraint vertices
- sieve_type::traits::coneSequence::iterator vConstraintBegin = vBegin;
- const int numSkip = 2*numBasis;
- for (int i=0; i < numSkip; ++i)
- ++vConstraintBegin;
-
- int iBasis = 0;
- for(sieve_type::traits::coneSequence::iterator v_iter=vConstraintBegin;
- v_iter != vEnd;
- ++v_iter, ++iBasis) {
- // Compute Jacobian and determinant of Jacobian at vertex
- double_array vertex(&verticesRef[iBasis*cohesiveDim], cohesiveDim);
- cellGeometry.jacobian(&jacobian, &jacobianDet, faceVertices, vertex);
-
- // Compute orientation
- cellGeometry.orientation(&vertexOrientation, jacobian, jacobianDet,
- upDir);
-
- // Update orientation
- _orientation->updateAddPoint(*v_iter, &vertexOrientation[0]);
- } // for
- } // for
-
- // Assemble orientation information
- ALE::Distribution<Mesh>::completeSection(mesh, _orientation);
-
- // Loop over vertices, make orientation information unit magnitude
- double_array vertexDir(orientationSize);
- for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
- v_iter != vertConstraintEnd;
- ++v_iter) {
- const real_section_type::value_type* vertexOrient =
- _orientation->restrictPoint(*v_iter);
-
- assert(spaceDim*spaceDim == orientationSize);
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- double mag = 0;
- for (int jDim=0, index=iDim*spaceDim; jDim < spaceDim; ++jDim)
- mag += pow(vertexOrient[index+jDim],2);
- mag = sqrt(mag);
- assert(mag > 0.0);
- for (int jDim=0, index=iDim*spaceDim; jDim < spaceDim; ++jDim)
- vertexDir[index+jDim] =
- vertexOrient[index+jDim] / mag;
- } // for
-
- _orientation->updatePoint(*v_iter, &vertexDir[0]);
- } // for
- PetscLogFlopsNoCheck(vertConstraintSize * orientationSize * 4);
-
- if (2 == cohesiveDim) {
- // Check orientation of first vertex, if dot product of fault
- // normal with preferred normal is negative, flip up/down dip direction.
- // If the user gives the correct normal direction, we should end
- // up with left-lateral-slip, reverse-slip, and fault-opening for
- // positive slip values.
-
- const real_section_type::value_type* vertexOrient =
- _orientation->restrictPoint(*vertConstraintBegin);
- if (vertConstraintBegin != vertConstraintEnd)
- assert(9 == _orientation->getFiberDimension(*vertConstraintBegin));
- double_array vertNormalDir(&vertexOrient[6], 3);
- const double dot =
- normalDir[0]*vertNormalDir[0] +
- normalDir[1]*vertNormalDir[1] +
- normalDir[2]*vertNormalDir[2];
- if (dot < 0.0)
- for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
- v_iter != vertConstraintEnd;
- ++v_iter) {
- const real_section_type::value_type* vertexOrient =
- _orientation->restrictPoint(*v_iter);
- assert(9 == _orientation->getFiberDimension(*v_iter));
- // Keep along-strike direction
- for (int iDim=0; iDim < 3; ++iDim)
- vertexDir[iDim] = vertexOrient[iDim];
- // Flip up-dip direction
- for (int iDim=3; iDim < 6; ++iDim)
- vertexDir[iDim] = -vertexOrient[iDim];
- // Keep normal direction
- for (int iDim=6; iDim < 9; ++iDim)
- vertexDir[iDim] = vertexOrient[iDim];
-
- // Update direction
- _orientation->updatePoint(*v_iter, &vertexDir[0]);
- } // for
- PetscLogFlopsNoCheck(5 + vertConstraintSize * 3);
- } // if
-
- _eqsrc->initialize(mesh, _faultMesh, _constraintVert, cs);
-
// Establish pairing between constraint vertices and first cell they
// appear in to prevent overlap in integrating Jacobian
- const int noCell = -1;
- _constraintCell = new int_section_type(mesh->comm(), mesh->debug());
- assert(!_constraintCell.isNull());
- for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
- v_iter != vertConstraintEnd;
- ++v_iter)
- _constraintCell->setFiberDimension(*v_iter, 1);
- mesh->allocate(_constraintCell);
- // Set values to noCell
- for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
- v_iter != vertConstraintEnd;
- ++v_iter)
- _constraintCell->updatePoint(*v_iter, &noCell);
+ _calcVertexCellPairs();
- for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
- c_iter != cellsCohesiveEnd;
- ++c_iter) {
- const ALE::Obj<sieve_type::traits::coneSequence>& cone =
- sieve->cone(*c_iter);
- assert(!cone.isNull());
- const sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
- const sieve_type::traits::coneSequence::iterator vEnd = cone->end();
- const int coneSize = cone->size();
- assert(coneSize % 3 == 0);
- sieve_type::traits::coneSequence::iterator v_iter = vBegin;
- // Skip over non-constraint vertices
- for (int i=0, numSkip=2*coneSize/3; i < numSkip; ++i)
- ++v_iter;
- // If haven't set cell-constraint pair, then set it for current
- // cell, otherwise move on.
- for(int i=0, numConstraintVert=coneSize/3;
- i < numConstraintVert;
- ++i, ++v_iter) {
- const int_section_type::value_type* curCell =
- _constraintCell->restrictPoint(*v_iter);
- if (noCell == *curCell) {
- int point = *c_iter;
- _constraintCell->updatePoint(*v_iter, &point);
- } // if
- } // for
- } // for
-
// Setup pseudo-stiffness of cohesive cells to improve conditioning
// of Jacobian matrix
- _pseudoStiffness = new real_section_type(mesh->comm(), mesh->debug());
- assert(!_pseudoStiffness.isNull());
- for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
- v_iter != vertConstraintEnd;
- ++v_iter)
- _pseudoStiffness->setFiberDimension(*v_iter, 1);
- mesh->allocate(_pseudoStiffness);
-
- matDB->open();
- const char* stiffnessVals[] = { "density", "vs" };
- const int numStiffnessVals = 2;
- matDB->queryVals(stiffnessVals, numStiffnessVals);
-
- double_array matprops(numStiffnessVals);
- for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
- v_iter != vertConstraintEnd;
- ++v_iter) {
- const real_section_type::value_type* vertexCoords =
- coordinates->restrictPoint(*v_iter);
- int err = matDB->query(&matprops[0], numStiffnessVals, vertexCoords,
- spaceDim, cs);
- if (err) {
- std::ostringstream msg;
- msg << "Could not find material properties at (";
- for (int i=0; i < spaceDim; ++i)
- msg << " " << vertexCoords[i];
- msg << ") using spatial database " << matDB->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
-
- 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);
- } // for
- PetscLogFlopsNoCheck(vertConstraintSize * 2);
+ _calcConditioning(cs, matDB);
+
+ // Compute orientation at vertices in fault mesh.
+ _calcOrientation(upDir, normalDir);
+
+ _eqsrc->initialize(_faultMesh, cs);
} // initialize
// ----------------------------------------------------------------------
@@ -405,41 +143,42 @@
if (!_useSolnIncr) {
// Compute slip field at current time step
assert(0 != _eqsrc);
- _slip = _eqsrc->slip(t, _constraintVert);
+ _slip = _eqsrc->slip(t, _faultMesh);
assert(!_slip.isNull());
} else {
// Compute increment of slip field at current time step
assert(0 != _eqsrc);
- _slip = _eqsrc->slipIncr(t-_dt, t, _constraintVert);
+ _slip = _eqsrc->slipIncr(t-_dt, t, _faultMesh);
assert(!_slip.isNull());
} // else
for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
c_iter != cellsCohesiveEnd;
++c_iter) {
+ const Mesh::point_type c_fault = _cohesiveToFault[*c_iter];
+
cellResidual = 0.0;
- // Get constraint/cell pairings (only valid at constraint vertices)
- mesh->restrict(_constraintCell, *c_iter, &cellConstraintCell[0],
- cellConstraintCell.size());
+ // Get Lagrange constraint / fault cell pairings.
+ _faultMesh->restrict(_faultVertexCell, c_fault, &cellConstraintCell[0],
+ cellConstraintCell.size());
- // Get orientations at cells vertices (only valid at constraint vertices)
- mesh->restrict(_orientation, *c_iter, &cellOrientation[0],
+ // Get orientations at fault cell's vertices.
+ _faultMesh->restrict(_orientation, c_fault, &cellOrientation[0],
cellOrientation.size());
- // Get pseudo stiffness at cells vertices (only valid at
- // constraint vertices)
- mesh->restrict(_pseudoStiffness, *c_iter, &cellStiffness[0],
- cellStiffness.size());
+ // Get pseudo stiffness at fault cell's vertices.
+ _faultMesh->restrict(_pseudoStiffness, c_fault, &cellStiffness[0],
+ cellStiffness.size());
- // Get slip at cells vertices (only valid at constraint vertices)
- mesh->restrict(_slip, *c_iter, &cellSlip[0], cellSlip.size());
+ // Get slip at fault cell's vertices.
+ _faultMesh->restrict(_slip, c_fault, &cellSlip[0], cellSlip.size());
- // Get solution at cells vertices (valid at all cohesive vertices)
+ // Get solution at cohesive cell's vertices.
mesh->restrict(solution, *c_iter, &cellSoln[0], cellSoln.size());
for (int iConstraint=0; iConstraint < numConstraintVert; ++iConstraint) {
// Skip setting values if they are set by another cell
- if (cellConstraintCell[iConstraint] != *c_iter)
+ if (cellConstraintCell[iConstraint] != c_fault)
continue;
// Blocks in cell matrix associated with normal cohesive
@@ -460,7 +199,8 @@
// Get orientation at constraint vertex
const real_section_type::value_type* constraintOrient =
&cellOrientation[iConstraint*orientationSize];
-
+ assert(0 != constraintOrient);
+
// Entries associated with constraint forces applied at node i
for (int iDim=0; iDim < spaceDim; ++iDim)
for (int kDim=0; kDim < spaceDim; ++kDim)
@@ -545,23 +285,24 @@
for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
c_iter != cellsCohesiveEnd;
++c_iter) {
+ const Mesh::point_type c_fault = _cohesiveToFault[*c_iter];
+
cellMatrix = 0.0;
- // Get orientations at cells vertices (only valid at constraint vertices)
- mesh->restrict(_orientation, *c_iter, &cellOrientation[0],
- cellOrientation.size());
+ // Get Lagrange constraint / fault cell pairings.
+ _faultMesh->restrict(_faultVertexCell, c_fault, &cellConstraintCell[0],
+ cellConstraintCell.size());
- // Get pseudo stiffness at cells vertices (only valid at
- // constraint vertices)
- mesh->restrict(_pseudoStiffness, *c_iter, &cellStiffness[0],
- cellStiffness.size());
-
- // Get constraint/cell pairings (only valid at constraint vertices)
- mesh->restrict(_constraintCell, *c_iter, &cellConstraintCell[0],
- cellConstraintCell.size());
+ // Get orientations at fault cell's vertices.
+ _faultMesh->restrict(_orientation, c_fault, &cellOrientation[0],
+ cellOrientation.size());
+ // Get pseudo stiffness at fault cell's vertices.
+ _faultMesh->restrict(_pseudoStiffness, c_fault, &cellStiffness[0],
+ cellStiffness.size());
+
for (int iConstraint=0; iConstraint < numConstraintVert; ++iConstraint) {
// Skip setting values if they are set by another cell
- if (cellConstraintCell[iConstraint] != *c_iter)
+ if (cellConstraintCell[iConstraint] != c_fault)
continue;
// Blocks in cell matrix associated with normal cohesive
@@ -573,6 +314,7 @@
// Get orientation at constraint vertex
const real_section_type::value_type* constraintOrient =
&cellOrientation[iConstraint*orientationSize];
+ assert(0 != constraintOrient);
const double pseudoStiffness = cellStiffness[iConstraint];
@@ -637,6 +379,7 @@
<< "'.";
throw std::runtime_error(msg.str());
} // if
+
const int numCorners = _quadrature->refGeometry().numCorners();
const ALE::Obj<ALE::Mesh::label_sequence>& cells =
mesh->getLabelStratum("material-id", id());
@@ -676,54 +419,32 @@
const int cohesiveDim = _faultMesh->getDimension();
if (0 == strcasecmp("slip", name)) {
- _allocateBufferVertexVector();
- assert(!_slip.isNull());
- _projectCohesiveVertexField(&_bufferVertexVector, _slip, mesh);
*fieldType = VECTOR_FIELD;
- return _bufferVertexVector;
+ return _slip;
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
- _allocateBufferVertexVector();
- const ALE::Obj<real_section_type>& strikeDir =
- _orientation->getFibration(0);
- assert(!strikeDir.isNull());
- _projectCohesiveVertexField(&_bufferVertexVector, strikeDir, mesh);
+ _bufferVertexVector = _orientation->getFibration(0);
*fieldType = VECTOR_FIELD;
return _bufferVertexVector;
} else if (2 == cohesiveDim && 0 == strcasecmp("dip_dir", name)) {
- _allocateBufferVertexVector();
- const ALE::Obj<real_section_type>& dipDir =
- _orientation->getFibration(1);
- assert(!dipDir.isNull());
- _projectCohesiveVertexField(&_bufferVertexVector, dipDir, mesh);
+ _bufferVertexVector = _orientation->getFibration(1);
*fieldType = VECTOR_FIELD;
return _bufferVertexVector;
} else if (0 == strcasecmp("normal_dir", name)) {
- _allocateBufferVertexVector();
const int space =
(0 == cohesiveDim) ? 0 : (1 == cohesiveDim) ? 1 : 2;
- const ALE::Obj<real_section_type>& normalDir =
- _orientation->getFibration(space);
- assert(!normalDir.isNull());
- _projectCohesiveVertexField(&_bufferVertexVector, normalDir, mesh);
+ _bufferVertexVector = _orientation->getFibration(space);
*fieldType = VECTOR_FIELD;
return _bufferVertexVector;
} else if (0 == strcasecmp("final_slip", name)) {
- _allocateBufferVertexVector();
- const ALE::Obj<real_section_type>& finalSlip = _eqsrc->finalSlip();
- assert(!finalSlip.isNull());
- _projectCohesiveVertexField(&_bufferVertexVector, finalSlip, mesh);
+ _bufferVertexVector = _eqsrc->finalSlip();
*fieldType = VECTOR_FIELD;
return _bufferVertexVector;
-
} else if (0 == strcasecmp("slip_time", name)) {
- _allocateBufferVertexScalar();
- const ALE::Obj<real_section_type>& slipTime = _eqsrc->slipTime();
- assert(!slipTime.isNull());
- _projectCohesiveVertexField(&_bufferVertexScalar, slipTime, mesh);
+ _bufferVertexScalar = _eqsrc->slipTime();
*fieldType = SCALAR_FIELD;
return _bufferVertexScalar;
@@ -766,10 +487,295 @@
throw std::runtime_error(msg.str());
// Return generic section to satisfy member function definition.
- return _bufferCellScalar;
+ return _bufferCellVector;
} // cellField
// ----------------------------------------------------------------------
+// Calculate orientation at fault vertices.
+void
+pylith::faults::FaultCohesiveKin::_calcOrientation(const double_array& upDir,
+ const double_array& normalDir)
+{ // _calcOrientation
+ assert(!_faultMesh.isNull());
+
+ // Get vertices in fault mesh
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ _faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ // Create orientation section for fault (constraint) vertices
+ const int cohesiveDim = _faultMesh->getDimension();
+ const int spaceDim = cohesiveDim + 1;
+ const int orientationSize = spaceDim*spaceDim;
+ _orientation = new real_section_type(_faultMesh->comm(),
+ _faultMesh->debug());
+ assert(!_orientation.isNull());
+ for (int iDim=0; iDim <= cohesiveDim; ++iDim)
+ _orientation->addSpace();
+ assert(cohesiveDim+1 == _orientation->getNumSpaces());
+ _orientation->setFiberDimension(vertices, orientationSize);
+ for (int iDim=0; iDim <= cohesiveDim; ++iDim)
+ _orientation->setFiberDimension(vertices, spaceDim, iDim);
+ _faultMesh->allocate(_orientation);
+
+ // Compute orientation of fault at constraint vertices
+
+ // Get section containing coordinates of vertices
+ const ALE::Obj<real_section_type>& coordinates =
+ _faultMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
+ // Set orientation function
+ assert(cohesiveDim == _quadrature->cellDim());
+ assert(spaceDim == _quadrature->spaceDim());
+
+ // Loop over cohesive cells, computing orientation at constraint vertices
+ const int numBasis = _quadrature->numBasis();
+ const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
+ const double_array& verticesRef = cellGeometry.vertices();
+ const int jacobianSize = (cohesiveDim > 0) ? spaceDim * cohesiveDim : 1;
+ double_array jacobian(jacobianSize);
+ double jacobianDet = 0;
+ double_array vertexOrientation(orientationSize);
+ double_array faceVertices(numBasis*spaceDim);
+
+ // Get fault cells (1 dimension lower than top-level cells)
+ const ALE::Obj<Mesh::label_sequence>& cells =
+ _faultMesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ const ALE::Obj<sieve_type>& sieve = _faultMesh->getSieve();
+ assert(!sieve.isNull());
+ const int faultDepth = _faultMesh->depth(); // depth of fault cells
+ typedef ALE::SieveAlg<Mesh> SieveAlg;
+
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ _faultMesh->restrict(coordinates, *c_iter,
+ &faceVertices[0], faceVertices.size());
+
+ const ALE::Obj<SieveAlg::coneArray>& cone =
+ SieveAlg::nCone(_faultMesh, *c_iter, faultDepth);
+ assert(!cone.isNull());
+ const SieveAlg::coneArray::iterator vBegin = cone->begin();
+ const SieveAlg::coneArray::iterator vEnd = cone->end();
+
+ int iBasis = 0;
+ for(SieveAlg::coneArray::iterator v_iter=vBegin;
+ v_iter != vEnd;
+ ++v_iter, ++iBasis) {
+ // Compute Jacobian and determinant of Jacobian at vertex
+ double_array vertex(&verticesRef[iBasis*cohesiveDim], cohesiveDim);
+ cellGeometry.jacobian(&jacobian, &jacobianDet, faceVertices, vertex);
+
+ // Compute orientation
+ cellGeometry.orientation(&vertexOrientation, jacobian, jacobianDet,
+ upDir);
+
+ // Update orientation
+ _orientation->updateAddPoint(*v_iter, &vertexOrientation[0]);
+ } // for
+ } // for
+
+ // Assemble orientation information
+ ALE::Distribution<Mesh>::completeSection(_faultMesh, _orientation);
+
+ // Loop over vertices, make orientation information unit magnitude
+ double_array vertexDir(orientationSize);
+ int count = 0;
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter, ++count) {
+ const real_section_type::value_type* vertexOrient =
+ _orientation->restrictPoint(*v_iter);
+ assert(0 != vertexOrient);
+
+ assert(spaceDim*spaceDim == orientationSize);
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ double mag = 0;
+ for (int jDim=0, index=iDim*spaceDim; jDim < spaceDim; ++jDim)
+ mag += pow(vertexOrient[index+jDim],2);
+ mag = sqrt(mag);
+ assert(mag > 0.0);
+ for (int jDim=0, index=iDim*spaceDim; jDim < spaceDim; ++jDim)
+ vertexDir[index+jDim] =
+ vertexOrient[index+jDim] / mag;
+ } // for
+
+ _orientation->updatePoint(*v_iter, &vertexDir[0]);
+ } // for
+ PetscLogFlopsNoCheck(count * orientationSize * 4);
+
+ if (2 == cohesiveDim) {
+ // Check orientation of first vertex, if dot product of fault
+ // normal with preferred normal is negative, flip up/down dip direction.
+ // If the user gives the correct normal direction, we should end
+ // up with left-lateral-slip, reverse-slip, and fault-opening for
+ // positive slip values.
+
+ const real_section_type::value_type* vertexOrient =
+ _orientation->restrictPoint(*vertices->begin());
+ assert(0 != vertexOrient);
+
+ double_array vertNormalDir(&vertexOrient[6], 3);
+ const double dot =
+ normalDir[0]*vertNormalDir[0] +
+ normalDir[1]*vertNormalDir[1] +
+ normalDir[2]*vertNormalDir[2];
+ if (dot < 0.0)
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter) {
+ const real_section_type::value_type* vertexOrient =
+ _orientation->restrictPoint(*v_iter);
+ assert(0 != vertexOrient);
+ assert(9 == _orientation->getFiberDimension(*v_iter));
+ // Keep along-strike direction
+ for (int iDim=0; iDim < 3; ++iDim)
+ vertexDir[iDim] = vertexOrient[iDim];
+ // Flip up-dip direction
+ for (int iDim=3; iDim < 6; ++iDim)
+ vertexDir[iDim] = -vertexOrient[iDim];
+ // Keep normal direction
+ for (int iDim=6; iDim < 9; ++iDim)
+ vertexDir[iDim] = vertexOrient[iDim];
+
+ // Update direction
+ _orientation->updatePoint(*v_iter, &vertexDir[0]);
+ } // for
+
+ PetscLogFlopsNoCheck(5 + count * 3);
+ } // if
+
+ //_orientation->view("ORIENTATION");
+} // _calcOrientation
+
+// ----------------------------------------------------------------------
+// Calculate conditioning field.
+void
+pylith::faults::FaultCohesiveKin::_calcVertexCellPairs(void)
+{ // _calcVertexCellPairs
+ assert(!_faultMesh.isNull());
+
+ // Get vertices in fault mesh
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ _faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ const int noCell = -1;
+ _faultVertexCell = new int_section_type(_faultMesh->comm(),
+ _faultMesh->debug());
+ assert(!_faultVertexCell.isNull());
+ _faultVertexCell->setFiberDimension(vertices, 1);
+ _faultMesh->allocate(_faultVertexCell);
+
+ // Set values to noCell
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter)
+ _faultVertexCell->updatePoint(*v_iter, &noCell);
+
+ const ALE::Obj<Mesh::label_sequence>& cells =
+ _faultMesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+
+ const ALE::Obj<sieve_type>& sieve = _faultMesh->getSieve();
+ assert(!sieve.isNull());
+ const int faultDepth = _faultMesh->depth(); // depth of fault cells
+ typedef ALE::SieveAlg<Mesh> SieveAlg;
+
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ const ALE::Obj<SieveAlg::coneArray>& cone =
+ SieveAlg::nCone(_faultMesh, *c_iter, faultDepth);
+ assert(!cone.isNull());
+ const SieveAlg::coneArray::iterator vBegin = cone->begin();
+ const SieveAlg::coneArray::iterator vEnd = cone->end();
+ const int coneSize = cone->size();
+
+ // If haven't set cell-constraint pair, then set it for current
+ // cell, otherwise move on.
+ SieveAlg::coneArray::iterator v_iter = vBegin;
+ for(int i=0; i < coneSize; ++i, ++v_iter) {
+ const int_section_type::value_type* curCell =
+ _faultVertexCell->restrictPoint(*v_iter);
+ assert(0 != curCell);
+ if (noCell == *curCell) {
+ int point = *c_iter;
+ _faultVertexCell->updatePoint(*v_iter, &point);
+ } // if
+ } // for
+ } // for
+
+ //_faultVertexCell->view("VERTEX/CELL PAIRINGS");
+} // _calcVertexCallPairs
+
+// ----------------------------------------------------------------------
+void
+pylith::faults::FaultCohesiveKin::_calcConditioning(
+ const spatialdata::geocoords::CoordSys* cs,
+ spatialdata::spatialdb::SpatialDB* matDB)
+{ // _calcConditioning
+ assert(0 != cs);
+ assert(0 != matDB);
+ assert(!_faultMesh.isNull());
+
+ const int spaceDim = cs->spaceDim();
+
+ // Get vertices in fault mesh
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ _faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ _pseudoStiffness = new real_section_type(_faultMesh->comm(),
+ _faultMesh->debug());
+ assert(!_pseudoStiffness.isNull());
+ _pseudoStiffness->setFiberDimension(vertices, 1);
+ _faultMesh->allocate(_pseudoStiffness);
+
+ matDB->open();
+ const char* stiffnessVals[] = { "density", "vs" };
+ const int numStiffnessVals = 2;
+ matDB->queryVals(stiffnessVals, numStiffnessVals);
+
+ // Get section containing coordinates of vertices
+ const ALE::Obj<real_section_type>& coordinates =
+ _faultMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
+ double_array matprops(numStiffnessVals);
+ int count = 0;
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter, ++count) {
+ const real_section_type::value_type* vertexCoords =
+ coordinates->restrictPoint(*v_iter);
+ assert(0 != vertexCoords);
+ int err = matDB->query(&matprops[0], numStiffnessVals, vertexCoords,
+ spaceDim, cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find material properties at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << vertexCoords[i];
+ msg << ") using spatial database " << matDB->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ 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);
+ } // for
+ PetscLogFlopsNoCheck(count * 2);
+} // _calcConditioning
+
+// ----------------------------------------------------------------------
// Allocate scalar field for output of vertex information.
void
pylith::faults::FaultCohesiveKin::_allocateBufferVertexScalar(void)
@@ -803,22 +809,6 @@
} // _allocateBufferVertexVector
// ----------------------------------------------------------------------
-// Allocate scalar field for output of cell information.
-void
-pylith::faults::FaultCohesiveKin::_allocateBufferCellScalar(void)
-{ // _allocateBufferCellScalar
- const int fiberDim = 1;
- if (_bufferCellScalar.isNull()) {
- _bufferCellScalar = new real_section_type(_faultMesh->comm(),
- _faultMesh->debug());
- const ALE::Obj<Mesh::label_sequence>& cells =
- _faultMesh->heightStratum(0);
- _bufferCellScalar->setFiberDimension(cells, fiberDim);
- _faultMesh->allocate(_bufferCellScalar);
- } // if
-} // _allocateBufferCellScalar
-
-// ----------------------------------------------------------------------
// Allocate vector field for output of cell information.
void
pylith::faults::FaultCohesiveKin::_allocateBufferCellVector(void)
@@ -835,51 +825,5 @@
} // if
} // _allocateBufferCellVector
-// ----------------------------------------------------------------------
-// Project field defined over cohesive cells to fault mesh.
-void
-pylith::faults::FaultCohesiveKin::_projectCohesiveVertexField(
- ALE::Obj<real_section_type>* fieldFault,
- const ALE::Obj<real_section_type>& fieldCohesive,
- const ALE::Obj<Mesh>& mesh)
-{ // _projectCohesiveVertexField
- assert(0 != _quadrature);
- // Get cohesive cells
- const ALE::Obj<ALE::Mesh::label_sequence>& cellsCohesive =
- mesh->getLabelStratum("material-id", id());
- assert(!cellsCohesive.isNull());
- const Mesh::label_sequence::iterator cellsCohesiveBegin =
- cellsCohesive->begin();
- const Mesh::label_sequence::iterator cellsCohesiveEnd =
- cellsCohesive->end();
- assert(!fieldCohesive.isNull());
-
- // Get fault cells
- const ALE::Obj<Mesh::label_sequence>& cellsFault =
- _faultMesh->heightStratum(0);
- assert(!cellsFault.isNull());
- const Mesh::label_sequence::iterator cellsFaultBegin =
- cellsFault->begin();
- const Mesh::label_sequence::iterator cellsFaultEnd =
- cellsFault->end();
- assert(!fieldFault->isNull());
-
- // Project field at constraint (Lagrange multiplier) vertices
- // defined over cohesive cells to fault mesh
- const ALE::Obj<Mesh::label_sequence>& vertices =
- _faultMesh->depthStratum(0);
- const int fiberDim = fieldCohesive->getFiberDimension(*vertices->begin());
- const int numBasis = _quadrature->numBasis();
- double_array cellData(numBasis*fiberDim);
- for (Mesh::label_sequence::iterator c_iter=cellsCohesive->begin(),
- f_iter=cellsFault->begin();
- c_iter != cellsCohesiveEnd;
- ++c_iter, ++f_iter) {
- mesh->restrict(fieldCohesive, *c_iter, &cellData[0], cellData.size());
- _faultMesh->update(*fieldFault, *f_iter, &cellData[0]);
- } // for
-} // _projectCohesiveVertexField
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-02-15 21:37:20 UTC (rev 11167)
@@ -40,6 +40,8 @@
#include "FaultCohesive.hh" // ISA FaultCohesive
#include "pylith/feassemble/Integrator.hh" // ISA Integrator
+#include <map> // HASA std::map
+
/// Namespace for pylith package
namespace pylith {
namespace faults {
@@ -165,30 +167,41 @@
// PRIVATE METHODS ////////////////////////////////////////////////////
private :
+ /** Calculate orientation at fault vertices.
+ *
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction; only applies to fault surfaces in a 3-D domain).
+ * @param normalDir General preferred direction for fault normal
+ * (used to pick which of two possible normal directions for
+ * interface; only applies to fault surfaces in a 3-D domain).
+ */
+ void _calcOrientation(const double_array& upDir,
+ const double_array& normalDir);
+
+ /** Calculate pairing between fault vertices and first cell they
+ * appear in to prevent double counting in integrating Jacobian.
+ */
+ void _calcVertexCellPairs(void);
+
+ /** Calculate conditioning field.
+ *
+ * @param cs Coordinate system for mesh
+ * @param matDB Database of bulk elastic properties for fault region
+ * (used to improve conditioning of Jacobian matrix)
+ */
+ void _calcConditioning(const spatialdata::geocoords::CoordSys* cs,
+ spatialdata::spatialdb::SpatialDB* matDB);
+
/// Allocate scalar field for output of vertex information.
void _allocateBufferVertexScalar(void);
/// Allocate vector field for output of vertex information.
void _allocateBufferVertexVector(void);
- /// Allocate scalar field for output of cell information.
- void _allocateBufferCellScalar(void);
-
/// Allocate vector field for output of cell information.
void _allocateBufferCellVector(void);
- /** Project field defined over cohesive cells to fault mesh.
- *
- * @param fieldFault Field defined over fault mesh.
- * @param fieldCohesive Field defined over cohesive cells.
- * @param mesh PETSc mesh for problem.
- */
- void _projectCohesiveVertexField(
- ALE::Obj<real_section_type>* fieldFault,
- const ALE::Obj<real_section_type>& fieldCohesive,
- const ALE::Obj<Mesh>& mesh);
-
-
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
@@ -203,35 +216,33 @@
EqKinSrc* _eqsrc; ///< Kinematic earthquake source information
- /// Pseudo-stiffness for scaling constraint information to improve
- /// conditioning of Jacobian matrix
+ /// Field over fault mesh vertices of pseudo-stiffness values for
+ /// scaling constraint information to improve conditioning of
+ /// Jacobian matrix.
ALE::Obj<real_section_type> _pseudoStiffness;
- /// Orientation of fault surface at vertices (fiber dimension is
- /// nonzero only at constraint vertices)
+ /// Field over the fault mesh vertices of orientation of fault
+ /// surface.
ALE::Obj<real_section_type> _orientation;
- /// Vector field of current slip or slip increment.
+ /// Field over the fault mesh vertices of vector field of current
+ /// slip or slip increment.
ALE::Obj<real_section_type> _slip;
- /// Fault vertices associated with constraints
- std::set<Mesh::point_type> _constraintVert;
-
- /// Label of cell used to compute Jacobian for each constraint vertex (must
+ /// Label of cell used to compute Jacobian for each fault vertex (must
/// prevent overlap so that only 1 cell will contribute for
/// each vertex).
- ALE::Obj<int_section_type> _constraintCell;
+ ALE::Obj<int_section_type> _faultVertexCell;
- /// Scalar field for output of vertex information.
+ std::map<Mesh::point_type, Mesh::point_type> _cohesiveToFault;
+
+ /// Scalar field for vertex information over fault mesh.
ALE::Obj<real_section_type> _bufferVertexScalar;
- /// Vector field for output of vertex information.
+ /// Vector field for vertex information over fault mesh.
ALE::Obj<real_section_type> _bufferVertexVector;
- /// Scalar field for output of cell information.
- ALE::Obj<real_section_type> _bufferCellScalar;
-
- /// Vector field for output of cell information.
+ /// Vector field for cell information over fault mesh.
ALE::Obj<real_section_type> _bufferCellVector;
}; // class FaultCohesiveKin
Modified: short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -18,8 +18,7 @@
// ----------------------------------------------------------------------
// Default constructor.
-pylith::faults::SlipTimeFn::SlipTimeFn(void) :
- _parameters(0)
+pylith::faults::SlipTimeFn::SlipTimeFn(void)
{ // constructor
} // constructor
@@ -27,15 +26,7 @@
// Destructor.
pylith::faults::SlipTimeFn::~SlipTimeFn(void)
{ // destructor
- delete _parameters; _parameters = 0;
} // destructor
-// ----------------------------------------------------------------------
-// Copy constructor.
-pylith::faults::SlipTimeFn::SlipTimeFn(const SlipTimeFn& f) :
- _parameters(0)
-{ // copy constructor
-} // copy constructor
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh 2008-02-15 21:37:20 UTC (rev 11167)
@@ -28,10 +28,6 @@
class SlipTimeFn;
class TestSlipTimeFn; // unit testing
} // faults
-
- namespace topology {
- class FieldsManager; // HOLDSA FieldsManager
- } // feassemble
} // pylith
/// Namespace for spatialdata package
@@ -58,38 +54,38 @@
/** Initialize slip time function.
*
- * @param mesh Finite-element mesh.
* @param faultMesh Finite-element mesh of fault.
- * @param vertices Vertices where slip will be prescribed.
* @param cs Coordinate system for mesh
*/
virtual
- void initialize(const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<Mesh>& faultMesh,
- const std::set<Mesh::point_type>& vertices,
+ void initialize(const ALE::Obj<Mesh>& faultMesh,
const spatialdata::geocoords::CoordSys* cs) = 0;
/** Get slip on fault surface at time t.
*
* @param t Time t.
- * @param vertices Vertices where slip will be prescribed.
+ * @param faultMesh Mesh over fault surface.
+ *
* @returns Slip vector as left-lateral/reverse/normal.
*/
virtual
- const ALE::Obj<real_section_type>& slip(const double t,
- const std::set<Mesh::point_type>& vertices) = 0;
-
+ const ALE::Obj<real_section_type>&
+ slip(const double t,
+ const ALE::Obj<Mesh>& faultMesh) = 0;
+
/** Get slip increment on fault surface between time t0 and t1.
*
* @param t0 Time t.
* @param t1 Time t+dt.
- * @param vertices Vertices where slip will be prescribed.
+ * @param faultMesh Mesh over fault surface.
+ *
* @returns Increment in slip vector as left-lateral/reverse/normal.
*/
virtual
- const ALE::Obj<real_section_type>& slipIncr(const double t0,
- const double t1,
- const std::set<Mesh::point_type>& vertices) = 0;
+ const ALE::Obj<real_section_type>&
+ slipIncr(const double t0,
+ const double t1,
+ const ALE::Obj<Mesh>& faultMesh) = 0;
/** Get final slip.
*
@@ -114,12 +110,6 @@
/// Not implemented
const SlipTimeFn& operator=(const SlipTimeFn& f);
- // PROTECTED MEMBERS //////////////////////////////////////////////////
-protected :
-
- /// Parameters for slip time function
- topology::FieldsManager* _parameters;
-
}; // class SlipTimeFn
#endif // pylith_faults_sliptimefn_hh
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am 2008-02-15 21:37:20 UTC (rev 11167)
@@ -38,6 +38,8 @@
TestBoundary.cc \
test_faults.cc
+
+
noinst_HEADERS = \
TestBruneSlipFn.hh \
TestEqKinSrc.hh \
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -211,78 +211,31 @@
void
pylith::faults::TestBruneSlipFn::testSlip(void)
{ // testSlip
- const char* meshFilename = "data/tri3.mesh";
- const char* faultLabel = "fault";
- const int faultId = 2;
- const char* finalSlipFilename = "data/tri3_finalslipB.spatialdb";
- const char* slipTimeFilename = "data/tri3_sliptime.spatialdb";
- const char* peakRateFilename = "data/tri3_peakrate.spatialdb";
- const int constraintPts[] = { 3, 4 };
const double finalSlipE[] = { 2.3, 0.1,
0.0, 0.0};
const double slipTimeE[] = { 1.2, 1.3 };
const double peakRateE[] = { 1.4, 1.5 };
- const int numConstraintPts = 2;
- ALE::Obj<Mesh> mesh;
- meshio::MeshIOAscii meshIO;
- meshIO.filename(meshFilename);
- meshIO.debug(false);
- meshIO.interpolate(false);
- meshIO.read(&mesh);
- CPPUNIT_ASSERT(!mesh.isNull());
- const int spaceDim = mesh->getDimension();
- spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(spaceDim);
-
- // Create fault mesh
ALE::Obj<Mesh> faultMesh;
- const bool useLagrangeConstraints = true;
- CohesiveTopology::create(&faultMesh, mesh,
- mesh->getIntSection(faultLabel),
- faultId);
- CPPUNIT_ASSERT(!faultMesh.isNull());
-
- // Create set of constraint vertices
- std::set<Mesh::point_type> eqsrcVertices;
- for (int i=0; i < numConstraintPts; ++i)
- eqsrcVertices.insert(constraintPts[i]);
- CPPUNIT_ASSERT_EQUAL(numConstraintPts, int(eqsrcVertices.size()));
+ BruneSlipFn slipfn;
+ _initialize(&faultMesh, &slipfn);
- // Setup databases
- spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
- spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
- ioFinalSlip.filename(finalSlipFilename);
- dbFinalSlip.ioHandler(&ioFinalSlip);
-
- spatialdata::spatialdb::SimpleDB dbSlipTime("slip time");
- spatialdata::spatialdb::SimpleIOAscii ioSlipTime;
- ioSlipTime.filename(slipTimeFilename);
- dbSlipTime.ioHandler(&ioSlipTime);
-
- spatialdata::spatialdb::SimpleDB dbPeakRate("peak rate");
- spatialdata::spatialdb::SimpleIOAscii ioPeakRate;
- ioPeakRate.filename(peakRateFilename);
- dbPeakRate.ioHandler(&ioPeakRate);
+ const int spaceDim = faultMesh->getDimension() + 1;
- // setup BruneSlipFn
- BruneSlipFn slipFn;
- slipFn.dbFinalSlip(&dbFinalSlip);
- slipFn.dbSlipTime(&dbSlipTime);
- slipFn.dbPeakRate(&dbPeakRate);
-
- slipFn.initialize(mesh, faultMesh, eqsrcVertices, &cs);
-
const double t = 2.134;
- const ALE::Obj<real_section_type>& slip = slipFn.slip(t, eqsrcVertices);
+ const ALE::Obj<real_section_type>& slip = slipfn.slip(t, faultMesh);
CPPUNIT_ASSERT(!slip.isNull());
- int iPoint = 0;
const double tolerance = 1.0e-06;
- typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
- const vert_iterator vBegin = eqsrcVertices.begin();
- const vert_iterator vEnd = eqsrcVertices.end();
- for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter, ++iPoint) {
+
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ int iPoint = 0;
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter, ++iPoint) {
double slipMag = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim)
slipMag += pow(finalSlipE[iPoint*spaceDim+iDim], 2);
@@ -295,6 +248,8 @@
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
const real_section_type::value_type* vals =
slip->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
for (int iDim=0; iDim < fiberDim; ++iDim) {
const double slipE = finalSlipE[iPoint*spaceDim+iDim] * slipNorm;
CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, vals[iDim], tolerance);
@@ -307,80 +262,32 @@
void
pylith::faults::TestBruneSlipFn::testSlipIncr(void)
{ // testSlipIncr
- const char* meshFilename = "data/tri3.mesh";
- const char* faultLabel = "fault";
- const int faultId = 2;
- const char* finalSlipFilename = "data/tri3_finalslipB.spatialdb";
- const char* slipTimeFilename = "data/tri3_sliptime.spatialdb";
- const char* peakRateFilename = "data/tri3_peakrate.spatialdb";
- const int constraintPts[] = { 3, 4 };
const double finalSlipE[] = { 2.3, 0.1,
0.0, 0.0};
const double slipTimeE[] = { 1.2, 1.3 };
const double peakRateE[] = { 1.4, 1.5 };
- const int numConstraintPts = 2;
- ALE::Obj<Mesh> mesh;
- meshio::MeshIOAscii meshIO;
- meshIO.filename(meshFilename);
- meshIO.debug(false);
- meshIO.interpolate(false);
- meshIO.read(&mesh);
- CPPUNIT_ASSERT(!mesh.isNull());
- const int spaceDim = mesh->getDimension();
- spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(spaceDim);
-
- // Create fault mesh
ALE::Obj<Mesh> faultMesh;
- const bool useLagrangeConstraints = true;
- CohesiveTopology::create(&faultMesh, mesh,
- mesh->getIntSection(faultLabel),
- faultId);
- CPPUNIT_ASSERT(!faultMesh.isNull());
+ BruneSlipFn slipfn;
+ _initialize(&faultMesh, &slipfn);
- // Create set of constraint vertices
- std::set<Mesh::point_type> eqsrcVertices;
- for (int i=0; i < numConstraintPts; ++i)
- eqsrcVertices.insert(constraintPts[i]);
- CPPUNIT_ASSERT_EQUAL(numConstraintPts, int(eqsrcVertices.size()));
-
- // Setup databases
- spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
- spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
- ioFinalSlip.filename(finalSlipFilename);
- dbFinalSlip.ioHandler(&ioFinalSlip);
-
- spatialdata::spatialdb::SimpleDB dbSlipTime("slip time");
- spatialdata::spatialdb::SimpleIOAscii ioSlipTime;
- ioSlipTime.filename(slipTimeFilename);
- dbSlipTime.ioHandler(&ioSlipTime);
-
- spatialdata::spatialdb::SimpleDB dbPeakRate("peak rate");
- spatialdata::spatialdb::SimpleIOAscii ioPeakRate;
- ioPeakRate.filename(peakRateFilename);
- dbPeakRate.ioHandler(&ioPeakRate);
+ const int spaceDim = faultMesh->getDimension() + 1;
- // setup BruneSlipFn
- BruneSlipFn slipFn;
- slipFn.dbFinalSlip(&dbFinalSlip);
- slipFn.dbSlipTime(&dbSlipTime);
- slipFn.dbPeakRate(&dbPeakRate);
-
- slipFn.initialize(mesh, faultMesh, eqsrcVertices, &cs);
-
const double t0 = 1.234;
const double t1 = 3.635;
- const ALE::Obj<real_section_type>& slip =
- slipFn.slipIncr(t0, t1, eqsrcVertices);
+ const ALE::Obj<real_section_type>& slip = slipfn.slipIncr(t0, t1, faultMesh);
CPPUNIT_ASSERT(!slip.isNull());
- int iPoint = 0;
const double tolerance = 1.0e-06;
- typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
- const vert_iterator vBegin = eqsrcVertices.begin();
- const vert_iterator vEnd = eqsrcVertices.end();
- for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter, ++iPoint) {
+
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ int iPoint = 0;
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter, ++iPoint) {
double slipMag = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim)
slipMag += pow(finalSlipE[iPoint*spaceDim+iDim], 2);
@@ -397,6 +304,8 @@
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
const real_section_type::value_type* vals =
slip->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
for (int iDim=0; iDim < fiberDim; ++iDim) {
const double slipE =
finalSlipE[iPoint*spaceDim+iDim] * (slipNorm1-slipNorm0);
@@ -417,19 +326,80 @@
const double tau = finalSlip / (exp(1.0) * peakRate);
const double slipE = finalSlip * (1.0 - exp(-t/tau) * (1.0 + t/tau));
- double slip = BruneSlipFn::_slip(t, finalSlip, peakRate);
+ double slip = BruneSlipFn::_slipFn(t, finalSlip, peakRate);
const double tolerance = 1.0e-06;
CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, slip, tolerance);
- slip = BruneSlipFn::_slip(-0.5, finalSlip, peakRate);
+ slip = BruneSlipFn::_slipFn(-0.5, finalSlip, peakRate);
CPPUNIT_ASSERT_EQUAL(0.0, slip);
- slip = BruneSlipFn::_slip(1.0e+10, finalSlip, peakRate);
+ slip = BruneSlipFn::_slipFn(1.0e+10, finalSlip, peakRate);
CPPUNIT_ASSERT_DOUBLES_EQUAL(finalSlip, slip, tolerance);
} // testSlipTH
// ----------------------------------------------------------------------
+// Initialize BruneSlipFn.
+void
+pylith::faults::TestBruneSlipFn::_initialize(ALE::Obj<Mesh>* faultMesh,
+ BruneSlipFn* slipfn)
+{ // _initialize
+ assert(0 != slipfn);
+
+ const char* meshFilename = "data/tri3.mesh";
+ const char* faultLabel = "fault";
+ const int faultId = 2;
+ const char* finalSlipFilename = "data/tri3_finalslipB.spatialdb";
+ const char* slipTimeFilename = "data/tri3_sliptime.spatialdb";
+ const char* peakRateFilename = "data/tri3_peakrate.spatialdb";
+
+ ALE::Obj<Mesh> mesh;
+ meshio::MeshIOAscii meshIO;
+ meshIO.filename(meshFilename);
+ meshIO.debug(false);
+ meshIO.interpolate(false);
+ meshIO.read(&mesh);
+ CPPUNIT_ASSERT(!mesh.isNull());
+ const int spaceDim = mesh->getDimension();
+ spatialdata::geocoords::CSCart cs;
+ cs.setSpaceDim(spaceDim);
+
+ // Create fault mesh
+ const bool useLagrangeConstraints = true;
+ CohesiveTopology::create(faultMesh, mesh,
+ mesh->getIntSection(faultLabel),
+ faultId);
+ CPPUNIT_ASSERT(!faultMesh->isNull());
+ // Need to copy coordinates from mesh to fault mesh since we are not
+ // using create() instead of createParallel().
+ (*faultMesh)->setRealSection("coordinates",
+ mesh->getRealSection("coordinates"));
+
+ // Setup databases
+ spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
+ spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
+ ioFinalSlip.filename(finalSlipFilename);
+ dbFinalSlip.ioHandler(&ioFinalSlip);
+
+ spatialdata::spatialdb::SimpleDB dbSlipTime("slip time");
+ spatialdata::spatialdb::SimpleIOAscii ioSlipTime;
+ ioSlipTime.filename(slipTimeFilename);
+ dbSlipTime.ioHandler(&ioSlipTime);
+
+ spatialdata::spatialdb::SimpleDB dbPeakRate("peak rate");
+ spatialdata::spatialdb::SimpleIOAscii ioPeakRate;
+ ioPeakRate.filename(peakRateFilename);
+ dbPeakRate.ioHandler(&ioPeakRate);
+
+ // setup BruneSlipFn
+ slipfn->dbFinalSlip(&dbFinalSlip);
+ slipfn->dbSlipTime(&dbSlipTime);
+ slipfn->dbPeakRate(&dbPeakRate);
+
+ slipfn->initialize(*faultMesh, &cs);
+} // _initialize
+
+// ----------------------------------------------------------------------
// Test initialize().
void
pylith::faults::TestBruneSlipFn::_testInitialize(const _TestBruneSlipFn::DataStruct& data)
@@ -455,13 +425,11 @@
mesh->getIntSection(data.faultLabel),
data.faultId);
CPPUNIT_ASSERT(!faultMesh.isNull());
+ // Need to copy coordinates from mesh to fault mesh since we are not
+ // using create() instead of createParallel().
+ faultMesh->setRealSection("coordinates",
+ mesh->getRealSection("coordinates"));
- // Create set of constraint vertices
- std::set<Mesh::point_type> eqsrcVertices;
- for (int i=0; i < data.numConstraintPts; ++i)
- eqsrcVertices.insert(data.constraintPts[i]);
- CPPUNIT_ASSERT_EQUAL(data.numConstraintPts, int(eqsrcVertices.size()));
-
// Setup databases
spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -479,56 +447,40 @@
dbPeakRate.ioHandler(&ioPeakRate);
// setup BruneSlipFn
- BruneSlipFn slipFn;
- slipFn.dbFinalSlip(&dbFinalSlip);
- slipFn.dbSlipTime(&dbSlipTime);
- slipFn.dbPeakRate(&dbPeakRate);
+ BruneSlipFn slipfn;
+ slipfn.dbFinalSlip(&dbFinalSlip);
+ slipfn.dbSlipTime(&dbSlipTime);
+ slipfn.dbPeakRate(&dbPeakRate);
- slipFn.initialize(mesh, faultMesh, eqsrcVertices, &cs);
+ slipfn.initialize(faultMesh, &cs);
- // Check parameter sections
const double tolerance = 1.0e-06;
- CPPUNIT_ASSERT(0 != slipFn._parameters);
- const ALE::Obj<real_section_type>& finalSlip =
- slipFn._parameters->getReal("final slip");
- const ALE::Obj<real_section_type>& slipTime =
- slipFn._parameters->getReal("slip time");
- const ALE::Obj<real_section_type>& peakRate =
- slipFn._parameters->getReal("peak rate");
- CPPUNIT_ASSERT(!finalSlip.isNull());
- CPPUNIT_ASSERT(!slipTime.isNull());
- CPPUNIT_ASSERT(!peakRate.isNull());
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+
int iPoint = 0;
- const vert_iterator vBegin = eqsrcVertices.begin();
- const vert_iterator vEnd = eqsrcVertices.end();
- for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter, ++iPoint) {
- { // final slip
- const int fiberDim = finalSlip->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const real_section_type::value_type* vals =
- finalSlip->restrictPoint(*v_iter);
- for (int iDim=0; iDim < fiberDim; ++iDim)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+iDim],
- vals[iDim],
- tolerance);
- } // final slip
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter, ++iPoint) {
+ const int fiberDim = slipfn._parameters->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(spaceDim+2, fiberDim);
- { // slip time
- const int fiberDim = slipTime->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(1, fiberDim);
- const real_section_type::value_type* vals =
- slipTime->restrictPoint(*v_iter);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint], vals[0], tolerance);
- } // slip time
+ const real_section_type::value_type* vals =
+ slipfn._parameters->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
- { // peak rate
- const int fiberDim = peakRate->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(1, fiberDim);
- const real_section_type::value_type* vals =
- peakRate->restrictPoint(*v_iter);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(data.peakRateE[iPoint], vals[0], tolerance);
- } // peak rate
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+iDim],
+ vals[iDim],
+ tolerance);
+
+ const double peakRate = vals[spaceDim ];
+ const double slipTime = vals[spaceDim+1];
+
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.peakRateE[iPoint], peakRate, tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint], slipTime, tolerance);
} // for
} // _testInitialize
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.hh 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.hh 2008-02-15 21:37:20 UTC (rev 11167)
@@ -21,12 +21,15 @@
#if !defined(pylith_faults_testbruneslipfn_hh)
#define pylith_faults_testbruneslipfn_hh
+#include "pylith/utils/sievetypes.hh" // USES Mesh
+
#include <cppunit/extensions/HelperMacros.h>
/// Namespace for pylith package
namespace pylith {
namespace faults {
class TestBruneSlipFn;
+ class BruneSlipFn;
namespace _TestBruneSlipFn {
struct DataStruct;
@@ -90,10 +93,20 @@
// PRIVATE METHODS ////////////////////////////////////////////////////
private :
+ /** Initialize BruneSlipFn.
+ *
+ * @param faultMesh Fault mesh.
+ * @param slipfn Brune slip function.
+ */
+ static
+ void _initialize(ALE::Obj<Mesh>* faultMesh,
+ BruneSlipFn* slipfn);
+
/** Test intialize().
*
* @param data Data for initialization and testing of BruneSlipFn.
*/
+ static
void _testInitialize(const _TestBruneSlipFn::DataStruct& data);
}; // class TestBruneSlipFn
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -54,73 +54,11 @@
void
pylith::faults::TestEqKinSrc::testInitialize(void)
{ // testInitialize
- const char* meshFilename = "data/tri3.mesh";
- const char* faultLabel = "fault";
- const int faultId = 2;
- const char* finalSlipFilename = "data/tri3_finalslip.spatialdb";
- const char* slipTimeFilename = "data/tri3_sliptime.spatialdb";
- const char* peakRateFilename = "data/tri3_peakrate.spatialdb";
- const int constraintPts[] = { 3, 4 };
- const double finalSlipE[] = { 2.3, 0.1,
- 2.4, 0.2};
- const double slipTimeE[] = { 1.2, 1.3 };
- const double peakRateE[] = { 1.4, 1.5 };
- const int numConstraintPts = 2;
-
- typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
-
- // Setup mesh
- ALE::Obj<Mesh> mesh;
- meshio::MeshIOAscii meshIO;
- meshIO.filename(meshFilename);
- meshIO.debug(false);
- meshIO.interpolate(false);
- meshIO.read(&mesh);
- CPPUNIT_ASSERT(!mesh.isNull());
- const int spaceDim = mesh->getDimension();
- spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(spaceDim);
-
- // Create fault mesh
ALE::Obj<Mesh> faultMesh;
- const bool useLagrangeConstraints = true;
- CohesiveTopology::create(&faultMesh, mesh,
- mesh->getIntSection(faultLabel),
- faultId);
- CPPUNIT_ASSERT(!faultMesh.isNull());
-
- // Create set of constraint vertices
- std::set<Mesh::point_type> eqsrcVertices;
- for (int i=0; i < numConstraintPts; ++i)
- eqsrcVertices.insert(constraintPts[i]);
- CPPUNIT_ASSERT_EQUAL(numConstraintPts, int(eqsrcVertices.size()));
-
- // Setup databases
- spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
- spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
- ioFinalSlip.filename(finalSlipFilename);
- dbFinalSlip.ioHandler(&ioFinalSlip);
-
- spatialdata::spatialdb::SimpleDB dbSlipTime("slip time");
- spatialdata::spatialdb::SimpleIOAscii ioSlipTime;
- ioSlipTime.filename(slipTimeFilename);
- dbSlipTime.ioHandler(&ioSlipTime);
-
- spatialdata::spatialdb::SimpleDB dbPeakRate("peak rate");
- spatialdata::spatialdb::SimpleIOAscii ioPeakRate;
- ioPeakRate.filename(peakRateFilename);
- dbPeakRate.ioHandler(&ioPeakRate);
-
- // setup EqKinSrc
- BruneSlipFn slipFn;
- slipFn.dbFinalSlip(&dbFinalSlip);
- slipFn.dbSlipTime(&dbSlipTime);
- slipFn.dbPeakRate(&dbPeakRate);
-
EqKinSrc eqsrc;
- eqsrc.slipfn(&slipFn);
- eqsrc.initialize(mesh, faultMesh, eqsrcVertices, &cs);
-
+ BruneSlipFn slipfn;
+ _initialize(&faultMesh, &eqsrc, &slipfn);
+
// Don't have access to details of slip time function, so we can't
// check parameters. Have to rely on test of slip() for verification
// of results.
@@ -131,80 +69,32 @@
void
pylith::faults::TestEqKinSrc::testSlip(void)
{ // testSlip
- const char* meshFilename = "data/tri3.mesh";
- const char* faultLabel = "fault";
- const int faultId = 2;
- const char* finalSlipFilename = "data/tri3_finalslip.spatialdb";
- const char* slipTimeFilename = "data/tri3_sliptime.spatialdb";
- const char* peakRateFilename = "data/tri3_peakrate.spatialdb";
- const int constraintPts[] = { 3, 4 };
const double finalSlipE[] = { 2.3, 0.1,
2.4, 0.2};
const double slipTimeE[] = { 1.2, 1.3 };
const double peakRateE[] = { 1.4, 1.5 };
- const int numConstraintPts = 2;
- ALE::Obj<Mesh> mesh;
- meshio::MeshIOAscii meshIO;
- meshIO.filename(meshFilename);
- meshIO.debug(false);
- meshIO.interpolate(false);
- meshIO.read(&mesh);
- CPPUNIT_ASSERT(!mesh.isNull());
- const int spaceDim = mesh->getDimension();
- spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(spaceDim);
-
- // Create fault mesh
ALE::Obj<Mesh> faultMesh;
- const bool useLagrangeConstraints = true;
- CohesiveTopology::create(&faultMesh, mesh,
- mesh->getIntSection(faultLabel),
- faultId);
- CPPUNIT_ASSERT(!faultMesh.isNull());
-
- // Create set of constraint vertices
- std::set<Mesh::point_type> eqsrcVertices;
- for (int i=0; i < numConstraintPts; ++i)
- eqsrcVertices.insert(constraintPts[i]);
- CPPUNIT_ASSERT_EQUAL(numConstraintPts, int(eqsrcVertices.size()));
-
- // Setup databases
- spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
- spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
- ioFinalSlip.filename(finalSlipFilename);
- dbFinalSlip.ioHandler(&ioFinalSlip);
-
- spatialdata::spatialdb::SimpleDB dbSlipTime("slip time");
- spatialdata::spatialdb::SimpleIOAscii ioSlipTime;
- ioSlipTime.filename(slipTimeFilename);
- dbSlipTime.ioHandler(&ioSlipTime);
-
- spatialdata::spatialdb::SimpleDB dbPeakRate("peak rate");
- spatialdata::spatialdb::SimpleIOAscii ioPeakRate;
- ioPeakRate.filename(peakRateFilename);
- dbPeakRate.ioHandler(&ioPeakRate);
-
- // setup EqKinSrc
- BruneSlipFn slipFn;
- slipFn.dbFinalSlip(&dbFinalSlip);
- slipFn.dbSlipTime(&dbSlipTime);
- slipFn.dbPeakRate(&dbPeakRate);
-
EqKinSrc eqsrc;
- eqsrc.slipfn(&slipFn);
- eqsrc.initialize(mesh, faultMesh, eqsrcVertices, &cs);
+ BruneSlipFn slipfn;
+ _initialize(&faultMesh, &eqsrc, &slipfn);
+ const int spaceDim = faultMesh->getDimension() + 1;
+
const double t = 2.134;
- const ALE::Obj<real_section_type>& slip = eqsrc.slip(t, eqsrcVertices);
+ const ALE::Obj<real_section_type>& slip = eqsrc.slip(t, faultMesh);
CPPUNIT_ASSERT(!slip.isNull());
- int iPoint = 0;
const double tolerance = 1.0e-06;
- typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
- const vert_iterator vBegin = eqsrcVertices.begin();
- const vert_iterator vEnd = eqsrcVertices.end();
- for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter, ++iPoint) {
+
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ int iPoint = 0;
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter, ++iPoint) {
double slipMag = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim)
slipMag += pow(finalSlipE[iPoint*spaceDim+iDim], 2);
@@ -217,6 +107,8 @@
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
const real_section_type::value_type* vals =
slip->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
for (int iDim=0; iDim < fiberDim; ++iDim) {
const double slipE = finalSlipE[iPoint*spaceDim+iDim] * slipNorm;
CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, vals[iDim], tolerance);
@@ -229,18 +121,72 @@
void
pylith::faults::TestEqKinSrc::testSlipIncr(void)
{ // testSlip
+ const double finalSlipE[] = { 2.3, 0.1,
+ 2.4, 0.2};
+ const double slipTimeE[] = { 1.2, 1.3 };
+ const double peakRateE[] = { 1.4, 1.5 };
+
+ ALE::Obj<Mesh> faultMesh;
+ EqKinSrc eqsrc;
+ BruneSlipFn slipfn;
+ _initialize(&faultMesh, &eqsrc, &slipfn);
+
+ const int spaceDim = faultMesh->getDimension() + 1;
+
+ const double t0 = 1.234;
+ const double t1 = 2.525;
+ const ALE::Obj<real_section_type>& slip =
+ eqsrc.slipIncr(t0, t1, faultMesh);
+ CPPUNIT_ASSERT(!slip.isNull());
+
+ const double tolerance = 1.0e-06;
+
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ int iPoint = 0;
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter, ++iPoint) {
+ double slipMag = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ slipMag += pow(finalSlipE[iPoint*spaceDim+iDim], 2);
+ slipMag = sqrt(slipMag);
+ const double tau = slipMag / (exp(1.0) * peakRateE[iPoint]);
+ const double tRef = slipTimeE[iPoint];
+ const double slipNorm0 =
+ (t0 > tRef) ? 1.0 - exp(-(t0-tRef)/tau) * (1.0 + (t0-tRef)/tau) : 0.0;
+ const double slipNorm1 =
+ (t1 > tRef) ? 1.0 - exp(-(t1-tRef)/tau) * (1.0 + (t1-tRef)/tau) : 0.0;
+
+ const int fiberDim = slip->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
+ const real_section_type::value_type* vals =
+ slip->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ const double slipE =
+ finalSlipE[iPoint*spaceDim+iDim] * (slipNorm1 - slipNorm0);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, vals[iDim], tolerance);
+ } // for
+ } // for
+} // testSlipIncr
+
+// ----------------------------------------------------------------------
+// Initialize EqKinSrc.
+void
+pylith::faults::TestEqKinSrc::_initialize(ALE::Obj<Mesh>* faultMesh,
+ EqKinSrc* eqsrc,
+ BruneSlipFn* slipfn)
+{ // _initialize
const char* meshFilename = "data/tri3.mesh";
const char* faultLabel = "fault";
const int faultId = 2;
const char* finalSlipFilename = "data/tri3_finalslip.spatialdb";
const char* slipTimeFilename = "data/tri3_sliptime.spatialdb";
const char* peakRateFilename = "data/tri3_peakrate.spatialdb";
- const int constraintPts[] = { 3, 4 };
- const double finalSlipE[] = { 2.3, 0.1,
- 2.4, 0.2};
- const double slipTimeE[] = { 1.2, 1.3 };
- const double peakRateE[] = { 1.4, 1.5 };
- const int numConstraintPts = 2;
ALE::Obj<Mesh> mesh;
meshio::MeshIOAscii meshIO;
@@ -254,19 +200,16 @@
cs.setSpaceDim(spaceDim);
// Create fault mesh
- ALE::Obj<Mesh> faultMesh;
const bool useLagrangeConstraints = true;
- CohesiveTopology::create(&faultMesh, mesh,
+ CohesiveTopology::create(faultMesh, mesh,
mesh->getIntSection(faultLabel),
faultId);
- CPPUNIT_ASSERT(!faultMesh.isNull());
+ CPPUNIT_ASSERT(!faultMesh->isNull());
+ // Need to copy coordinates from mesh to fault mesh since we are not
+ // using create() instead of createParallel().
+ (*faultMesh)->setRealSection("coordinates",
+ mesh->getRealSection("coordinates"));
- // Create set of constraint vertices
- std::set<Mesh::point_type> eqsrcVertices;
- for (int i=0; i < numConstraintPts; ++i)
- eqsrcVertices.insert(constraintPts[i]);
- CPPUNIT_ASSERT_EQUAL(numConstraintPts, int(eqsrcVertices.size()));
-
// Setup databases
spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
@@ -284,49 +227,13 @@
dbPeakRate.ioHandler(&ioPeakRate);
// setup EqKinSrc
- BruneSlipFn slipFn;
- slipFn.dbFinalSlip(&dbFinalSlip);
- slipFn.dbSlipTime(&dbSlipTime);
- slipFn.dbPeakRate(&dbPeakRate);
+ slipfn->dbFinalSlip(&dbFinalSlip);
+ slipfn->dbSlipTime(&dbSlipTime);
+ slipfn->dbPeakRate(&dbPeakRate);
- EqKinSrc eqsrc;
- eqsrc.slipfn(&slipFn);
- eqsrc.initialize(mesh, faultMesh, eqsrcVertices, &cs);
-
- const double t0 = 1.234;
- const double t1 = 2.525;
- const ALE::Obj<real_section_type>& slip =
- eqsrc.slipIncr(t0, t1, eqsrcVertices);
- CPPUNIT_ASSERT(!slip.isNull());
+ eqsrc->slipfn(slipfn);
+ eqsrc->initialize(*faultMesh, &cs);
+} // _initialize
- int iPoint = 0;
- const double tolerance = 1.0e-06;
- typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
- const vert_iterator vBegin = eqsrcVertices.begin();
- const vert_iterator vEnd = eqsrcVertices.end();
- for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter, ++iPoint) {
- double slipMag = 0.0;
- for (int iDim=0; iDim < spaceDim; ++iDim)
- slipMag += pow(finalSlipE[iPoint*spaceDim+iDim], 2);
- slipMag = sqrt(slipMag);
- const double tau = slipMag / (exp(1.0) * peakRateE[iPoint]);
- const double tRef = slipTimeE[iPoint];
- const double slipNorm0 =
- (t0 > tRef) ? 1.0 - exp(-(t0-tRef)/tau) * (1.0 + (t0-tRef)/tau) : 0.0;
- const double slipNorm1 =
- (t1 > tRef) ? 1.0 - exp(-(t1-tRef)/tau) * (1.0 + (t1-tRef)/tau) : 0.0;
- const int fiberDim = slip->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const real_section_type::value_type* vals =
- slip->restrictPoint(*v_iter);
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- const double slipE =
- finalSlipE[iPoint*spaceDim+iDim] * (slipNorm1 - slipNorm0);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, vals[iDim], tolerance);
- } // for
- } // for
-} // testSlipIncr
-
-
// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.hh 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.hh 2008-02-15 21:37:20 UTC (rev 11167)
@@ -21,16 +21,20 @@
#if !defined(pylith_faults_testeqkinsrc_hh)
#define pylith_faults_testeqkinsrc_hh
+#include "pylith/utils/sievetypes.hh" // USES Mesh
+
#include <cppunit/extensions/HelperMacros.h>
/// Namespace for pylith package
namespace pylith {
namespace faults {
class TestEqKinSrc;
+ class EqKinSrc;
+ class BruneSlipFn;
namespace _TestEqKinSrc {
struct DataStruct;
- } // _BruneSlipTimeFn
+ } // _TestEqKinSrc
} // faults
} // pylith
@@ -70,6 +74,20 @@
/// slip().
void testSlipIncr(void);
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /** Initialize EqKinSrc.
+ *
+ * @param faultMesh Fault mesh.
+ * @param eqsrc Earthquake source.
+ * @param slipfn Slip time function.
+ */
+ static
+ void _initialize(ALE::Obj<Mesh>* faultMesh,
+ EqKinSrc* eqsrc,
+ BruneSlipFn* slipfn);
+
}; // class TestEqKinSrc
#endif // pylith_faults_testeqkinsrc_hh
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -116,33 +116,32 @@
FaultCohesiveKin fault;
_initialize(&mesh, &fault);
- // Check set of constraint vertices
- CPPUNIT_ASSERT_EQUAL(_data->numConstraintVert,
- int(fault._constraintVert.size()));
- const std::set<Mesh::point_type>::const_iterator vertConstraintBegin =
- fault._constraintVert.begin();
- const std::set<Mesh::point_type>::const_iterator vertConstraintEnd =
- fault._constraintVert.end();
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ fault._faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
int iVertex = 0;
- for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
- v_iter != vertConstraintEnd;
- ++v_iter, ++iVertex)
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter, ++iVertex) {
CPPUNIT_ASSERT_EQUAL(_data->constraintVertices[iVertex],
*v_iter);
+ } // for
+ CPPUNIT_ASSERT_EQUAL(_data->numConstraintVert, iVertex);
// Check orientation
- iVertex = 0;
const int cellDim = _data->cellDim;
const int spaceDim = _data->spaceDim;
const int orientationSize = spaceDim*spaceDim;
- for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
- v_iter != vertConstraintEnd;
+ iVertex = 0;
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
++v_iter, ++iVertex) {
const int fiberDim = fault._orientation->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(orientationSize, fiberDim);
const real_section_type::value_type* vertexOrient =
fault._orientation->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != vertexOrient);
+
const double tolerance = 1.0e-06;
for (int i=0; i < orientationSize; ++i) {
const int index = iVertex*orientationSize+i;
@@ -153,21 +152,21 @@
// Check pairing of constraint vertices with cells
iVertex = 0;
- for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
- v_iter != vertConstraintEnd;
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
++v_iter, ++iVertex) {
- const int fiberDim = fault._constraintCell->getFiberDimension(*v_iter);
+ const int fiberDim = fault._faultVertexCell->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(1, fiberDim);
const int_section_type::value_type* vertexCell =
- fault._constraintCell->restrictPoint(*v_iter);
+ fault._faultVertexCell->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != vertexCell);
CPPUNIT_ASSERT_EQUAL(_data->constraintCells[iVertex], vertexCell[0]);
} // for
// Check pseudoStiffness
iVertex = 0;
- for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
- v_iter != vertConstraintEnd;
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
++v_iter, ++iVertex) {
const int fiberDim = fault._pseudoStiffness->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(1, fiberDim);
@@ -238,9 +237,9 @@
CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
const real_section_type::value_type* vals =
residual->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
- const bool isConstraint =
- (fault._constraintVert.end() != fault._constraintVert.find(*v_iter));
+ const bool isConstraint = _isConstraintVertex(*v_iter);
const double pseudoStiffness =
(!isConstraint) ? _data->pseudoStiffness : 1.0;
for (int i=0; i < fiberDimE; ++i) {
@@ -273,9 +272,9 @@
CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
const real_section_type::value_type* vals =
residual->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
- const bool isConstraint =
- (fault._constraintVert.end() != fault._constraintVert.find(*v_iter));
+ const bool isConstraint = _isConstraintVertex(*v_iter);
const double pseudoStiffness =
(!isConstraint) ? _data->pseudoStiffness : 1.0;
for (int i=0; i < fiberDimE; ++i) {
@@ -464,5 +463,22 @@
} // catch
} // _initialize
+// ----------------------------------------------------------------------
+// Determine if vertex is a Lagrange multiplier constraint vertex.
+bool
+pylith::faults::TestFaultCohesiveKin::_isConstraintVertex(const int vertex) const
+{ // _isConstraintVertex
+ assert(0 != _data);
+ const int numConstraintVert = _data->numConstraintVert;
+ bool isFound = false;
+ for (int i=0; i < _data->numConstraintVert; ++i)
+ if (_data->constraintVertices[i] == vertex) {
+ isFound = true;
+ break;
+ } // if
+ return isFound;
+} // _isConstraintVertex
+
+
// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh 2008-02-15 21:37:20 UTC (rev 11167)
@@ -106,6 +106,15 @@
void _initialize(ALE::Obj<ALE::Mesh>* mesh,
FaultCohesiveKin* const fault) const;
+ /** Determine if vertex is a Lagrange multiplier constraint vertex.
+ *
+ * @param vertex Label of vertex.
+ *
+ * @returns True if vertex is a constraint vertex, false otherwise.
+ */
+ bool _isConstraintVertex(const int vertex) const;
+
+
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataHex8.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataHex8.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -141,7 +141,7 @@
};
const int pylith::faults::CohesiveKinDataHex8::_constraintCells[] = {
- 22, 22, 22, 22
+ 24, 24, 24, 24
};
const double pylith::faults::CohesiveKinDataHex8::_valsResidual[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -89,7 +89,7 @@
};
const int pylith::faults::CohesiveKinDataLine2::_constraintCells[] = {
- 7
+ 9
};
const double pylith::faults::CohesiveKinDataLine2::_valsResidualIncr[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -118,7 +118,7 @@
};
const int pylith::faults::CohesiveKinDataQuad4::_constraintCells[] = {
- 12, 12
+ 14, 14
};
const double pylith::faults::CohesiveKinDataQuad4::_valsResidual[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4e.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4e.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4e.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -143,7 +143,7 @@
};
const int pylith::faults::CohesiveKinDataQuad4e::_constraintCells[] = {
- 19, 19, 20
+ 23, 23, 24
};
const double pylith::faults::CohesiveKinDataQuad4e::_valsResidual[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -107,7 +107,7 @@
};
const int pylith::faults::CohesiveKinDataTet4::_constraintCells[] = {
- 13, 13, 13
+ 15, 15, 15
};
const double pylith::faults::CohesiveKinDataTet4::_valsResidual[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4e.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4e.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4e.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -111,7 +111,7 @@
};
const int pylith::faults::CohesiveKinDataTet4e::_constraintCells[] = {
- 18, 18, 18, 19
+ 22, 22, 22, 23
};
const double pylith::faults::CohesiveKinDataTet4e::_valsResidual[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4f.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -107,7 +107,7 @@
};
const int pylith::faults::CohesiveKinDataTet4f::_constraintCells[] = {
- 13, 13, 13
+ 15, 15, 15
};
const double pylith::faults::CohesiveKinDataTet4f::_valsResidual[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -129,7 +129,7 @@
};
const int pylith::faults::CohesiveKinDataTri3::_constraintCells[] = {
- 10, 10
+ 12, 12
};
const double pylith::faults::CohesiveKinDataTri3::_valsResidual[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3d.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3d.cc 2008-02-15 21:36:44 UTC (rev 11166)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3d.cc 2008-02-15 21:37:20 UTC (rev 11167)
@@ -139,7 +139,7 @@
};
const int pylith::faults::CohesiveKinDataTri3d::_constraintCells[] = {
- 16, 16, 17
+ 20, 20, 21
};
const double pylith::faults::CohesiveKinDataTri3d::_valsResidual[] = {
Added: short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichletBoundary.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichletBoundary.py (rev 0)
+++ short/3D/PyLith/trunk/unittests/pytests/bc/TestDirichletBoundary.py 2008-02-15 21:37:20 UTC (rev 11167)
@@ -0,0 +1,201 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file unittests/pytests/bc/TestDirichletBoundary.py
+
+## @brief Unit testing of DirichletBoundary object.
+
+import unittest
+
+from pylith.bc.DirichletBoundary import DirichletBoundary
+
+# ----------------------------------------------------------------------
+class TestDirichletBoundary(unittest.TestCase):
+ """
+ Unit testing of DirichletBoundary object.
+ """
+
+ def test_implementsConstraint(self):
+ """
+ Test to make sure DirichletBoundary satisfies constraint requirements.
+ """
+ bc = DirichletBoundary()
+ from pylith.feassemble.Constraint import implementsConstraint
+ self.failUnless(implementsConstraint(bc))
+ return
+
+
+ def test_constructor(self):
+ """
+ Test constructor.
+ """
+ from pylith.bc.DirichletBoundary import DirichletBoundary
+ bc = DirichletBoundary()
+ return
+
+
+ def test_initialize(self):
+ """
+ Test initialize().
+
+ WARNING: This is not a rigorous test of initialize() because we
+ don't verify the results.
+ """
+
+ (mesh, bc, fields) = self._initialize()
+
+ self.assertNotEqual(None, bc.cppHandle)
+
+ # We should really add something here to check to make sure things
+ # actually initialized correctly
+ return
+
+
+ def test_setConstraintSizes(self):
+ """
+ Test setConstraintSizes().
+
+ WARNING: This is not a rigorous test of setConstraintSizes() because we
+ don't verify the results.
+ """
+
+ (mesh, bc, fields) = self._initialize()
+ field = fields.getReal("field")
+ bc.setConstraintSizes(field)
+
+ # We should really add something here to check to make sure things
+ # actually initialized correctly
+ return
+
+
+ def test_setConstraints(self):
+ """
+ Test setConstraints().
+
+ WARNING: This is not a rigorous test of setConstraints() because we
+ don't verify the results.
+ """
+
+ (mesh, bc, fields) = self._initialize()
+ field = fields.getReal("field")
+ bc.setConstraintSizes(field)
+ mesh.allocateRealSection(field)
+ bc.setConstraints(field)
+
+ # We should really add something here to check to make sure things
+ # actually initialized correctly
+ return
+
+
+ def test_useSolnIncr(self):
+ """
+ Test useSolnIncr().
+ """
+ (mesh, bc, fields) = self._initialize()
+ bc.useSolnIncr(True)
+ return
+
+
+ def test_setField(self):
+ """
+ Test setField().
+
+ WARNING: This is not a rigorous test of setField() because we
+ don't verify the results.
+ """
+
+ (mesh, bc, fields) = self._initialize()
+ field = fields.getReal("field")
+ bc.setConstraintSizes(field)
+ mesh.allocateRealSection(field)
+ bc.setConstraints(field)
+ from pyre.units.time import second
+ t = 1.0*second
+ bc.setField(t, field)
+
+ # We should really add something here to check to make sure things
+ # actually initialized correctly
+ return
+
+
+ def test_finalize(self):
+ """
+ Test finalize().
+
+ WARNING: This is not a rigorous test of finalize() because we
+ neither set the input fields or verify the results.
+ """
+ (mesh, bc, fields) = self._initialize()
+ bc.finalize()
+
+ # We should really add something here to check to make sure things
+ # actually initialized correctly.
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _initialize(self):
+ """
+ Initialize DirichletBoundary boundary condition.
+ """
+ from pylith.bc.DirichletBoundary import DirichletBoundary
+ bc = DirichletBoundary()
+ bc._configure()
+ bc.output._configure()
+ bc.output.writer._configure()
+ bc.label = "bc"
+ bc.fixedDOF = [1]
+
+ from spatialdata.spatialdb.SimpleDB import SimpleDB
+ db = SimpleDB()
+ db._configure()
+ db.label = "TestDirichletBoundary tri3"
+ db.iohandler.filename = "data/tri3.spatialdb"
+ db.initialize()
+ bc.db = db
+
+ from pylith.bc.FixedDOFDB import FixedDOFDB
+ dbRate = FixedDOFDB()
+ dbRate._configure()
+ dbRate.label = "TestDirichletBoundary rate tri3"
+ dbRate.initialize()
+ bc.dbRate = dbRate
+
+ from spatialdata.geocoords.CSCart import CSCart
+ cs = CSCart()
+ cs.spaceDim = 2
+
+ from pylith.meshio.MeshIOAscii import MeshIOAscii
+ importer = MeshIOAscii()
+ importer.filename = "data/tri3.mesh"
+ importer.coordsys = cs
+ mesh = importer.read(debug=False, interpolate=False)
+
+ bc.preinitialize(mesh)
+ from pyre.units.time import second
+ bc.initialize(totalTime=0.0*second, numTimeSteps=1)
+
+ # Setup fields
+ from pylith.topology.FieldsManager import FieldsManager
+ fields = FieldsManager(mesh)
+ fields.addReal("field")
+ fields.setFiberDimension("field", cs.spaceDim)
+ fields.allocate("field")
+
+ import pylith.topology.topology as bindings
+ bindings.zeroRealSection(fields.getReal("field"))
+
+ return (mesh, bc, fields)
+
+
+# End of file
More information about the cig-commits
mailing list