[cig-commits] r13292 - in short/3D/PyLith/trunk: libsrc libsrc/faults modulesrc/faults pylith pylith/faults unittests/libtests/faults unittests/libtests/faults/data
brad at geodynamics.org
brad at geodynamics.org
Tue Nov 11 14:36:31 PST 2008
Author: brad
Date: 2008-11-11 14:36:30 -0800 (Tue, 11 Nov 2008)
New Revision: 13292
Added:
short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc
short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh
short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.icc
short/3D/PyLith/trunk/pylith/faults/LiuCosSlipFn.py
short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.hh
short/3D/PyLith/trunk/unittests/libtests/faults/data/line2_risetime.spatialdb
short/3D/PyLith/trunk/unittests/libtests/faults/data/tet4_risetime.spatialdb
short/3D/PyLith/trunk/unittests/libtests/faults/data/tri3_risetime.spatialdb
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/faults/Makefile.am
short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src
short/3D/PyLith/trunk/pylith/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/faults/data/Makefile.am
Log:
Added Liu, Archuleta, and Hartzell slip time function.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2008-11-11 22:24:06 UTC (rev 13291)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2008-11-11 22:36:30 UTC (rev 13292)
@@ -35,6 +35,7 @@
faults/FaultCohesive.cc \
faults/FaultCohesiveKin.cc \
faults/FaultCohesiveDyn.cc \
+ faults/LiuCosSlipFn.cc \
faults/SlipTimeFn.cc \
faults/StepSlipFn.cc \
feassemble/CellGeometry.cc \
Added: short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc 2008-11-11 22:36:30 UTC (rev 13292)
@@ -0,0 +1,306 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "LiuCosSlipFn.hh" // implementation of object methods
+
+#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
+#include "pylith/utils/array.hh" // USES double_array
+
+#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+
+#include <assert.h> // USES assert()
+#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
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::faults::LiuCosSlipFn::LiuCosSlipFn(void) :
+ _dbFinalSlip(0),
+ _dbSlipTime(0),
+ _dbRiseTime(0),
+ _spaceDim(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::faults::LiuCosSlipFn::~LiuCosSlipFn(void)
+{ // destructor
+ _dbFinalSlip = 0;
+ _dbSlipTime = 0;
+ _dbRiseTime = 0;
+} // destructor
+
+// ----------------------------------------------------------------------
+// Initialize slip time function.
+void
+pylith::faults::LiuCosSlipFn::initialize(
+ const ALE::Obj<Mesh>& faultMesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ 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;
+
+ // 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 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());
+
+ // Open databases and set query values
+ _dbFinalSlip->open();
+ switch (spaceDim)
+ { // switch
+ case 1 : {
+ const char* slipValues[] = {"fault-opening"};
+ _dbFinalSlip->queryVals(slipValues, 1);
+ break;
+ } // case 1
+ case 2 : {
+ const char* slipValues[] = {"left-lateral-slip", "fault-opening"};
+ _dbFinalSlip->queryVals(slipValues, 2);
+ break;
+ } // case 2
+ case 3 : {
+ const char* slipValues[] = {"left-lateral-slip", "reverse-slip",
+ "fault-opening"};
+ _dbFinalSlip->queryVals(slipValues, 3);
+ break;
+ } // case 3
+ default :
+ assert(0);
+ } // switch
+
+ _dbSlipTime->open();
+ const char* slipTimeValues[] = {"slip-time"};
+ _dbSlipTime->queryVals(slipTimeValues, 1);
+
+ _dbRiseTime->open();
+ const char* peakRateValues[] = {"rise-time"};
+ _dbRiseTime->queryVals(peakRateValues, 1);
+
+ // Get coordinates of vertices
+ const ALE::Obj<real_section_type>& coordinates =
+ faultMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
+ 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* coordsVertex =
+ coordinates->restrictPoint(*v_iter);
+ assert(0 != coordsVertex);
+
+ 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 << " " << coordsVertex[i];
+ msg << ") using spatial database " << _dbFinalSlip->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ err = _dbRiseTime->query(¶msVertex[indexRiseTime], 1,
+ coordsVertex, spaceDim, cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find rise time at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << coordsVertex[i];
+ msg << ") using spatial database " << _dbRiseTime->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ err = _dbSlipTime->query(¶msVertex[indexSlipTime], 1,
+ coordsVertex, spaceDim, cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find slip initiation time at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << coordsVertex[i];
+ msg << ") using spatial database " << _dbSlipTime->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ // add origin time to rupture time
+ paramsVertex[indexSlipTime] += originTime;
+
+ _parameters->updatePoint(*v_iter, ¶msVertex[0]);
+ } // for
+
+ // Close databases
+ _dbFinalSlip->close();
+ _dbSlipTime->close();
+ _dbRiseTime->close();
+} // initialize
+
+// ----------------------------------------------------------------------
+// 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)
+{ // slip
+ assert(!_parameters.isNull());
+ assert(!slipField.isNull());
+ assert(!faultMesh.isNull());
+
+ 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();
+
+ for (Mesh::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);
+
+ 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 = sqrt(finalSlipMag);
+
+ const double slip = _slipFn(t-slipTime, finalSlipMag, riseTime);
+ const double scale = finalSlipMag > 0.0 ? slip / finalSlipMag : 0.0;
+ for (int i=0; i < spaceDim; ++i)
+ slipValues[i] = scale * finalSlip[i];
+
+ // Update field
+ slipField->updateAddPoint(*v_iter, &slipValues[0]);
+ } // for
+
+ PetscLogFlops(numVertices * (2+8 + 3*spaceDim));
+} // 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,
+ const double t0,
+ const double t1,
+ const ALE::Obj<Mesh>& faultMesh)
+{ // slipIncr
+ assert(!_parameters.isNull());
+ assert(!slipField.isNull());
+ assert(!faultMesh.isNull());
+
+ 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();
+
+ for (Mesh::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);
+
+ 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 = sqrt(finalSlipMag);
+
+ const double slip0 = _slipFn(t0-slipTime, finalSlipMag, riseTime);
+ const double slip1 = _slipFn(t1-slipTime, finalSlipMag, riseTime);
+ const double scale = finalSlipMag > 0.0 ?
+ (slip1 - slip0) / finalSlipMag : 0.0;
+ for (int i=0; i < spaceDim; ++i)
+ slipValues[i] = scale * finalSlip[i];
+
+ // Update field
+ slipField->updateAddPoint(*v_iter, &slipValues[0]);
+ } // for
+
+ PetscLogFlops(numVertices * (3+2*8 + 3*spaceDim));
+} // slipIncr
+
+// ----------------------------------------------------------------------
+// Get final slip.
+ALE::Obj<pylith::real_section_type>
+pylith::faults::LiuCosSlipFn::finalSlip(void)
+{ // finalSlip
+ return _parameters->getFibration(0);
+} // finalSlip
+
+// ----------------------------------------------------------------------
+// Get time when slip begins at each point.
+ALE::Obj<pylith::real_section_type>
+pylith::faults::LiuCosSlipFn::slipTime(void)
+{ // slipTime
+ return _parameters->getFibration(2);
+} // slipTime
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh 2008-11-11 22:36:30 UTC (rev 13292)
@@ -0,0 +1,176 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/faults/LiuCosSlipFn.hh
+ *
+ * @brief C++ implementation of LiuCos slip time function.
+ *
+ * Sine/cosine slip time function from Liu, Archuleta, and Hartzell,
+ * BSSA, 2006 (doi:10.1785/0120060036) which has a rapid rise and then
+ * a gradual falloff with a finite duration.
+ */
+
+#if !defined(pylith_faults_liucosslipfn_hh)
+#define pylith_faults_liucosslipfn_hh
+
+#include "SlipTimeFn.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace faults {
+ class LiuCosSlipFn;
+ class TestLiuCosSlipFn; // unit testing
+ } // faults
+} // pylith
+
+/// Namespace for spatialdata package
+namespace spatialdata {
+ namespace spatialdb {
+ class SpatialDB;
+ } // spatialdb
+} // spatialdata
+
+/// C++ implementation of LiuCos slip time function.
+class pylith::faults::LiuCosSlipFn : public SlipTimeFn
+{ // class LiuCosSlipFn
+ friend class TestLiuCosSlipFn; // unit testing
+
+// PUBLIC METHODS ///////////////////////////////////////////////////////
+public :
+
+ /// Default constructor.
+ LiuCosSlipFn(void);
+
+ /// Destructor.
+ virtual
+ ~LiuCosSlipFn(void);
+
+ /** Set spatial database for final slip.
+ *
+ * @param db Spatial database
+ */
+ void dbFinalSlip(spatialdata::spatialdb::SpatialDB* const db);
+
+ /** Set spatial database for slip initiation time.
+ *
+ * @param db Spatial database
+ */
+ void dbSlipTime(spatialdata::spatialdb::SpatialDB* const db);
+
+ /** Set spatial database for rise time. The rise time is the time it
+ * takes for the slip to increase from 0.0 to 0.95 of the final
+ * value.
+ *
+ * @param db Spatial database
+ */
+ void dbRiseTime(spatialdata::spatialdb::SpatialDB* const db);
+
+ /** Initialize slip time function.
+ *
+ * @param faultMesh Finite-element mesh of fault.
+ * @param cs Coordinate system for mesh.
+ * @param originTime Origin time for earthquake source.
+ */
+ void initialize(const ALE::Obj<Mesh>& faultMesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ const double originTime =0.0);
+
+ /** Get slip on fault surface at time t.
+ *
+ * @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);
+
+ /** 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,
+ const double t0,
+ const double t1,
+ const ALE::Obj<Mesh>& faultMesh);
+
+
+ /** Get final slip.
+ *
+ * @returns Final slip.
+ */
+ ALE::Obj<real_section_type> finalSlip(void);
+
+ /** Get time when slip begins at each point.
+ *
+ * @returns Time when slip begins.
+ */
+ ALE::Obj<real_section_type> slipTime(void);
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ LiuCosSlipFn(const LiuCosSlipFn& m);
+
+ /// Not implemented
+ const LiuCosSlipFn& operator=(const LiuCosSlipFn& f);
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /** Compute slip using slip time function.
+ *
+ * @param t Time relative to slip starting time at point
+ * @param finalSlip Final slip at point
+ * @param riseTime Rise time (t95) at point
+ *
+ * @returns Slip at point at time t
+ */
+ static
+ double _slipFn(const double t,
+ const double finalSlip,
+ const double riseTime);
+
+// 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;
+
+ /// Spatial database for final slip.
+ spatialdata::spatialdb::SpatialDB* _dbFinalSlip;
+
+ /// Spatial database for slip time.
+ spatialdata::spatialdb::SpatialDB* _dbSlipTime;
+
+ /// Spatial database for rise time.
+ spatialdata::spatialdb::SpatialDB* _dbRiseTime;
+
+ int _spaceDim; ///< Spatial dimension for slip field.
+
+}; // class LiuCosSlipFn
+
+#include "LiuCosSlipFn.icc" // inline methods
+
+#endif // pylith_faults_liucosslipfn_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.icc 2008-11-11 22:36:30 UTC (rev 13292)
@@ -0,0 +1,76 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_faults_liucosslipfn_hh)
+#error "LiuCosSlipFn.icc can only be included from LiuCosSlipFn.hh"
+#endif
+
+#include <math.h> // USES sin(), cos(), M_PI
+
+// Set spatial database for final slip.
+inline
+void
+pylith::faults::LiuCosSlipFn::dbFinalSlip(spatialdata::spatialdb::SpatialDB* const db) {
+ _dbFinalSlip = db;
+} // dbFinalSlip
+
+// Set spatial database for slip initiation time.
+inline
+void
+pylith::faults::LiuCosSlipFn::dbSlipTime(spatialdata::spatialdb::SpatialDB* const db) {
+ _dbSlipTime = db;
+} // dbSlipTime
+
+// Set spatial database for rise time.
+inline
+void
+pylith::faults::LiuCosSlipFn::dbRiseTime(spatialdata::spatialdb::SpatialDB* const db) {
+ _dbRiseTime = db;
+} // dbRiseTime
+
+// Compute slip using slip time function.
+inline
+double
+pylith::faults::LiuCosSlipFn::_slipFn(const double t,
+ const double finalSlip,
+ const double riseTime) {
+ double slip = 0.0;
+ if (t > 0.0) {
+ const double tau = riseTime * 1.525;
+ const double tau1 = 0.13 * tau;
+ const double tau2 = tau - tau1;
+ const double Cn =
+ M_PI / (1.4 * M_PI * tau1 + 1.2 * tau1 + 0.3 * M_PI * tau2);
+
+ if (t <= tau1) {
+ slip = 0.7*t - 0.7*tau1/M_PI*sin(M_PI*t/tau1)
+ - 0.6*tau1/(0.5*M_PI)*(cos(0.5*M_PI*t/tau1) - 1.0);
+ slip *= Cn;
+ } else if (t <= 2.0*tau1) {
+ slip = 1.0*t - 0.7*tau1/M_PI*sin(M_PI*t/tau1) +
+ 0.3*tau2/M_PI*sin(M_PI*(t-tau1)/tau2) +
+ 1.2*tau1/M_PI - 0.3*tau1;
+ slip *= Cn;
+ } else if (t <= tau) {
+ slip = 0.3*t + 0.3*tau2/M_PI*sin(M_PI*(t-tau1)/tau2) +
+ 1.1*tau1 + 1.2*tau1/M_PI;
+ slip *= Cn;
+ } else
+ slip = 1.0;
+ } // if
+ slip *= finalSlip;
+
+ return slip;
+} // _slip
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/Makefile.am 2008-11-11 22:24:06 UTC (rev 13291)
+++ short/3D/PyLith/trunk/libsrc/faults/Makefile.am 2008-11-11 22:36:30 UTC (rev 13292)
@@ -27,6 +27,8 @@
FaultCohesiveDyn.icc \
FaultCohesiveKin.hh \
FaultCohesiveKin.icc \
+ LiuCosSlipFn.hh \
+ LiuCosSlipFn.icc \
SlipTimeFn.hh \
StepSlipFn.hh \
StepSlipFn.icc
Modified: short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src 2008-11-11 22:24:06 UTC (rev 13291)
+++ short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src 2008-11-11 22:36:30 UTC (rev 13292)
@@ -17,6 +17,7 @@
#include "pylith/faults/EqKinSrc.hh"
#include "pylith/faults/SlipTimeFn.hh"
#include "pylith/faults/BruneSlipFn.hh"
+#include "pylith/faults/LiuCosSlipFn.hh"
#include "pylith/faults/ConstRateSlipFn.hh"
#include "pylith/faults/StepSlipFn.hh"
@@ -1171,6 +1172,131 @@
# ----------------------------------------------------------------------
+cdef class LiuCosSlipFn(SlipTimeFn):
+
+ def __init__(self):
+ """
+ Constructor.
+ """
+ # create shim for constructor
+ #embed{ void* LiuCosSlipFn_constructor()
+ void* result = 0;
+ try {
+ result = (void*)(new pylith::faults::LiuCosSlipFn);
+ assert(0 != result);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ return result;
+ #}embed
+
+ SlipTimeFn.__init__(self)
+ self.thisptr = LiuCosSlipFn_constructor()
+ self.handle = self._createHandle()
+ return
+
+
+ property dbFinalSlip:
+ def __set__(self, value):
+ """
+ Set database for final slip.
+ """
+ # create shim for method 'dbFinalSlip'
+ #embed{ void LiuCosSlipFn_dbFinalSlip_set(void* objVptr, void* dbVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != dbVptr);
+ spatialdata::spatialdb::SpatialDB* db =
+ (spatialdata::spatialdb::SpatialDB*) dbVptr;
+ ((pylith::faults::LiuCosSlipFn*) objVptr)->dbFinalSlip(db);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if not value.name == "spatialdata_spatialdb_SpatialDB":
+ raise TypeError, \
+ "Argument must be extension module type " \
+ "'spatialdata::spatialdb::SpatialDB'."
+ LiuCosSlipFn_dbFinalSlip_set(self.thisptr, ptrFromHandle(value))
+
+
+ property dbSlipTime:
+ def __set__(self, value):
+ """
+ Set database for slip initiation time.
+ """
+ # create shim for method 'dbSlipTime'
+ #embed{ void LiuCosSlipFn_dbSlipTime_set(void* objVptr, void* dbVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != dbVptr);
+ spatialdata::spatialdb::SpatialDB* db =
+ (spatialdata::spatialdb::SpatialDB*) dbVptr;
+ ((pylith::faults::LiuCosSlipFn*) objVptr)->dbSlipTime(db);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if not value.name == "spatialdata_spatialdb_SpatialDB":
+ raise TypeError, \
+ "Argument must be extension module type " \
+ "'spatialdata::spatialdb::SpatialDB'."
+ LiuCosSlipFn_dbSlipTime_set(self.thisptr, ptrFromHandle(value))
+
+
+ property dbRiseTime:
+ def __set__(self, value):
+ """
+ Set database for rise time.
+ """
+ # create shim for method 'dbRiseTime'
+ #embed{ void LiuCosSlipFn_dbRiseTime_set(void* objVptr, void* dbVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != dbVptr);
+ spatialdata::spatialdb::SpatialDB* db =
+ (spatialdata::spatialdb::SpatialDB*) dbVptr;
+ ((pylith::faults::LiuCosSlipFn*) objVptr)->dbRiseTime(db);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if not value.name == "spatialdata_spatialdb_SpatialDB":
+ raise TypeError, \
+ "Argument must be extension module type " \
+ "'spatialdata::spatialdb::SpatialDB'."
+ LiuCosSlipFn_dbRiseTime_set(self.thisptr, ptrFromHandle(value))
+
+
+# ----------------------------------------------------------------------
cdef class ConstRateSlipFn(SlipTimeFn):
def __init__(self):
Modified: short/3D/PyLith/trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/pylith/Makefile.am 2008-11-11 22:24:06 UTC (rev 13291)
+++ short/3D/PyLith/trunk/pylith/Makefile.am 2008-11-11 22:36:30 UTC (rev 13292)
@@ -27,6 +27,7 @@
faults/Fault.py \
faults/FaultCohesive.py \
faults/FaultCohesiveKin.py \
+ faults/LiuCosSlipFn.py \
faults/SingleRupture.py \
faults/SlipTimeFn.py \
faults/StepSlipFn.py \
Added: short/3D/PyLith/trunk/pylith/faults/LiuCosSlipFn.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/LiuCosSlipFn.py (rev 0)
+++ short/3D/PyLith/trunk/pylith/faults/LiuCosSlipFn.py 2008-11-11 22:36:30 UTC (rev 13292)
@@ -0,0 +1,131 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/faults/BruneSlipFn.py
+##
+## @brief Sine/cosine slip time function from Liu, Archuleta, and Hartzell,
+## BSSA, 2006 (doi:10.1785/0120060036) which has a rapid rise and then
+## a gradual falloff with a finite duration.
+##
+## Factory: slip_time_fn
+
+from SlipTimeFn import SlipTimeFn
+
+# LiuCosSlipFn class
+class LiuCosSlipFn(SlipTimeFn):
+ """
+ Sine/cosine slip time function from Liu, Archuleta, and Hartzell,
+ BSSA, 2006 (doi:10.1785/0120060036) which has a rapid rise and then
+ a gradual falloff with a finite duration.
+
+ Factory: slip_time_fn
+ """
+
+ # INVENTORY //////////////////////////////////////////////////////////
+
+ class Inventory(SlipTimeFn.Inventory):
+ """
+ Python object for managing LiuCosSlipFn facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing LiuCosSlipFn facilities and properties.
+ ##
+ ## \b Properties
+ ## @li None
+ ##
+ ## \b Facilities
+ ## @li \b slip Spatial database of final slip.
+ ## @li \b slip_time Spatial database of slip initiation time.
+ ## @li \b rise_time Spatial database of rise time (t95).
+
+ import pyre.inventory
+
+ from spatialdata.spatialdb.SimpleDB import SimpleDB
+
+ slip = pyre.inventory.facility("slip", family="spatial_database",
+ factory=SimpleDB)
+ slip.meta['tip'] = "Spatial database of slip."
+
+ slipTime = pyre.inventory.facility("slip_time", family="spatial_database",
+ factory=SimpleDB)
+ slipTime.meta['tip'] = "Spatial database of slip initiation time."
+
+ riseTime = pyre.inventory.facility("rise_time", family="spatial_database",
+ factory=SimpleDB)
+ riseTime.meta['tip'] = "Spatial database of rise time (t95)."
+
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="bruneslipfn"):
+ """
+ Constructor.
+ """
+ SlipTimeFn.__init__(self, name)
+ self._loggingPrefix = "BrSF "
+ return
+
+
+ def initialize(self):
+ """
+ Initialize.
+ """
+ logEvent = "%sinit" % self._loggingPrefix
+ self._logger.eventBegin(logEvent)
+
+ self.slip.initialize()
+ self.slipTime.initialize()
+ self.riseTime.initialize()
+ assert(None != self.cppHandle)
+
+ self.cppHandle.dbFinalSlip = self.slip.cppHandle
+ self.cppHandle.dbSlipTime = self.slipTime.cppHandle
+ self.cppHandle.dbRiseTime = self.riseTime.cppHandle
+
+ self._logger.eventEnd(logEvent)
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Setup members using inventory.
+ """
+ SlipTimeFn._configure(self)
+ self.slip = self.inventory.slip
+ self.slipTime = self.inventory.slipTime
+ self.riseTime = self.inventory.riseTime
+ return
+
+
+ def _createCppHandle(self):
+ """
+ Create handle to C++ object.
+ """
+ if None == self.cppHandle:
+ import pylith.faults.faults as bindings
+ self.cppHandle = bindings.LiuCosSlipFn()
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def slip_time_fn():
+ """
+ Factory associated with LiuCosSlipFn.
+ """
+ return LiuCosSlipFn()
+
+
+# End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am 2008-11-11 22:24:06 UTC (rev 13291)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am 2008-11-11 22:36:30 UTC (rev 13292)
@@ -22,6 +22,7 @@
# Primary source files
testfaults_SOURCES = \
TestBruneSlipFn.cc \
+ TestLiuCosSlipFn.cc \
TestConstRateSlipFn.cc \
TestStepSlipFn.cc \
TestEqKinSrc.cc \
@@ -50,6 +51,7 @@
noinst_HEADERS = \
TestBruneSlipFn.hh \
+ TestLiuCosSlipFn.hh \
TestConstRateSlipFn.hh \
TestStepSlipFn.hh \
TestEqKinSrc.hh \
Added: short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.cc 2008-11-11 22:36:30 UTC (rev 13292)
@@ -0,0 +1,546 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestLiuCosSlipFn.hh" // Implementation of class methods
+
+#include "pylith/faults/LiuCosSlipFn.hh" // USES LiuCosSlipFn
+
+#include "pylith/faults/CohesiveTopology.hh" // USES CohesiveTopology
+#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
+#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
+
+#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
+#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestLiuCosSlipFn );
+
+// ----------------------------------------------------------------------
+namespace pylith {
+ namespace faults {
+ namespace _TestLiuCosSlipFn {
+ struct DataStruct {
+ const char* meshFilename;
+ const char* faultLabel;
+ const int faultId;
+ const char* finalSlipFilename;
+ const char* slipTimeFilename;
+ const char* riseTimeFilename;
+ const int* constraintPts;
+ const double* finalSlipE;
+ const double* slipTimeE;
+ const double* riseTimeE;
+ const int numConstraintPts;
+ }; // DataStruct
+ } // _TestLiuCosSlipFn
+ } // faults
+} // pylith
+
+// ----------------------------------------------------------------------
+// Test constructor.
+void
+pylith::faults::TestLiuCosSlipFn::testConstructor(void)
+{ // testConstructor
+ LiuCosSlipFn slipfn;
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test dbFinalSlip().
+void
+pylith::faults::TestLiuCosSlipFn::testDbFinalSlip(void)
+{ // testDbFinalSlip
+ const char* label = "database ABC";
+ LiuCosSlipFn slipfn;
+
+ spatialdata::spatialdb::SimpleDB db(label);
+ slipfn.dbFinalSlip(&db);
+
+ CPPUNIT_ASSERT(0 != slipfn._dbFinalSlip);
+ CPPUNIT_ASSERT_EQUAL(std::string(label),
+ std::string(slipfn._dbFinalSlip->label()));
+ CPPUNIT_ASSERT(0 == slipfn._dbSlipTime);
+ CPPUNIT_ASSERT(0 == slipfn._dbRiseTime);
+} // testDbFinalSlip
+
+// ----------------------------------------------------------------------
+// Test dbSlipTime().
+void
+pylith::faults::TestLiuCosSlipFn::testDbSlipTime(void)
+{ // testDbSlipTime
+ const char* label = "database ABCD";
+ LiuCosSlipFn slipfn;
+
+ spatialdata::spatialdb::SimpleDB db(label);
+ slipfn.dbSlipTime(&db);
+
+ CPPUNIT_ASSERT(0 != slipfn._dbSlipTime);
+ CPPUNIT_ASSERT_EQUAL(std::string(label),
+ std::string(slipfn._dbSlipTime->label()));
+ CPPUNIT_ASSERT(0 == slipfn._dbFinalSlip);
+ CPPUNIT_ASSERT(0 == slipfn._dbRiseTime);
+} // testDbSlipTime
+
+// ----------------------------------------------------------------------
+// Test dbRiseTime().
+void
+pylith::faults::TestLiuCosSlipFn::testDbRiseTime(void)
+{ // testDbRiseTime
+ const char* label = "database ABCDE";
+ LiuCosSlipFn slipfn;
+
+ spatialdata::spatialdb::SimpleDB db(label);
+ slipfn.dbRiseTime(&db);
+
+ CPPUNIT_ASSERT(0 != slipfn._dbRiseTime);
+ CPPUNIT_ASSERT_EQUAL(std::string(label),
+ std::string(slipfn._dbRiseTime->label()));
+ CPPUNIT_ASSERT(0 == slipfn._dbFinalSlip);
+ CPPUNIT_ASSERT(0 == slipfn._dbSlipTime);
+} // testDbRiseTime
+
+// ----------------------------------------------------------------------
+// Test initialize() in 1-D.
+void
+pylith::faults::TestLiuCosSlipFn::testInitialize1D(void)
+{ // testInitialize1D
+ const char* meshFilename = "data/line2.mesh";
+ const char* faultLabel = "fault";
+ const int faultId = 2;
+ const char* finalSlipFilename = "data/line2_finalslip.spatialdb";
+ const char* slipTimeFilename = "data/line2_sliptime.spatialdb";
+ const char* riseTimeFilename = "data/line2_risetime.spatialdb";
+ const int constraintPts[] = { 3 };
+ const double finalSlipE[] = { 2.3 };
+ const double slipTimeE[] = { 1.2 };
+ const double riseTimeE[] = { 1.4 };
+ const int numConstraintPts = 1;
+
+ _TestLiuCosSlipFn::DataStruct data = {meshFilename,
+ faultLabel,
+ faultId,
+ finalSlipFilename,
+ slipTimeFilename,
+ riseTimeFilename,
+ constraintPts,
+ finalSlipE,
+ slipTimeE,
+ riseTimeE,
+ numConstraintPts};
+ _testInitialize(data);
+} // testInitialize1D
+
+// ----------------------------------------------------------------------
+// Test initialize() in 2-D.
+void
+pylith::faults::TestLiuCosSlipFn::testInitialize2D(void)
+{ // testInitialize2D
+ 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* riseTimeFilename = "data/tri3_risetime.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 riseTimeE[] = { 1.4, 1.5 };
+ const int numConstraintPts = 2;
+
+ _TestLiuCosSlipFn::DataStruct data = {meshFilename,
+ faultLabel,
+ faultId,
+ finalSlipFilename,
+ slipTimeFilename,
+ riseTimeFilename,
+ constraintPts,
+ finalSlipE,
+ slipTimeE,
+ riseTimeE,
+ numConstraintPts};
+ _testInitialize(data);
+} // testInitialize2D
+
+// ----------------------------------------------------------------------
+// Test initialize() in 3-D.
+void
+pylith::faults::TestLiuCosSlipFn::testInitialize3D(void)
+{ // testInitialize3D
+ const char* meshFilename = "data/tet4.mesh";
+ const char* faultLabel = "fault";
+ const int faultId = 2;
+ const char* finalSlipFilename = "data/tet4_finalslip.spatialdb";
+ const char* slipTimeFilename = "data/tet4_sliptime.spatialdb";
+ const char* riseTimeFilename = "data/tet4_risetime.spatialdb";
+ const int constraintPts[] = { 3, 4, 5 };
+ const double finalSlipE[] = { 2.3, -0.7, 0.1,
+ 2.4, -0.8, 0.2,
+ 2.5, -0.9, 0.3 };
+ const double slipTimeE[] = { 1.2, 1.3, 1.4 };
+ const double riseTimeE[] = { 1.5, 1.6, 1.7 };
+ const int numConstraintPts = 3;
+
+ _TestLiuCosSlipFn::DataStruct data = {meshFilename,
+ faultLabel,
+ faultId,
+ finalSlipFilename,
+ slipTimeFilename,
+ riseTimeFilename,
+ constraintPts,
+ finalSlipE,
+ slipTimeE,
+ riseTimeE,
+ numConstraintPts};
+ _testInitialize(data);
+} // testInitialize3D
+
+// ----------------------------------------------------------------------
+// Test slip().
+void
+pylith::faults::TestLiuCosSlipFn::testSlip(void)
+{ // testSlip
+ const double finalSlipE[] = { 2.3, 0.1,
+ 0.0, 0.0};
+ const double slipTimeE[] = { 1.2, 1.3 };
+ const double riseTimeE[] = { 1.4, 1.5 };
+ const double originTime = 5.064;
+
+ ALE::Obj<Mesh> faultMesh;
+ LiuCosSlipFn slipfn;
+ _initialize(&faultMesh, &slipfn, originTime);
+
+ const int spaceDim = faultMesh->getDimension() + 1;
+
+ const ALE::Obj<Mesh::label_sequence>& vertices = faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+ ALE::Obj<real_section_type> slip =
+ new real_section_type(faultMesh->comm(), faultMesh->debug());
+ slip->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(),
+ vertices->end()),
+ *std::max_element(vertices->begin(), vertices->end())+1));
+ slip->setFiberDimension(vertices, spaceDim);
+ faultMesh->allocate(slip);
+ CPPUNIT_ASSERT(!slip.isNull());
+
+ const double t = 2.134;
+ slipfn.slip(slip, originTime+t, faultMesh);
+
+ const double tolerance = 1.0e-06;
+ 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 slipNorm =
+ _slipFn(t - slipTimeE[iPoint], slipMag, riseTimeE[iPoint]) / slipMag;
+
+ 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] * slipNorm;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, vals[iDim], tolerance);
+ } // for
+ } // for
+} // testSlip
+
+// ----------------------------------------------------------------------
+// Test slipIncr().
+void
+pylith::faults::TestLiuCosSlipFn::testSlipIncr(void)
+{ // testSlipIncr
+ const double finalSlipE[] = { 2.3, 0.1,
+ 0.0, 0.0};
+ const double slipTimeE[] = { 1.2, 1.3 };
+ const double riseTimeE[] = { 1.4, 1.5 };
+ const double originTime = 1.064;
+
+ ALE::Obj<Mesh> faultMesh;
+ LiuCosSlipFn slipfn;
+ _initialize(&faultMesh, &slipfn, originTime);
+
+ const int spaceDim = faultMesh->getDimension() + 1;
+
+ const ALE::Obj<Mesh::label_sequence>& vertices = faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+ ALE::Obj<real_section_type> slip =
+ new real_section_type(faultMesh->comm(), faultMesh->debug());
+ slip->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(),
+ vertices->end()),
+ *std::max_element(vertices->begin(), vertices->end())+1));
+ slip->setFiberDimension(vertices, spaceDim);
+ faultMesh->allocate(slip);
+ CPPUNIT_ASSERT(!slip.isNull());
+
+ const double t0 = 1.234;
+ const double t1 = 3.635;
+ slipfn.slipIncr(slip, originTime+t0, originTime+t1, faultMesh);
+
+ const double tolerance = 1.0e-06;
+ 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 slipNorm0 =
+ _slipFn(t0 - slipTimeE[iPoint], slipMag, riseTimeE[iPoint]) / slipMag;
+ const double slipNorm1 =
+ _slipFn(t1 - slipTimeE[iPoint], slipMag, riseTimeE[iPoint]) / slipMag;
+
+ 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
+
+// ----------------------------------------------------------------------
+// Test _slip().
+void
+pylith::faults::TestLiuCosSlipFn::testSlipTH(void)
+{ // testSlipTH
+ const double t = 0.734;
+ const double finalSlip = 4.64;
+ const double riseTime = 3.23;
+
+ const double slipE = _slipFn(t, finalSlip, riseTime);
+
+ double slip = LiuCosSlipFn::_slipFn(t, finalSlip, riseTime);
+
+ const double tolerance = 1.0e-06;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(slipE, slip, tolerance);
+
+ slip = LiuCosSlipFn::_slipFn(-0.5, finalSlip, riseTime);
+ CPPUNIT_ASSERT_EQUAL(0.0, slip);
+
+ slip = LiuCosSlipFn::_slipFn(1.0e+10, finalSlip, riseTime);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(finalSlip, slip, tolerance);
+} // testSlipTH
+
+// ----------------------------------------------------------------------
+// Initialize LiuCosSlipFn.
+void
+pylith::faults::TestLiuCosSlipFn::_initialize(ALE::Obj<Mesh>* faultMesh,
+ LiuCosSlipFn* slipfn,
+ const double originTime)
+{ // _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* riseTimeFilename = "data/tri3_risetime.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;
+ (*faultMesh) = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
+ ALE::Obj<ALE::Mesh> faultBd = NULL;
+ CohesiveTopology::createFault(*faultMesh, faultBd,
+ mesh,
+ mesh->getIntSection(faultLabel));
+ CohesiveTopology::create(*faultMesh, faultBd, mesh,
+ mesh->getIntSection(faultLabel),
+ faultId,
+ useLagrangeConstraints);
+ 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 dbRiseTime("rise time");
+ spatialdata::spatialdb::SimpleIOAscii ioRiseTime;
+ ioRiseTime.filename(riseTimeFilename);
+ dbRiseTime.ioHandler(&ioRiseTime);
+
+ // setup LiuCosSlipFn
+ slipfn->dbFinalSlip(&dbFinalSlip);
+ slipfn->dbSlipTime(&dbSlipTime);
+ slipfn->dbRiseTime(&dbRiseTime);
+
+ slipfn->initialize(*faultMesh, &cs, originTime);
+} // _initialize
+
+// ----------------------------------------------------------------------
+// Test initialize().
+void
+pylith::faults::TestLiuCosSlipFn::_testInitialize(const _TestLiuCosSlipFn::DataStruct& data)
+{ // _testInitialize
+ typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
+
+ // Setup mesh
+ ALE::Obj<Mesh> mesh;
+ meshio::MeshIOAscii meshIO;
+ meshIO.filename(data.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 = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
+ ALE::Obj<ALE::Mesh> faultBd = NULL;
+ const bool useLagrangeConstraints = true;
+ CohesiveTopology::createFault(faultMesh, faultBd,
+ mesh,
+ mesh->getIntSection(data.faultLabel));
+ CohesiveTopology::create(faultMesh, faultBd, mesh,
+ mesh->getIntSection(data.faultLabel),
+ data.faultId,
+ useLagrangeConstraints);
+ 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(data.finalSlipFilename);
+ dbFinalSlip.ioHandler(&ioFinalSlip);
+
+ spatialdata::spatialdb::SimpleDB dbSlipTime("slip time");
+ spatialdata::spatialdb::SimpleIOAscii ioSlipTime;
+ ioSlipTime.filename(data.slipTimeFilename);
+ dbSlipTime.ioHandler(&ioSlipTime);
+
+ spatialdata::spatialdb::SimpleDB dbRiseTime("rise time");
+ spatialdata::spatialdb::SimpleIOAscii ioRiseTime;
+ ioRiseTime.filename(data.riseTimeFilename);
+ dbRiseTime.ioHandler(&ioRiseTime);
+
+ // setup LiuCosSlipFn
+ LiuCosSlipFn slipfn;
+ slipfn.dbFinalSlip(&dbFinalSlip);
+ slipfn.dbSlipTime(&dbSlipTime);
+ slipfn.dbRiseTime(&dbRiseTime);
+
+ const double originTime = 5.353;
+
+ slipfn.initialize(faultMesh, &cs, originTime);
+
+ 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) {
+ const int fiberDim = slipfn._parameters->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(spaceDim+2, fiberDim);
+
+ const real_section_type::value_type* vals =
+ slipfn._parameters->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.finalSlipE[iPoint*spaceDim+iDim],
+ vals[iDim],
+ tolerance);
+
+ const double riseTime = vals[spaceDim ];
+ const double slipTime = vals[spaceDim+1];
+
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.riseTimeE[iPoint], riseTime, tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.slipTimeE[iPoint]+originTime,
+ slipTime, tolerance);
+ } // for
+} // _testInitialize
+
+// ----------------------------------------------------------------------
+// Slip time function.
+double
+pylith::faults::TestLiuCosSlipFn::_slipFn(const double t,
+ const double finalSlip,
+ const double riseTime)
+{ // _slipFn
+ const float tau = riseTime * 1.525;
+ const float tau1 = 0.13 * tau;
+ const float tau2 = tau - tau1;
+ const float Cn =
+ M_PI / (1.4 * M_PI * tau1 + 1.2 * tau1 + 0.3 * M_PI * tau2);
+
+ double slip = 0.0;
+ if (t <= tau1) {
+ slip = 0.7*t - 0.7*tau1/M_PI*sin(M_PI*t/tau1)
+ - 0.6*tau1/(0.5*M_PI)*(cos(0.5*M_PI*t/tau1) - 1.0);
+ slip *= Cn;
+ } else if (t <= 2.0*tau1) {
+ slip = 1.0*t - 0.7*tau1/M_PI*sin(M_PI*t/tau1) +
+ 0.3*tau2/M_PI*sin(M_PI*(t-tau1)/tau2) +
+ 1.2*tau1/M_PI - 0.3*tau1;
+ slip *= Cn;
+ } else if (t <= tau) {
+ slip = 0.3*t + 0.3*tau2/M_PI*sin(M_PI*(t-tau1)/tau2) +
+ 1.1*tau1 + 1.2*tau1/M_PI;
+ slip *= Cn;
+ } else
+ slip = 1.0;
+ slip *= finalSlip;
+
+ return slip;
+} // _slipFn
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestLiuCosSlipFn.hh 2008-11-11 22:36:30 UTC (rev 13292)
@@ -0,0 +1,130 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/faults/TestLiuCosSlipFn.hh
+ *
+ * @brief C++ TestLiuCosSlipFn object
+ *
+ * C++ unit testing for LiuCosSlipFn.
+ */
+
+#if !defined(pylith_faults_testliucosslipfn_hh)
+#define pylith_faults_testliucosslipfn_hh
+
+#include "pylith/utils/sievetypes.hh" // USES Mesh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace faults {
+ class TestLiuCosSlipFn;
+ class LiuCosSlipFn;
+
+ namespace _TestLiuCosSlipFn {
+ struct DataStruct;
+ } // _LiuCosSlipTimeFn
+ } // faults
+} // pylith
+
+/// C++ unit testing for LiuCosSlipFn
+class pylith::faults::TestLiuCosSlipFn : public CppUnit::TestFixture
+{ // class TestLiuCosSlipFn
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestLiuCosSlipFn );
+
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testDbFinalSlip );
+ CPPUNIT_TEST( testDbSlipTime );
+ CPPUNIT_TEST( testDbRiseTime );
+ CPPUNIT_TEST( testInitialize1D );
+ CPPUNIT_TEST( testInitialize2D );
+ CPPUNIT_TEST( testInitialize3D );
+ CPPUNIT_TEST( testSlip );
+ CPPUNIT_TEST( testSlipIncr );
+ CPPUNIT_TEST( testSlipTH );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Test constructor.
+ void testConstructor(void);
+
+ /// Test dbFinalSlip().
+ void testDbFinalSlip(void);
+
+ /// Test dbSlipTime().
+ void testDbSlipTime(void);
+
+ /// Test dbRiseTime().
+ void testDbRiseTime(void);
+
+ /// Test initialize() in 1-D.
+ void testInitialize1D(void);
+
+ /// Test initialize() in 2-D.
+ void testInitialize2D(void);
+
+ /// Test initialize() in 3-D.
+ void testInitialize3D(void);
+
+ /// Test slip().
+ void testSlip(void);
+
+ /// Test slipIncr().
+ void testSlipIncr(void);
+
+ /// Test _slip().
+ void testSlipTH(void);
+
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /** Initialize LiuCosSlipFn.
+ *
+ * @param faultMesh Fault mesh.
+ * @param slipfn LiuCos slip function.
+ * @param originTime Origin time for earthquake rupture.
+ */
+ static
+ void _initialize(ALE::Obj<Mesh>* faultMesh,
+ LiuCosSlipFn* slipfn,
+ const double originTime);
+
+ /** Test intialize().
+ *
+ * @param data Data for initialization and testing of LiuCosSlipFn.
+ */
+ static
+ void _testInitialize(const _TestLiuCosSlipFn::DataStruct& data);
+
+ /** Slip time function.
+ *
+ * @param t Time relative to when slip begins.
+ * @param finalSlip Final slip.
+ * @param riseTime Rise time (t95).
+ */
+ static
+ double _slipFn(const double t,
+ const double finalSlip,
+ const double riseTime);
+
+}; // class TestLiuCosSlipFn
+
+#endif // pylith_faults_testliucosslipfn_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/Makefile.am 2008-11-11 22:24:06 UTC (rev 13291)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/Makefile.am 2008-11-11 22:36:30 UTC (rev 13292)
@@ -19,6 +19,7 @@
line2_sliptime.spatialdb \
line2_peakrate.spatialdb \
line2_sliprate.spatialdb \
+ line2_risetime.spatialdb \
tri3.mesh \
tri3b.mesh \
tri3c.mesh \
@@ -30,6 +31,7 @@
tri3_sliptime.spatialdb \
tri3_peakrate.spatialdb \
tri3_sliprate.spatialdb \
+ tri3_risetime.spatialdb \
tri3d_finalslip.spatialdb \
tri3d_sliptime.spatialdb \
tri3d_peakrate.spatialdb \
@@ -61,6 +63,7 @@
tet4_sliptime.spatialdb \
tet4_peakrate.spatialdb \
tet4_sliprate.spatialdb \
+ tet4_risetime.spatialdb \
tet4e_finalslip.spatialdb \
tet4e_sliptime.spatialdb \
tet4e_peakrate.spatialdb \
Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/line2_risetime.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/line2_risetime.spatialdb (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/line2_risetime.spatialdb 2008-11-11 22:36:30 UTC (rev 13292)
@@ -0,0 +1,14 @@
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 1
+ value-names = rise-time
+ value-units = s
+ num-locs = 1
+ data-dim = 0
+ space-dim = 1
+ cs-data = cartesian {
+ to-meters = 1.0
+ space-dim = 1
+ }
+}
+0.0 1.4
Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/tet4_risetime.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/tet4_risetime.spatialdb (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/tet4_risetime.spatialdb 2008-11-11 22:36:30 UTC (rev 13292)
@@ -0,0 +1,16 @@
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 1
+ value-names = rise-time
+ value-units = s
+ num-locs = 3
+ data-dim = 2
+ space-dim = 3
+ cs-data = cartesian {
+ to-meters = 1.0
+ space-dim = 3
+ }
+}
+0.0 -1.0 0.0 1.5
+0.0 0.0 1.0 1.6
+0.0 1.0 0.0 1.7
Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/tri3_risetime.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/tri3_risetime.spatialdb (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/tri3_risetime.spatialdb 2008-11-11 22:36:30 UTC (rev 13292)
@@ -0,0 +1,15 @@
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 1
+ value-names = rise-time
+ value-units = s
+ num-locs = 2
+ data-dim = 1
+ space-dim = 2
+ cs-data = cartesian {
+ to-meters = 1.0
+ space-dim = 2
+ }
+}
+0.0 +1.0 1.4
+0.0 -1.0 1.5
More information about the CIG-COMMITS
mailing list