[cig-commits] r14596 - in short/3D/PyLith/branches/pylith-swig/libsrc: . faults
brad at geodynamics.org
brad at geodynamics.org
Sat Apr 4 20:11:48 PDT 2009
Author: brad
Date: 2009-04-04 20:11:47 -0700 (Sat, 04 Apr 2009)
New Revision: 14596
Modified:
short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/EqKinSrc.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/EqKinSrc.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/SlipTimeFn.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/faultsfwd.hh
Log:
More work on updating fault stuff.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-04-05 03:11:47 UTC (rev 14596)
@@ -32,6 +32,9 @@
faults/SlipTimeFn.cc \
faults/StepSlipFn.cc \
faults/ConstRateSlipFn.cc \
+ faults/BruneSlipFn.cc \
+ faults/LiuCosSlipFn.cc \
+ faults/EqKinSrc.cc \
feassemble/CellGeometry.cc \
feassemble/Constraint.cc \
feassemble/GeometryPoint1D.cc \
@@ -94,22 +97,11 @@
utils/TestArray.cc
# faults/TopologyOps.cc \
-# faults/BruneSlipFn.cc \
# faults/CohesiveTopology.cc \
-# faults/EqKinSrc.cc \
# faults/FaultCohesive.cc \
# faults/FaultCohesiveKin.cc \
# faults/FaultCohesiveDyn.cc \
-# faults/LiuCosSlipFn.cc \
-# faults/StepSlipFn.cc \
# materials/GenMaxwellIsotropic3D.cc \
-# meshio/CellFilter.cc \
-# meshio/CellFilterAvg.cc \
-# meshio/DataWriter.cc \
-# meshio/DataWriterVTK.cc \
-# meshio/OutputManager.cc \
-# meshio/VertexFilter.cc \
-# meshio/VertexFilterVecNorm.cc \
# meshio/UCDFaultFile.cc \
# topology/Distributor.cc \
# topology/MeshRefiner.cc \
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.cc 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.cc 2009-04-05 03:11:47 UTC (rev 14596)
@@ -14,8 +14,9 @@
#include "BruneSlipFn.hh" // implementation of object methods
-#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
-#include "pylith/utils/array.hh" // USES double_array
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/Field.hh" // USES Field
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -25,22 +26,19 @@
#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
+// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
+typedef pylith::topology::SubMesh::SieveMesh::label_sequence label_sequence;
+typedef pylith::topology::SubMesh::RealSection RealSection;
// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::BruneSlipFn::BruneSlipFn(void) :
+ _slipTimeVertex(0),
+ _riseTimeVertex(0),
_dbFinalSlip(0),
_dbSlipTime(0),
- _dbPeakRate(0),
- _spaceDim(0)
+ _dbRiseTime(0)
{ // constructor
} // constructor
@@ -48,50 +46,68 @@
// Destructor.
pylith::faults::BruneSlipFn::~BruneSlipFn(void)
{ // destructor
- _dbFinalSlip = 0;
- _dbSlipTime = 0;
- _dbPeakRate = 0;
+ delete _parameters; _parameters = 0;
+ _dbFinalSlip = 0; // :TODO: Use shared pointer.
+ _dbSlipTime = 0; // :TODO: Use shared pointer.
+ _dbRiseTime = 0; // :TODO: Use shared pointer.
} // destructor
// ----------------------------------------------------------------------
// Initialize slip time function.
void
pylith::faults::BruneSlipFn::initialize(
- const ALE::Obj<Mesh>& faultMesh,
- const spatialdata::geocoords::CoordSys* cs,
- const spatialdata::units::Nondimensional& normalizer,
- const double originTime)
+ const topology::SubMesh& faultMesh,
+ const spatialdata::units::Nondimensional& normalizer,
+ const double originTime)
{ // initialize
- assert(!faultMesh.isNull());
- assert(0 != cs);
assert(0 != _dbFinalSlip);
assert(0 != _dbSlipTime);
- assert(0 != _dbPeakRate);
+ assert(0 != _dbRiseTime);
- _spaceDim = cs->spaceDim();
- const int spaceDim = _spaceDim;
- const int indexFinalSlip = 0;
- const int indexPeakRate = spaceDim + _BruneSlipFn::offsetPeakRate;
- const int indexSlipTime = spaceDim + _BruneSlipFn::offsetSlipTime;
+ const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
+ assert(0 != cs);
+ const int spaceDim = cs->spaceDim();
+ const double lengthScale = normalizer.lengthScale();
+ const double timeScale = normalizer.timeScale();
+
// 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<SieveMesh>& sieveMesh = faultMesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<label_sequence>& vertices = sieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const label_sequence::iterator verticesEnd = vertices->end();
- 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->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(), vertices->end()), *std::max_element(vertices->begin(), vertices->end())+1));
- _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());
+ delete _parameters; _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(faultMesh);
+ assert(0 != _parameters);
+ _parameters->add("final slip", "final slip");
+ topology::Field<topology::SubMesh>& finalSlip =
+ _parameters->get("final slip");
+ finalSlip.newSection(vertices, spaceDim);
+ finalSlip.allocate();
+ finalSlip.scale(lengthScale);
+ finalSlip.vectorFieldType(topology::FieldBase::VECTOR);
+ const ALE::Obj<RealSection>& finalSlipSection = finalSlip.section();
+ assert(!finalSlipSection.isNull());
+ _parameters->add("slip time", "slip time");
+ topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
+ slipTime.newSection(finalSlipSection->getChart(), 1);
+ slipTime.allocate();
+ slipTime.scale(timeScale);
+ slipTime.vectorFieldType(topology::FieldBase::SCALAR);
+ const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
+ assert(!slipTimeSection.isNull());
+
+ _parameters->add("rise time", "rise time");
+ topology::Field<topology::SubMesh>& riseTime = _parameters->get("rise time");
+ riseTime.newSection(slipTime);
+ riseTime.allocate();
+ riseTime.scale(timeScale);
+ riseTime.vectorFieldType(topology::FieldBase::SCALAR);
+ const ALE::Obj<RealSection>& riseTimeSection = riseTime.section();
+ assert(!riseTimeSection.isNull());
+
// Open databases and set query values
_dbFinalSlip->open();
switch (spaceDim)
@@ -120,199 +136,211 @@
const char* slipTimeValues[] = {"slip-time"};
_dbSlipTime->queryVals(slipTimeValues, 1);
- _dbPeakRate->open();
- const char* peakRateValues[] = {"slip-rate"};
- _dbPeakRate->queryVals(peakRateValues, 1);
+ _dbRiseTime->open();
+ const char* riseTimeValues[] = {"rise-time"};
+ _dbSlipTime->queryVals(slipTimeValues, 1);
// Get coordinates of vertices
- const ALE::Obj<real_section_type>& coordinates =
- faultMesh->getRealSection("coordinates");
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- const double lengthScale = normalizer.lengthScale();
- const double timeScale = normalizer.timeScale();
- const double velocityScale =
- normalizer.lengthScale() / normalizer.timeScale();
-
- double_array paramsVertex(fiberDim);
+ _slipVertex.resize(spaceDim);
double_array vCoordsGlobal(spaceDim);
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ for (label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter) {
coordinates->restrictPoint(*v_iter,
&vCoordsGlobal[0], vCoordsGlobal.size());
normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(),
lengthScale);
-
- int err = _dbFinalSlip->query(¶msVertex[indexFinalSlip], spaceDim,
- &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
+
+ int err = _dbFinalSlip->query(&_slipVertex[0], _slipVertex.size(),
+ &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
- msg << "Could not find final slip at (";
+ msg << "Could not find slip rate at (";
for (int i=0; i < spaceDim; ++i)
msg << " " << vCoordsGlobal[i];
msg << ") using spatial database " << _dbFinalSlip->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(¶msVertex[indexFinalSlip], spaceDim,
+ normalizer.nondimensionalize(&_slipVertex[0], _slipVertex.size(),
lengthScale);
- err = _dbPeakRate->query(¶msVertex[indexPeakRate], 1,
+
+ err = _dbSlipTime->query(&_slipTimeVertex, 1,
&vCoordsGlobal[0], vCoordsGlobal.size(), 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 << " " << vCoordsGlobal[i];
- msg << ") using spatial database " << _dbPeakRate->label() << ".";
+ msg << ") using spatial database " << _dbSlipTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(¶msVertex[indexPeakRate], 1,
- velocityScale);
+ normalizer.nondimensionalize(&_slipTimeVertex, 1, timeScale);
+ // add origin time to rupture time
+ _slipTimeVertex += originTime;
- err = _dbSlipTime->query(¶msVertex[indexSlipTime], 1,
+ err = _dbRiseTime->query(&_riseTimeVertex, 1,
&vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
- msg << "Could not find slip initiation time at (";
+ msg << "Could not find rise time at (";
for (int i=0; i < spaceDim; ++i)
msg << " " << vCoordsGlobal[i];
- msg << ") using spatial database " << _dbSlipTime->label() << ".";
+ msg << ") using spatial database " << _dbRiseTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(¶msVertex[indexSlipTime], 1,
- timeScale);
+ normalizer.nondimensionalize(&_riseTimeVertex, 1, timeScale);
- // add origin time to rupture time
- paramsVertex[indexSlipTime] += originTime;
-
- _parameters->updatePoint(*v_iter, ¶msVertex[0]);
+ finalSlipSection->updatePoint(*v_iter, &_slipVertex[0]);
+ slipTimeSection->updatePoint(*v_iter, &_slipTimeVertex);
+ riseTimeSection->updatePoint(*v_iter, &_riseTimeVertex);
} // for
// Close databases
_dbFinalSlip->close();
_dbSlipTime->close();
- _dbPeakRate->close();
+ _dbRiseTime->close();
} // initialize
// ----------------------------------------------------------------------
// Get slip on fault surface at time t.
void
-pylith::faults::BruneSlipFn::slip(const ALE::Obj<pylith::real_section_type>& slipField,
- const double t,
- const ALE::Obj<Mesh>& faultMesh)
+pylith::faults::BruneSlipFn::slip(topology::Field<topology::SubMesh>* slip,
+ const double t)
{ // slip
- assert(!_parameters.isNull());
- assert(!slipField.isNull());
- assert(!faultMesh.isNull());
+ assert(0 != slip);
+ assert(0 != _parameters);
- 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 vertices in fault mesh
- const ALE::Obj<Mesh::label_sequence>& vertices = faultMesh->depthStratum(0);
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- const int numVertices = vertices->size();
+ const ALE::Obj<SieveMesh>& sieveMesh = slip->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<label_sequence>& vertices = sieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const label_sequence::iterator verticesEnd = vertices->end();
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ // Get sections
+ const topology::Field<topology::SubMesh>& finalSlip =
+ _parameters->get("final slip");
+ const ALE::Obj<RealSection>& finalSlipSection = finalSlip.section();
+ assert(!finalSlipSection.isNull());
+ const topology::Field<topology::SubMesh>& slipTime =
+ _parameters->get("slip time");
+ const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
+ assert(!slipTimeSection.isNull());
+ const topology::Field<topology::SubMesh>& riseTime =
+ _parameters->get("rise time");
+ const ALE::Obj<RealSection>& riseTimeSection = riseTime.section();
+ assert(!riseTimeSection.isNull());
+ const ALE::Obj<RealSection>& slipSection = slip->section();
+ assert(!slipSection.isNull());
+
+ const int spaceDim = _slipVertex.size();
+ for (label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter) {
- const real_section_type::value_type* paramsVertex =
- _parameters->restrictPoint(*v_iter);
- assert(0 != paramsVertex);
+ finalSlipSection->restrictPoint(*v_iter, &_slipVertex[0],
+ _slipVertex.size());
+ slipTimeSection->restrictPoint(*v_iter, &_slipTimeVertex, 1);
+ riseTimeSection->restrictPoint(*v_iter, &_riseTimeVertex, 1);
- 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 += _slipVertex[i]*_slipVertex[i];
finalSlipMag = sqrt(finalSlipMag);
- const double slip = _slipFn(t-slipTime, finalSlipMag, peakRate);
+ const double slip = _slipFn(t-_slipTimeVertex, finalSlipMag,
+ _riseTimeVertex);
const double scale = finalSlipMag > 0.0 ? slip / finalSlipMag : 0.0;
- for (int i=0; i < spaceDim; ++i)
- slipValues[i] = scale * finalSlip[i];
-
+ _slipVertex *= scale;
+
// Update field
- slipField->updateAddPoint(*v_iter, &slipValues[0]);
+ slipSection->updateAddPoint(*v_iter, &_slipVertex[0]);
} // for
- PetscLogFlops(numVertices * (2+8 + 3*spaceDim));
+ PetscLogFlops(vertices->size() * (2+8 + 3*_slipVertex.size()));
} // slip
// ----------------------------------------------------------------------
// Get increment of slip on fault surface between time t0 and t1.
void
-pylith::faults::BruneSlipFn::slipIncr(const ALE::Obj<pylith::real_section_type>& slipField,
+pylith::faults::BruneSlipFn::slipIncr(
+ topology::Field<topology::SubMesh>* slip,
const double t0,
- const double t1,
- const ALE::Obj<Mesh>& faultMesh)
+ const double t1)
{ // slipIncr
- assert(!_parameters.isNull());
- assert(!slipField.isNull());
- assert(!faultMesh.isNull());
+ assert(0 != slip);
+ assert(0 != _parameters);
- 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 vertices in fault mesh
- const ALE::Obj<Mesh::label_sequence>& vertices = faultMesh->depthStratum(0);
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- const int numVertices = vertices->size();
+ const ALE::Obj<SieveMesh>& sieveMesh = slip->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<label_sequence>& vertices = sieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const label_sequence::iterator verticesEnd = vertices->end();
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ // Get sections
+ const topology::Field<topology::SubMesh>& finalSlip =
+ _parameters->get("final slip");
+ const ALE::Obj<RealSection>& finalSlipSection = finalSlip.section();
+ assert(!finalSlipSection.isNull());
+ const topology::Field<topology::SubMesh>& slipTime =
+ _parameters->get("slip time");
+ const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
+ assert(!slipTimeSection.isNull());
+ const topology::Field<topology::SubMesh>& riseTime =
+ _parameters->get("rise time");
+ const ALE::Obj<RealSection>& riseTimeSection = riseTime.section();
+ assert(!riseTimeSection.isNull());
+ const ALE::Obj<RealSection>& slipSection = slip->section();
+ assert(!slipSection.isNull());
+
+ const int spaceDim = _slipVertex.size();
+ for (label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter) {
- const real_section_type::value_type* paramsVertex =
- _parameters->restrictPoint(*v_iter);
- assert(0 != paramsVertex);
+ finalSlipSection->restrictPoint(*v_iter, &_slipVertex[0],
+ _slipVertex.size());
+ slipTimeSection->restrictPoint(*v_iter, &_slipTimeVertex, 1);
+ riseTimeSection->restrictPoint(*v_iter, &_riseTimeVertex, 1);
- 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 += _slipVertex[i]*_slipVertex[i];
finalSlipMag = sqrt(finalSlipMag);
- const double slip0 = _slipFn(t0-slipTime, finalSlipMag, peakRate);
- const double slip1 = _slipFn(t1-slipTime, finalSlipMag, peakRate);
+ const double slip0 = _slipFn(t0-_slipTimeVertex, finalSlipMag,
+ _riseTimeVertex);
+ const double slip1 = _slipFn(t1-_slipTimeVertex, finalSlipMag,
+ _riseTimeVertex);
const double scale = finalSlipMag > 0.0 ?
(slip1 - slip0) / finalSlipMag : 0.0;
- for (int i=0; i < spaceDim; ++i)
- slipValues[i] = scale * finalSlip[i];
+ _slipVertex *= scale;
+
// Update field
- slipField->updateAddPoint(*v_iter, &slipValues[0]);
+ slipSection->updateAddPoint(*v_iter, &_slipVertex[0]);
} // for
- PetscLogFlops(numVertices * (3+2*8 + 3*spaceDim));
+ PetscLogFlops(vertices->size() * (3+2*8 + 3*_slipVertex.size()));
} // slipIncr
// ----------------------------------------------------------------------
// Get final slip.
-ALE::Obj<pylith::real_section_type>
+const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::BruneSlipFn::finalSlip(void)
{ // finalSlip
- return _parameters->getFibration(0);
+ return _parameters->get("final slip");
} // finalSlip
// ----------------------------------------------------------------------
// Get time when slip begins at each point.
-ALE::Obj<pylith::real_section_type>
+const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::BruneSlipFn::slipTime(void)
{ // slipTime
- return _parameters->getFibration(2);
+ return _parameters->get("slip time");
} // slipTime
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.hh 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.hh 2009-04-05 03:11:47 UTC (rev 14596)
@@ -14,7 +14,7 @@
*
* @brief C++ implementation of Brune slip time function.
*
- * Slip time function follows the integral of Brune's (1907) far-field
+ * Slip time function follows the integral of Brune's (1970) far-field
* time function.
*
* Normalize slip = 1 - exp(-t/tau)(1 + t/tau),
@@ -31,7 +31,7 @@
#include "pylith/utils/array.hh" // HASA double_array
-// SlipTimeFn -----------------------------------------------------------
+// BruneSlipFn ----------------------------------------------------------
class pylith::faults::BruneSlipFn : public SlipTimeFn
{ // class BruneSlipFn
friend class TestBruneSlipFn; // unit testing
@@ -43,7 +43,6 @@
BruneSlipFn(void);
/// Destructor.
- virtual
~BruneSlipFn(void);
/** Set spatial database for final slip.
@@ -136,7 +135,7 @@
double _slipTimeVertex; ///< Slip time at a vertex.
double _riseTimeVertex; ///< Rise time at a vertex.
- double_array _finalSlipVertex; ///< Final slip at a vertex.
+ double_array _slipVertex; ///< Slip at a vertex.
/// Parameters for Brune slip time function, final slip (vector),
/// rise time (scalar), slip time (scalar).
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.cc 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.cc 2009-04-05 03:11:47 UTC (rev 14596)
@@ -46,8 +46,8 @@
pylith::faults::ConstRateSlipFn::~ConstRateSlipFn(void)
{ // destructor
delete _parameters; _parameters = 0;
- _dbSlipRate = 0;
- _dbSlipTime = 0;
+ _dbSlipRate = 0; // :TODO: Use shared pointer.
+ _dbSlipTime = 0; // :TODO: Use shared pointer.
} // destructor
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.hh 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.hh 2009-04-05 03:11:47 UTC (rev 14596)
@@ -30,7 +30,7 @@
#include "pylith/utils/array.hh" // HASA double_array
-// SlipTimeFn -----------------------------------------------------------
+// ConstRateTimeFn ------------------------------------------------------
class pylith::faults::ConstRateSlipFn : public SlipTimeFn
{ // class ConstRateSlipFn
friend class TestConstRateSlipFn; // unit testing
@@ -42,7 +42,6 @@
ConstRateSlipFn(void);
/// Destructor.
- virtual
~ConstRateSlipFn(void);
/** Set spatial database for slip rate.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/EqKinSrc.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/EqKinSrc.cc 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/EqKinSrc.cc 2009-04-05 03:11:47 UTC (rev 14596)
@@ -32,7 +32,7 @@
// Destructor.
pylith::faults::EqKinSrc::~EqKinSrc(void)
{ // destructor
- _slipfn = 0; // Don't manage memory for slip fn
+ _slipfn = 0; // :TODO: Use shared pointer.
} // destructor
// ----------------------------------------------------------------------
@@ -56,49 +56,49 @@
void
pylith::faults::EqKinSrc::slipfn(SlipTimeFn* slipfn)
{ // slipfn
- _slipfn = slipfn; // Don't manage memory for slip fn
+ _slipfn = slipfn; // :TODO: Use shared pointer.
} // slipfn
// ----------------------------------------------------------------------
// Initialize slip time function.
void
pylith::faults::EqKinSrc::initialize(
- const ALE::Obj<Mesh>& faultMesh,
- const spatialdata::geocoords::CoordSys* cs,
- const spatialdata::units::Nondimensional& normalizer)
+ const topology::SubMesh& faultMesh,
+ const spatialdata::units::Nondimensional& normalizer)
{ // initialize
+ // :TODO: Normalize slip time in Python?
normalizer.nondimensionalize(&_originTime, 1, normalizer.timeScale());
assert(0 != _slipfn);
- _slipfn->initialize(faultMesh, cs, normalizer, _originTime);
+ _slipfn->initialize(faultMesh, normalizer, _originTime);
} // initialize
// ----------------------------------------------------------------------
// Get slip on fault surface at time t.
void
-pylith::faults::EqKinSrc::slip(const ALE::Obj<pylith::real_section_type>& slipField,
- const double t,
- const ALE::Obj<Mesh>& faultMesh)
+pylith::faults::EqKinSrc::slip(
+ topology::Field<topology::SubMesh>* const slipField,
+ const double t)
{ // slip
assert(0 != _slipfn);
- _slipfn->slip(slipField, t, faultMesh);
+ _slipfn->slip(slipField, t);
} // slip
// ----------------------------------------------------------------------
// Get slip increment on fault surface from time t0 to 1.
void
-pylith::faults::EqKinSrc::slipIncr(const ALE::Obj<pylith::real_section_type>& slipField,
- const double t0,
- const double t1,
- const ALE::Obj<Mesh>& faultMesh)
+pylith::faults::EqKinSrc::slipIncr(
+ topology::Field<topology::SubMesh>* const slipField,
+ const double t0,
+ const double t1)
{ // slip
assert(0 != _slipfn);
- _slipfn->slipIncr(slipField, t0, t1, faultMesh);
+ _slipfn->slipIncr(slipField, t0, t1);
} // slip
// ----------------------------------------------------------------------
// Get final slip.
-ALE::Obj<pylith::real_section_type>
-pylith::faults::EqKinSrc::finalSlip(void)
+const pylith::topology::Field<pylith::topology::SubMesh>&
+pylith::faults::EqKinSrc::finalSlip(void) const
{ // finalSlip
assert(0 != _slipfn);
return _slipfn->finalSlip();
@@ -106,8 +106,8 @@
// ----------------------------------------------------------------------
// Get time when slip begins at each point.
-ALE::Obj<pylith::real_section_type>
-pylith::faults::EqKinSrc::slipTime(void)
+const pylith::topology::Field<pylith::topology::SubMesh>&
+pylith::faults::EqKinSrc::slipTime(void) const
{ // slipTime
assert(0 != _slipfn);
return _slipfn->slipTime();
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/EqKinSrc.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/EqKinSrc.hh 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/EqKinSrc.hh 2009-04-05 03:11:47 UTC (rev 14596)
@@ -22,32 +22,14 @@
#if !defined(pylith_faults_eqkinsrc_hh)
#define pylith_faults_eqkinsrc_hh
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+// Include directives ---------------------------------------------------
+#include "faultsfwd.hh" // forward declarations
-/// Namespace for pylith package
-namespace pylith {
- namespace faults {
- class EqKinSrc;
- class TestEqKinSrc; // unit testing
+#include "pylith/topology/topologyfwd.hh" // USES Field<SubMesh>
- class SlipTimeFn; // HOLDSA SlipTimeFn
- } // faults
-} // pylith
+#include "spatialdata/units/unitsfwd.hh" // USES Nondimensional
-/// Namespace for spatialdata package
-namespace spatialdata {
- namespace spatialdb {
- class SpatialDB;
- } // spatialdb
- namespace geocoords {
- class CoordSys;
- } // geocoord
- namespace units {
- class Nondimensional;
- } // units
-} // spatialdata
-
-/// C++ oject for managing parameters for a kinematic earthquake source.
+// EqKinSrc -------------------------------------------------------------
class pylith::faults::EqKinSrc
{ // class EqKinSrc
friend class TestEqKinSrc; // unit testing
@@ -82,55 +64,47 @@
/** Initialize slip time function.
*
* @param faultMesh Finite-element mesh of fault.
- * @param cs Coordinate system for mesh
* @param normalizer Nondimensionalization of scales.
*/
- void initialize(const ALE::Obj<Mesh>& faultMesh,
- const spatialdata::geocoords::CoordSys* cs,
+ void initialize(const topology::SubMesh& faultMesh,
const spatialdata::units::Nondimensional& normalizer);
/** Get slip on fault surface at time t.
*
* @param slipField Slip field over fault mesh.
* @param t Time t.
- * @param faultMesh Finite-element mesh of fault.
*/
- void slip(const ALE::Obj<real_section_type>& slipField,
- const double t,
- const ALE::Obj<Mesh>& faultMesh);
+ void slip(topology::Field<topology::SubMesh>* const slipField,
+ const double t);
/** Get increment of slip on fault surface between time t0 and t1.
*
- * @param slipField Slip field over fault mesh.
- * @param t Time t.
- * @param faultMesh Finite-element mesh of fault.
+ * @param slipField Slip increment field over fault mesh.
+ * @param t0 Time for start of slip increment.
+ * @param t1 Time for end of slip increment.
*/
- void slipIncr(const ALE::Obj<real_section_type>& slipField,
+ void slipIncr(topology::Field<topology::SubMesh>* const slipField,
const double t0,
- const double t1,
- const ALE::Obj<Mesh>& faultMesh);
+ const double t1);
/** Get final slip.
*
* @returns Final slip.
*/
- ALE::Obj<real_section_type> finalSlip(void);
+ const topology::Field<topology::SubMesh>& finalSlip(void) const;
/** Get time when slip begins at each point.
*
* @returns Time when slip begins.
*/
- ALE::Obj<real_section_type> slipTime(void);
+ const topology::Field<topology::SubMesh>& slipTime(void) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
- /// Not implemented
- EqKinSrc(const EqKinSrc& s);
+ EqKinSrc(const EqKinSrc&); ///< Not implemented
+ const EqKinSrc& operator=(const EqKinSrc&); ///< Not implemented
- /// Not implemented
- const EqKinSrc& operator=(const EqKinSrc& s);
-
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.cc 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.cc 2009-04-05 03:11:47 UTC (rev 14596)
@@ -14,8 +14,9 @@
#include "LiuCosSlipFn.hh" // implementation of object methods
-#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
-#include "pylith/utils/array.hh" // USES double_array
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/Field.hh" // USES Field
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -25,22 +26,19 @@
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
-namespace pylith {
- namespace faults {
- namespace _LiuCosSlipFn {
- const int offsetRiseTime = 0;
- const int offsetSlipTime = 1;
- } // _LiuCosSlipFn
- } // faults
-} // pylith
+// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
+typedef pylith::topology::SubMesh::SieveMesh::label_sequence label_sequence;
+typedef pylith::topology::SubMesh::RealSection RealSection;
// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::LiuCosSlipFn::LiuCosSlipFn(void) :
+ _slipTimeVertex(0),
+ _riseTimeVertex(0),
_dbFinalSlip(0),
_dbSlipTime(0),
- _dbRiseTime(0),
- _spaceDim(0)
+ _dbRiseTime(0)
{ // constructor
} // constructor
@@ -48,50 +46,68 @@
// Destructor.
pylith::faults::LiuCosSlipFn::~LiuCosSlipFn(void)
{ // destructor
- _dbFinalSlip = 0;
- _dbSlipTime = 0;
- _dbRiseTime = 0;
+ delete _parameters; _parameters = 0;
+ _dbFinalSlip = 0; // :TODO: Use shared pointer.
+ _dbSlipTime = 0; // :TODO: Use shared pointer.
+ _dbRiseTime = 0; // :TODO: Use shared pointer.
} // destructor
// ----------------------------------------------------------------------
// Initialize slip time function.
void
pylith::faults::LiuCosSlipFn::initialize(
- const ALE::Obj<Mesh>& faultMesh,
- const spatialdata::geocoords::CoordSys* cs,
- const spatialdata::units::Nondimensional& normalizer,
- const double originTime)
+ const topology::SubMesh& faultMesh,
+ const spatialdata::units::Nondimensional& normalizer,
+ const double originTime)
{ // initialize
- assert(!faultMesh.isNull());
- assert(0 != cs);
assert(0 != _dbFinalSlip);
assert(0 != _dbSlipTime);
assert(0 != _dbRiseTime);
- _spaceDim = cs->spaceDim();
- const int spaceDim = _spaceDim;
- const int indexFinalSlip = 0;
- const int indexRiseTime = spaceDim + _LiuCosSlipFn::offsetRiseTime;
- const int indexSlipTime = spaceDim + _LiuCosSlipFn::offsetSlipTime;
+ const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
+ assert(0 != cs);
+ const int spaceDim = cs->spaceDim();
+ const double lengthScale = normalizer.lengthScale();
+ const double timeScale = normalizer.timeScale();
+
// 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<SieveMesh>& sieveMesh = faultMesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<label_sequence>& vertices = sieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const label_sequence::iterator verticesEnd = vertices->end();
- const int fiberDim = spaceDim + 2;
- _parameters = new real_section_type(faultMesh->comm(), faultMesh->debug());
- _parameters->addSpace(); // final slip
- _parameters->addSpace(); // rise time
- _parameters->addSpace(); // slip time
- assert(3 == _parameters->getNumSpaces());
- _parameters->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(), vertices->end()), *std::max_element(vertices->begin(), vertices->end())+1));
- _parameters->setFiberDimension(vertices, fiberDim);
- _parameters->setFiberDimension(vertices, spaceDim, 0); // final slip
- _parameters->setFiberDimension(vertices, 1, 1); // rise time
- _parameters->setFiberDimension(vertices, 1, 2); // slip time
- faultMesh->allocate(_parameters);
- assert(!_parameters.isNull());
+ delete _parameters; _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(faultMesh);
+ assert(0 != _parameters);
+ _parameters->add("final slip", "final slip");
+ topology::Field<topology::SubMesh>& finalSlip =
+ _parameters->get("final slip");
+ finalSlip.newSection(vertices, spaceDim);
+ finalSlip.allocate();
+ finalSlip.scale(lengthScale);
+ finalSlip.vectorFieldType(topology::FieldBase::VECTOR);
+ const ALE::Obj<RealSection>& finalSlipSection = finalSlip.section();
+ assert(!finalSlipSection.isNull());
+ _parameters->add("slip time", "slip time");
+ topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
+ slipTime.newSection(finalSlipSection->getChart(), 1);
+ slipTime.allocate();
+ slipTime.scale(timeScale);
+ slipTime.vectorFieldType(topology::FieldBase::SCALAR);
+ const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
+ assert(!slipTimeSection.isNull());
+
+ _parameters->add("rise time", "rise time");
+ topology::Field<topology::SubMesh>& riseTime = _parameters->get("rise time");
+ riseTime.newSection(slipTime);
+ riseTime.allocate();
+ riseTime.scale(timeScale);
+ riseTime.vectorFieldType(topology::FieldBase::SCALAR);
+ const ALE::Obj<RealSection>& riseTimeSection = riseTime.section();
+ assert(!riseTimeSection.isNull());
+
// Open databases and set query values
_dbFinalSlip->open();
switch (spaceDim)
@@ -121,69 +137,66 @@
_dbSlipTime->queryVals(slipTimeValues, 1);
_dbRiseTime->open();
- const char* peakRateValues[] = {"rise-time"};
- _dbRiseTime->queryVals(peakRateValues, 1);
+ const char* riseTimeValues[] = {"rise-time"};
+ _dbSlipTime->queryVals(slipTimeValues, 1);
// Get coordinates of vertices
- const ALE::Obj<real_section_type>& coordinates =
- faultMesh->getRealSection("coordinates");
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- const double lengthScale = normalizer.lengthScale();
- const double timeScale = normalizer.timeScale();
-
- double_array paramsVertex(fiberDim);
+ _slipVertex.resize(spaceDim);
double_array vCoordsGlobal(spaceDim);
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ for (label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter) {
coordinates->restrictPoint(*v_iter,
&vCoordsGlobal[0], vCoordsGlobal.size());
normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(),
lengthScale);
-
- int err = _dbFinalSlip->query(¶msVertex[indexFinalSlip], spaceDim,
- &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
+
+ int err = _dbFinalSlip->query(&_slipVertex[0], _slipVertex.size(),
+ &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
- msg << "Could not find final slip at (";
+ msg << "Could not find slip rate at (";
for (int i=0; i < spaceDim; ++i)
msg << " " << vCoordsGlobal[i];
msg << ") using spatial database " << _dbFinalSlip->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(¶msVertex[indexFinalSlip], spaceDim,
+ normalizer.nondimensionalize(&_slipVertex[0], _slipVertex.size(),
lengthScale);
- err = _dbRiseTime->query(¶msVertex[indexRiseTime], 1,
+ err = _dbSlipTime->query(&_slipTimeVertex, 1,
&vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
- msg << "Could not find rise time at (";
+ msg << "Could not find slip initiation time at (";
for (int i=0; i < spaceDim; ++i)
msg << " " << vCoordsGlobal[i];
- msg << ") using spatial database " << _dbRiseTime->label() << ".";
+ msg << ") using spatial database " << _dbSlipTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(¶msVertex[indexRiseTime], spaceDim,
- timeScale);
+ normalizer.nondimensionalize(&_slipTimeVertex, 1, timeScale);
+ // add origin time to rupture time
+ _slipTimeVertex += originTime;
- err = _dbSlipTime->query(¶msVertex[indexSlipTime], 1,
+ err = _dbRiseTime->query(&_riseTimeVertex, 1,
&vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
- msg << "Could not find slip initiation time at (";
+ msg << "Could not find rise time at (";
for (int i=0; i < spaceDim; ++i)
msg << " " << vCoordsGlobal[i];
- msg << ") using spatial database " << _dbSlipTime->label() << ".";
+ msg << ") using spatial database " << _dbRiseTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(¶msVertex[indexSlipTime], 1,
- timeScale);
- // add origin time to rupture time
- paramsVertex[indexSlipTime] += originTime;
+ normalizer.nondimensionalize(&_riseTimeVertex, 1, timeScale);
- _parameters->updatePoint(*v_iter, ¶msVertex[0]);
+ finalSlipSection->updatePoint(*v_iter, &_slipVertex[0]);
+ slipTimeSection->updatePoint(*v_iter, &_slipTimeVertex);
+ riseTimeSection->updatePoint(*v_iter, &_riseTimeVertex);
} // for
// Close databases
@@ -195,122 +208,138 @@
// ----------------------------------------------------------------------
// Get slip on fault surface at time t.
void
-pylith::faults::LiuCosSlipFn::slip(const ALE::Obj<pylith::real_section_type>& slipField,
- const double t,
- const ALE::Obj<Mesh>& faultMesh)
+pylith::faults::LiuCosSlipFn::slip(topology::Field<topology::SubMesh>* slip,
+ const double t)
{ // slip
- assert(!_parameters.isNull());
- assert(!slipField.isNull());
- assert(!faultMesh.isNull());
+ assert(0 != slip);
+ assert(0 != _parameters);
- const int spaceDim = _spaceDim;
- const int indexFinalSlip = 0;
- const int indexRiseTime = spaceDim + _LiuCosSlipFn::offsetRiseTime;
- const int indexSlipTime = spaceDim + _LiuCosSlipFn::offsetSlipTime;
-
- double_array slipValues(spaceDim);
-
// 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 numVertices = vertices->size();
+ const ALE::Obj<SieveMesh>& sieveMesh = slip->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<label_sequence>& vertices = sieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const label_sequence::iterator verticesEnd = vertices->end();
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ // Get sections
+ const topology::Field<topology::SubMesh>& finalSlip =
+ _parameters->get("final slip");
+ const ALE::Obj<RealSection>& finalSlipSection = finalSlip.section();
+ assert(!finalSlipSection.isNull());
+ const topology::Field<topology::SubMesh>& slipTime =
+ _parameters->get("slip time");
+ const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
+ assert(!slipTimeSection.isNull());
+ const topology::Field<topology::SubMesh>& riseTime =
+ _parameters->get("rise time");
+ const ALE::Obj<RealSection>& riseTimeSection = riseTime.section();
+ assert(!riseTimeSection.isNull());
+ const ALE::Obj<RealSection>& slipSection = slip->section();
+ assert(!slipSection.isNull());
+
+ const int spaceDim = _slipVertex.size();
+ for (label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter) {
- const real_section_type::value_type* paramsVertex =
- _parameters->restrictPoint(*v_iter);
- assert(0 != paramsVertex);
+ finalSlipSection->restrictPoint(*v_iter, &_slipVertex[0],
+ _slipVertex.size());
+ slipTimeSection->restrictPoint(*v_iter, &_slipTimeVertex, 1);
+ riseTimeSection->restrictPoint(*v_iter, &_riseTimeVertex, 1);
- const double* finalSlip = ¶msVertex[indexFinalSlip];
- const double riseTime = paramsVertex[indexRiseTime];
- const double slipTime = paramsVertex[indexSlipTime];
-
double finalSlipMag = 0.0;
for (int i=0; i < spaceDim; ++i)
- finalSlipMag += finalSlip[i]*finalSlip[i];
+ finalSlipMag += _slipVertex[i]*_slipVertex[i];
finalSlipMag = sqrt(finalSlipMag);
- const double slip = _slipFn(t-slipTime, finalSlipMag, riseTime);
+ const double slip = _slipFn(t-_slipTimeVertex, finalSlipMag,
+ _riseTimeVertex);
const double scale = finalSlipMag > 0.0 ? slip / finalSlipMag : 0.0;
- for (int i=0; i < spaceDim; ++i)
- slipValues[i] = scale * finalSlip[i];
-
+ _slipVertex *= scale;
+
// Update field
- slipField->updateAddPoint(*v_iter, &slipValues[0]);
+ slipSection->updateAddPoint(*v_iter, &_slipVertex[0]);
} // for
- PetscLogFlops(numVertices * (2+8 + 3*spaceDim));
+ PetscLogFlops(vertices->size() * (2+28 + 3*_slipVertex.size()));
} // slip
// ----------------------------------------------------------------------
// Get increment of slip on fault surface between time t0 and t1.
void
-pylith::faults::LiuCosSlipFn::slipIncr(const ALE::Obj<pylith::real_section_type>& slipField,
+pylith::faults::LiuCosSlipFn::slipIncr( topology::Field<topology::SubMesh>* slip,
const double t0,
- const double t1,
- const ALE::Obj<Mesh>& faultMesh)
+ const double t1)
{ // slipIncr
- assert(!_parameters.isNull());
- assert(!slipField.isNull());
- assert(!faultMesh.isNull());
+ assert(0 != slip);
+ assert(0 != _parameters);
- const int spaceDim = _spaceDim;
- const int indexFinalSlip = 0;
- const int indexRiseTime = spaceDim + _LiuCosSlipFn::offsetRiseTime;
- const int indexSlipTime = spaceDim + _LiuCosSlipFn::offsetSlipTime;
-
- double_array slipValues(spaceDim);
-
// 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 numVertices = vertices->size();
+ const ALE::Obj<SieveMesh>& sieveMesh = slip->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<label_sequence>& vertices = sieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const label_sequence::iterator verticesEnd = vertices->end();
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ // Get sections
+ const topology::Field<topology::SubMesh>& finalSlip =
+ _parameters->get("final slip");
+ const ALE::Obj<RealSection>& finalSlipSection = finalSlip.section();
+ assert(!finalSlipSection.isNull());
+ const topology::Field<topology::SubMesh>& slipTime =
+ _parameters->get("slip time");
+ const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
+ assert(!slipTimeSection.isNull());
+ const topology::Field<topology::SubMesh>& riseTime =
+ _parameters->get("rise time");
+ const ALE::Obj<RealSection>& riseTimeSection = riseTime.section();
+ assert(!riseTimeSection.isNull());
+ const ALE::Obj<RealSection>& slipSection = slip->section();
+ assert(!slipSection.isNull());
+
+ const int spaceDim = _slipVertex.size();
+ for (label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter) {
- const real_section_type::value_type* paramsVertex =
- _parameters->restrictPoint(*v_iter);
- assert(0 != paramsVertex);
+ finalSlipSection->restrictPoint(*v_iter, &_slipVertex[0],
+ _slipVertex.size());
+ slipTimeSection->restrictPoint(*v_iter, &_slipTimeVertex, 1);
+ riseTimeSection->restrictPoint(*v_iter, &_riseTimeVertex, 1);
- const double* finalSlip = ¶msVertex[indexFinalSlip];
- const double riseTime = paramsVertex[indexRiseTime];
- const double slipTime = paramsVertex[indexSlipTime];
-
double finalSlipMag = 0.0;
for (int i=0; i < spaceDim; ++i)
- finalSlipMag += finalSlip[i]*finalSlip[i];
+ finalSlipMag += _slipVertex[i]*_slipVertex[i];
finalSlipMag = sqrt(finalSlipMag);
- const double slip0 = _slipFn(t0-slipTime, finalSlipMag, riseTime);
- const double slip1 = _slipFn(t1-slipTime, finalSlipMag, riseTime);
+ const double slip0 = _slipFn(t0-_slipTimeVertex, finalSlipMag,
+ _riseTimeVertex);
+ const double slip1 = _slipFn(t1-_slipTimeVertex, finalSlipMag,
+ _riseTimeVertex);
const double scale = finalSlipMag > 0.0 ?
(slip1 - slip0) / finalSlipMag : 0.0;
- for (int i=0; i < spaceDim; ++i)
- slipValues[i] = scale * finalSlip[i];
+ _slipVertex *= scale;
+
// Update field
- slipField->updateAddPoint(*v_iter, &slipValues[0]);
+ slipSection->updateAddPoint(*v_iter, &_slipVertex[0]);
} // for
- PetscLogFlops(numVertices * (3+2*8 + 3*spaceDim));
+ PetscLogFlops(vertices->size() * (3+2*28 + 3*_slipVertex.size()));
} // slipIncr
// ----------------------------------------------------------------------
// Get final slip.
-ALE::Obj<pylith::real_section_type>
+const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::LiuCosSlipFn::finalSlip(void)
{ // finalSlip
- return _parameters->getFibration(0);
+ return _parameters->get("final slip");
} // finalSlip
// ----------------------------------------------------------------------
// Get time when slip begins at each point.
-ALE::Obj<pylith::real_section_type>
+const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::LiuCosSlipFn::slipTime(void)
{ // slipTime
- return _parameters->getFibration(2);
+ return _parameters->get("slip time");
} // slipTime
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.hh 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.hh 2009-04-05 03:11:47 UTC (rev 14596)
@@ -22,27 +22,14 @@
#if !defined(pylith_faults_liucosslipfn_hh)
#define pylith_faults_liucosslipfn_hh
+// Include directives ---------------------------------------------------
#include "SlipTimeFn.hh"
-/// Namespace for pylith package
-namespace pylith {
- namespace faults {
- class LiuCosSlipFn;
- class TestLiuCosSlipFn; // unit testing
- } // faults
-} // pylith
+#include "pylith/topology/topologyfwd.hh" // USES Fields<Field<SubMesh> >
-/// Namespace for spatialdata package
-namespace spatialdata {
- namespace spatialdb {
- class SpatialDB;
- } // spatialdb
- namespace units {
- class Nondimensional;
- } // units
-} // spatialdata
+#include "pylith/utils/array.hh" // HASA double_array
-/// C++ implementation of LiuCos slip time function.
+// LiuCosSlipFn ---------------------------------------------------------
class pylith::faults::LiuCosSlipFn : public SlipTimeFn
{ // class LiuCosSlipFn
friend class TestLiuCosSlipFn; // unit testing
@@ -54,7 +41,6 @@
LiuCosSlipFn(void);
/// Destructor.
- virtual
~LiuCosSlipFn(void);
/** Set spatial database for final slip.
@@ -80,12 +66,11 @@
/** Initialize slip time function.
*
* @param faultMesh Finite-element mesh of fault.
- * @param cs Coordinate system for mesh.
- * @param originTime Origin time for earthquake source.
+ * @param cs Coordinate system for mesh
* @param normalizer Nondimensionalization of scales.
+ * @param originTime Origin time for earthquake source.
*/
- void initialize(const ALE::Obj<Mesh>& faultMesh,
- const spatialdata::geocoords::CoordSys* cs,
+ void initialize(const topology::SubMesh& faultMesh,
const spatialdata::units::Nondimensional& normalizer,
const double originTime =0.0);
@@ -93,50 +78,42 @@
*
* @param slipField Slip field over fault surface.
* @param t Time t.
- * @param faultMesh Mesh over fault surface.
*
* @returns Slip vector as left-lateral/reverse/normal.
*/
- void slip(const ALE::Obj<real_section_type>& slipField,
- const double t,
- const ALE::Obj<Mesh>& faultMesh);
+ void slip(topology::Field<topology::SubMesh>* const slipField,
+ const double t);
/** Get slip increment on fault surface between time t0 and t1.
*
* @param slipField Slip field over fault surface.
* @param t0 Time t.
* @param t1 Time t+dt.
- * @param faultMesh Mesh over fault surface.
*
* @returns Increment in slip vector as left-lateral/reverse/normal.
*/
- void slipIncr(const ALE::Obj<real_section_type>& slipField,
+ void slipIncr(topology::Field<topology::SubMesh>* slipField,
const double t0,
- const double t1,
- const ALE::Obj<Mesh>& faultMesh);
+ const double t1);
-
/** Get final slip.
*
* @returns Final slip.
*/
- ALE::Obj<real_section_type> finalSlip(void);
+ const topology::Field<topology::SubMesh>& finalSlip(void);
/** Get time when slip begins at each point.
*
* @returns Time when slip begins.
*/
- ALE::Obj<real_section_type> slipTime(void);
+ const topology::Field<topology::SubMesh>& slipTime(void);
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
- /// Not implemented
- LiuCosSlipFn(const LiuCosSlipFn& m);
+ LiuCosSlipFn(const LiuCosSlipFn&); ///< Not implemented
+ const LiuCosSlipFn& operator=(const LiuCosSlipFn&); ///< Not implemented
- /// Not implemented
- const LiuCosSlipFn& operator=(const LiuCosSlipFn& f);
-
// PRIVATE METHODS //////////////////////////////////////////////////////
private :
@@ -156,10 +133,14 @@
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
- /// Parameters for LiuCos slip time function.
- /// Final slip (vector), peak slip rate (scalar), slip time (scalar).
- ALE::Obj<real_section_type> _parameters;
+ double _slipTimeVertex; ///< Slip time at a vertex.
+ double _riseTimeVertex; ///< Rise time at a vertex.
+ double_array _slipVertex; ///< Slip at a vertex.
+ /// Parameters for Liu cosine/sine slip time function, final slip
+ /// (vector), slip time (scalar), rise time (scalar).
+ topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
+
/// Spatial database for final slip.
spatialdata::spatialdb::SpatialDB* _dbFinalSlip;
@@ -169,8 +150,6 @@
/// Spatial database for rise time.
spatialdata::spatialdb::SpatialDB* _dbRiseTime;
- int _spaceDim; ///< Spatial dimension for slip field.
-
}; // class LiuCosSlipFn
#include "LiuCosSlipFn.icc" // inline methods
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/SlipTimeFn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/SlipTimeFn.hh 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/SlipTimeFn.hh 2009-04-05 03:11:47 UTC (rev 14596)
@@ -21,9 +21,12 @@
#define pylith_faults_sliptimefn_hh
// Include directives ---------------------------------------------------
-#include "Fault.hh" // ISA Fault
+#include "faultsfwd.hh" // forward declarations
+#include "pylith/topology/topologyfwd.hh" // USES Field<SubMesh>
+
#include "spatialdata/units/unitsfwd.hh" // USES Nondimensional
+#include "spatialdata/spatialdb/spatialdbfwd.hh" // USES SpatialDB
// SlipTimeFn -----------------------------------------------------------
class pylith::faults::SlipTimeFn
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.hh 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.hh 2009-04-05 03:11:47 UTC (rev 14596)
@@ -29,7 +29,7 @@
#include "pylith/utils/array.hh" // HASA double_array
-// SlipTimeFn -----------------------------------------------------------
+// StepSlipFn -----------------------------------------------------------
class pylith::faults::StepSlipFn : public SlipTimeFn
{ // class StepSlipFn
friend class TestStepSlipFn; // unit testing
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/faultsfwd.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/faultsfwd.hh 2009-04-04 23:31:36 UTC (rev 14595)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/faultsfwd.hh 2009-04-05 03:11:47 UTC (rev 14596)
@@ -30,9 +30,9 @@
class FaultCohesiveDyn;
class FaultCohesiveKin;
- class EqSinSrc;
+ class EqKinSrc;
class SlipTimeFn;
- class BrundSlipFn;
+ class BruneSlipFn;
class ConstRateSlipFn;
class LiuCosSlipFn;
class StepSlipFn;
More information about the CIG-COMMITS
mailing list