[cig-commits] r14595 - in short/3D/PyLith/branches/pylith-swig: . doc/install libsrc libsrc/faults libsrc/topology
brad at geodynamics.org
brad at geodynamics.org
Sat Apr 4 16:31:38 PDT 2009
Author: brad
Date: 2009-04-04 16:31:36 -0700 (Sat, 04 Apr 2009)
New Revision: 14595
Added:
short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/faultsfwd.hh
Modified:
short/3D/PyLith/branches/pylith-swig/Makefile.am
short/3D/PyLith/branches/pylith-swig/doc/install/MacBookPro_Aagaard.txt
short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.icc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.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/ConstRateSlipFn.icc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.icc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/faults/SlipTimeFn.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.icc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Fields.hh
Log:
Worked on updating fault stuff.
Modified: short/3D/PyLith/branches/pylith-swig/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/Makefile.am 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/Makefile.am 2009-04-04 23:31:36 UTC (rev 14595)
@@ -10,7 +10,7 @@
# ----------------------------------------------------------------------
#
-ACLOCAL_AMFLAGS = -I ./m4
+ACLOCAL_AMFLAGS = -I m4
SUBDIRS = \
libsrc \
Modified: short/3D/PyLith/branches/pylith-swig/doc/install/MacBookPro_Aagaard.txt
===================================================================
--- short/3D/PyLith/branches/pylith-swig/doc/install/MacBookPro_Aagaard.txt 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/doc/install/MacBookPro_Aagaard.txt 2009-04-04 23:31:36 UTC (rev 14595)
@@ -96,8 +96,8 @@
To create libmpich.dylib:
cd ${TOOLS_DIR}/mpich2-1.0.5/${TOOLS_FORMAT}/lib
- mkdir lib
- cd lib
+ mkdir tmp
+ cd tmp
ar -x ../libmpich.a
gcc -dynamiclib -single_module -undefined dynamic_lookup -o ${TOOLS_DIR}/mpich2-1.0.5/${TOOLS_FORMAT}/lib/libmpich.dylib *.o
rm *.o
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-04-04 23:31:36 UTC (rev 14595)
@@ -28,6 +28,10 @@
bc/DirichletBoundary.cc \
bc/Neumann.cc \
bc/AbsorbingDampers.cc \
+ faults/Fault.cc \
+ faults/SlipTimeFn.cc \
+ faults/StepSlipFn.cc \
+ faults/ConstRateSlipFn.cc \
feassemble/CellGeometry.cc \
feassemble/Constraint.cc \
feassemble/GeometryPoint1D.cc \
@@ -89,16 +93,14 @@
utils/EventLogger.cc \
utils/TestArray.cc
+# faults/TopologyOps.cc \
# faults/BruneSlipFn.cc \
-# faults/ConstRateSlipFn.cc \
# faults/CohesiveTopology.cc \
# faults/EqKinSrc.cc \
-# faults/Fault.cc \
# faults/FaultCohesive.cc \
# faults/FaultCohesiveKin.cc \
# faults/FaultCohesiveDyn.cc \
# faults/LiuCosSlipFn.cc \
-# faults/SlipTimeFn.cc \
# faults/StepSlipFn.cc \
# materials/GenMaxwellIsotropic3D.cc \
# meshio/CellFilter.cc \
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.hh 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.hh 2009-04-04 23:31:36 UTC (rev 14595)
@@ -24,27 +24,14 @@
#if !defined(pylith_faults_bruneslipfn_hh)
#define pylith_faults_bruneslipfn_hh
+// Include directives ---------------------------------------------------
#include "SlipTimeFn.hh"
-/// Namespace for pylith package
-namespace pylith {
- namespace faults {
- class BruneSlipFn;
- class TestBruneSlipFn; // 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 Brune slip time function.
+// SlipTimeFn -----------------------------------------------------------
class pylith::faults::BruneSlipFn : public SlipTimeFn
{ // class BruneSlipFn
friend class TestBruneSlipFn; // unit testing
@@ -71,21 +58,20 @@
*/
void dbSlipTime(spatialdata::spatialdb::SpatialDB* const db);
- /** Set spatial database for peak slip rate.
+ /** Set spatial database for rise time (0 -> 0.95 final slip).
*
* @param db Spatial database
*/
- void dbPeakRate(spatialdata::spatialdb::SpatialDB* const db);
+ 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 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,84 +79,78 @@
*
* @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
- BruneSlipFn(const BruneSlipFn& m);
+ BruneSlipFn(const BruneSlipFn&); ///< Not implemented
+ const BruneSlipFn& operator=(const BruneSlipFn&); ///< Not implemented
- /// Not implemented
- const BruneSlipFn& operator=(const BruneSlipFn& 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 peakRate Peak slip rate at point
+ * @param t Time relative to slip starting time at point.
+ * @param finalSlip Final slip at point.
+ * @param riseTime Peak slip rate at point.
*
* @returns Slip at point at time t
*/
static
double _slipFn(const double t,
const double finalSlip,
- const double peakRate);
+ const double riseTime);
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
- /// Parameters for Brune 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 _finalSlipVertex; ///< Final slip at a vertex.
- /// Spatial database for final slip
+ /// Parameters for Brune slip time function, final slip (vector),
+ /// rise time (scalar), slip time (scalar).
+ topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
+
+ /// Spatial database for final slip.
spatialdata::spatialdb::SpatialDB* _dbFinalSlip;
- /// Spatial database for slip time
+ /// Spatial database for slip time.
spatialdata::spatialdb::SpatialDB* _dbSlipTime;
- /// Spatial database for peak slip rate
- spatialdata::spatialdb::SpatialDB* _dbPeakRate;
+ /// Spatial database for rise time (0 -> 0.95 final slip).
+ spatialdata::spatialdb::SpatialDB* _dbRiseTime;
- int _spaceDim; ///< Spatial dimension for slip field.
-
}; // class BruneSlipFn
#include "BruneSlipFn.icc" // inline methods
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.icc 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.icc 2009-04-04 23:31:36 UTC (rev 14595)
@@ -31,11 +31,11 @@
_dbSlipTime = db;
} // dbSlipTime
-// Set spatial database for peak slip rate.
+// Set spatial database for rise time.
inline
void
-pylith::faults::BruneSlipFn::dbPeakRate(spatialdata::spatialdb::SpatialDB* const db) {
- _dbPeakRate = db;
+pylith::faults::BruneSlipFn::dbRiseTime(spatialdata::spatialdb::SpatialDB* const db) {
+ _dbRiseTime = db;
} // dbPeakRate
// Compute slip using slip time function.
@@ -43,10 +43,11 @@
double
pylith::faults::BruneSlipFn::_slipFn(const double t,
const double finalSlip,
- const double peakRate) {
+ const double riseTime) {
double slip = 0.0;
if (t > 0.0) {
- assert(peakRate > 0.0);
+ assert(riseTime > 0.0);
+ const double peakRate = finalSlip / riseTime * 1.745;
const double tau =
// prevent 0 == tau when 0 == finalSlip
(finalSlip > 0.0) ? finalSlip / (exp(1.0) * peakRate) : 0.1;
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.hh 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.hh 2009-04-04 23:31:36 UTC (rev 14595)
@@ -18,146 +18,18 @@
#if !defined(pylith_faults_cohesivetopology_hh)
#define pylith_faults_cohesivetopology_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 CohesiveTopology;
- } // faults
-} // pylith
-
-/// C++ object to manage creation of cohesive cells.
+// CohesiveTopology -----------------------------------------------------
class pylith::faults::CohesiveTopology
{ // class CohesiveTopology
public :
- typedef std::set<Mesh::point_type> PointSet;
- typedef std::vector<sieve_type::point_type> PointArray;
+ typedef std::set<SieveMesh::point_type> PointSet;
+ typedef std::vector<sieve_type::point_type> PointArray;
typedef std::pair<sieve_type::point_type, int> oPoint_type;
- typedef std::vector<oPoint_type> oPointArray;
+ typedef std::vector<oPoint_type> oPointArray;
-protected:
- template<typename Sieve, typename Renumbering>
- class ReplaceVisitor {
- public:
- typedef typename Sieve::point_type point_type;
- protected:
- Renumbering& renumbering;
- const int size;
- int i;
- const int debug;
- point_type *points;
- bool mapped;
- public:
- ReplaceVisitor(Renumbering& r, const int size, const int debug = 0) : renumbering(r), size(size), i(0), debug(debug) {
- this->points = new point_type[this->size];
- this->mapped = false;
- };
- ~ReplaceVisitor() {delete [] this->points;};
- void visitPoint(const point_type& point) {
- if (i >= this->size) {throw ALE::Exception("Too many points for ReplaceVisitor");}
- if (this->renumbering.find(point) != this->renumbering.end()) {
- if (debug) std::cout << " point " << this->renumbering[point] << std::endl;
- points[i] = this->renumbering[point];
- this->mapped = true;
- } else {
- if (debug) std::cout << " point " << point << std::endl;
- points[i] = point;
- }
- ++i;
- };
- void visitArrow(const typename Sieve::arrow_type&) {};
- public:
- const point_type *getPoints() {return this->points;};
- bool mappedPoint() {return this->mapped;};
- void clear() {this->i = 0; this->mapped = false;};
- };
- template<typename Sieve>
- class ClassifyVisitor {
- public:
- typedef typename Sieve::point_type point_type;
- protected:
- const Sieve& sieve;
- const PointSet& replaceCells;
- const PointSet& noReplaceCells;
- const point_type firstCohesiveCell;
- const int faceSize;
- const int debug;
- PointSet vReplaceCells;
- PointSet vNoReplaceCells;
- bool modified;
- bool setupMode;
- int size;
- ALE::ISieveVisitor::PointRetriever<Sieve> pR;
- public:
- ClassifyVisitor(const Sieve& s, const PointSet& rC, const PointSet& nrC, const point_type& fC, const int fS, const int debug = 0) : sieve(s), replaceCells(rC), noReplaceCells(nrC), firstCohesiveCell(fC), faceSize(fS), debug(debug), modified(false), setupMode(true), size(0) {
- pR.setSize(s.getMaxConeSize());
- };
- ~ClassifyVisitor() {};
- void visitPoint(const point_type& point) {
- if (this->setupMode) {
- if (replaceCells.find(point) != replaceCells.end()) vReplaceCells.insert(point);
- if (noReplaceCells.find(point) != noReplaceCells.end()) vNoReplaceCells.insert(point);
- if (point >= firstCohesiveCell) return;
- this->modified = true;
- this->size++;
- return;
- }
- bool classified = false;
-
- if (debug) {std::cout << "Checking neighbor " << point << std::endl;}
- if (vReplaceCells.find(point) != vReplaceCells.end()) {
- if (debug) {std::cout << " already in replaceCells" << std::endl;}
- return;
- }
- if (vNoReplaceCells.find(point) != vNoReplaceCells.end()) {
- if (debug) {std::cout << " already in noReplaceCells" << std::endl;}
- return;
- }
- if (point >= firstCohesiveCell) {
- if (debug) {std::cout << " already a cohesive cell" << std::endl;}
- return;
- }
- // If neighbor shares a face with anyone in replaceCells, then add
- for(PointSet::const_iterator c_iter = vReplaceCells.begin(); c_iter != vReplaceCells.end(); ++c_iter) {
- sieve.meet(*c_iter, point, pR);
-
- if (pR.getSize() == faceSize) {
- if (debug) {std::cout << " Scheduling " << point << " for replacement" << std::endl;}
- vReplaceCells.insert(point);
- modified = true;
- classified = true;
- pR.clear();
- break;
- }
- pR.clear();
- }
- if (classified) return;
- // It is unclear whether taking out the noReplace cells will speed this up
- for(PointSet::const_iterator c_iter = vNoReplaceCells.begin(); c_iter != vNoReplaceCells.end(); ++c_iter) {
- sieve.meet(*c_iter, point, pR);
-
- if (pR.getSize() == faceSize) {
- if (debug) {std::cout << " Scheduling " << point << " for no replacement" << std::endl;}
- vNoReplaceCells.insert(point);
- modified = true;
- classified = true;
- pR.clear();
- break;
- }
- pR.clear();
- }
- };
- void visitArrow(const typename Sieve::arrow_type&) {};
- public:
- const PointSet& getReplaceCells() const {return this->vReplaceCells;};
- const PointSet& getNoReplaceCells() const {return this->vNoReplaceCells;};
- const bool getModified() const {return this->modified;};
- const int getSize() const {return this->size;};
- void setMode(const bool isSetup) {this->setupMode = isSetup;};
- void reset() {this->modified = false;};
- };
-
// PUBLIC METHODS /////////////////////////////////////////////////////
public :
/** Create the fault mesh.
@@ -168,10 +40,10 @@
* fault surface
*/
static
- void createFault(Obj<SubMesh>& ifault,
- Obj<ALE::Mesh>& faultBd,
- const Obj<Mesh>& mesh,
- const Obj<Mesh::int_section_type>& groupField,
+ void createFault(topology::SubMesh* ifault,
+ ALE::Obj<ALE::Mesh>& faultBd,
+ const topology::Mesh& mesh,
+ const ALE::Obj<topology::Mesh::IntSection>& groupField,
const bool flipFault =false);
/** Create cohesive cells.
@@ -183,10 +55,10 @@
* Lagrange multipliers that require extra vertices, false otherwise
*/
static
- void create(Obj<SubMesh>& ifault,
- const Obj<ALE::Mesh>& faultBd,
- const Obj<Mesh>& mesh,
- const Obj<Mesh::int_section_type>& groupField,
+ void create(topology::SubMesh* ifault,
+ const ALE::Obj<ALE::Mesh>& faultBd,
+ const topology::Mesh& mesh,
+ const ALE::Obj<topology::Mesh::IntSection>& groupField,
const int materialId,
const bool constraintCell =false);
@@ -200,107 +72,12 @@
* Lagrange multipliers that require extra vertices, false otherwise.
*/
static
- void createParallel(ALE::Obj<SubMesh>* ifault,
+ void createParallel(topology::SubMesh* ifault,
std::map<Mesh::point_type, Mesh::point_type>* cohesiveToFault,
- const ALE::Obj<Mesh>& mesh,
+ const topology::Mesh& mesh,
const int materialId,
const bool constraintCell =false);
- // PRIVATE METHODS ////////////////////////////////////////////////////
-private :
- /** Get number of vertices on face.
- *
- * @param cell Finite-element cell
- * @param mesh Finite-element mesh
- *
- * @returns Number of vertices on cell face
- */
- static
- unsigned int _numFaceVertices(const Mesh::point_type& cell,
- const ALE::Obj<Mesh>& mesh,
- const int depth =-1);
-
- /** Determine a face orientation
- * We should really have an interpolated mesh, instead of
- * calculating this on the fly.
- *
- * @param cell Finite-element cell
- * @param mesh Finite-element mesh
- *
- * @returns True for positive orientation, otherwise false
- */
- static
- bool _faceOrientation(const Mesh::point_type& cell,
- const ALE::Obj<Mesh>& mesh,
- const int numCorners,
- const int indices[],
- const int oppositeVertex,
- PointArray *origVertices,
- PointArray *faceVertices);
-
- template<typename FaceType>
- static
- bool _getOrientedFace(const ALE::Obj<Mesh>& mesh,
- const Mesh::point_type& cell,
- FaceType face,
- const int numCorners,
- int indices[],
- PointArray *origVertices,
- PointArray *faceVertices);
-
- template<class InputPoints>
- static
- bool _compatibleOrientation(const ALE::Obj<Mesh>& mesh,
- const Mesh::point_type& p,
- const Mesh::point_type& q,
- const int numFaultCorners,
- const int faultFaceSize,
- const int faultDepth,
- const Obj<InputPoints>& points,
- int indices[],
- PointArray *origVertices,
- PointArray *faceVertices,
- PointArray *neighborVertices);
-
- static
- void _computeCensoredDepth(const ALE::Obj<Mesh::label_type>& depth,
- const ALE::Obj<Mesh::sieve_type>& sieve,
- const Mesh::point_type& firstCohesiveCell);
-
- static
- void classifyCells(const ALE::Obj<Mesh::sieve_type>& sieve,
- const Mesh::point_type& vertex,
- const int depth,
- const int faceSize,
- const Mesh::point_type& firstCohesiveCell,
- PointSet& replaceCells,
- PointSet& noReplaceCells,
- const int debug);
-
- static
- void createFaultSieveFromVertices(const int dim,
- const int firstCell,
- const PointSet& faultVertices,
- const Obj<Mesh>& mesh,
- const Obj<ALE::Mesh::arrow_section_type>& orientation,
- const Obj<ALE::Mesh::sieve_type>& faultSieve,
- const bool flipFault);
-
- static
- void createFaultSieveFromFaces(const int dim,
- const int firstCell,
- const int numFaces,
- const int faultVertices[],
- const int faultCells[],
- const Obj<Mesh>& mesh,
- const Obj<ALE::Mesh::arrow_section_type>& orientation,
- const Obj<ALE::Mesh::sieve_type>& faultSieve);
-
- static
- void orientFaultSieve(const int dim,
- const Obj<Mesh>& mesh,
- const Obj<ALE::Mesh::arrow_section_type>& orientation,
- const Obj<ALE::Mesh>& fault);
}; // class CohesiveTopology
#endif // pylith_faults_cohesivetopology_hh
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.cc 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.cc 2009-04-04 23:31:36 UTC (rev 14595)
@@ -14,8 +14,9 @@
#include "ConstRateSlipFn.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,20 +26,18 @@
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
-namespace pylith {
- namespace faults {
- namespace _ConstRateSlipFn {
- const int offsetSlipTime = 0;
- } // _ConstRateSlipFn
- } // 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::ConstRateSlipFn::ConstRateSlipFn(void) :
+ _slipTimeVertex(0),
+ _parameters(0),
_dbSlipRate(0),
- _dbSlipTime(0),
- _spaceDim(0)
+ _dbSlipTime(0)
{ // constructor
} // constructor
@@ -46,6 +45,7 @@
// Destructor.
pylith::faults::ConstRateSlipFn::~ConstRateSlipFn(void)
{ // destructor
+ delete _parameters; _parameters = 0;
_dbSlipRate = 0;
_dbSlipTime = 0;
} // destructor
@@ -54,37 +54,49 @@
// Initialize slip time function.
void
pylith::faults::ConstRateSlipFn::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 != _dbSlipRate);
assert(0 != _dbSlipTime);
- _spaceDim = cs->spaceDim();
- const int spaceDim = _spaceDim;
- const int indexSlipRate = 0;
- const int indexSlipTime = spaceDim + _ConstRateSlipFn::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();
+ const double velocityScale =
+ normalizer.lengthScale() / 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 + 1;
- _parameters = new real_section_type(faultMesh->comm(), faultMesh->debug());
- _parameters->addSpace(); // slip rate
- _parameters->addSpace(); // slip time
- assert(2 == _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); // slip rate
- _parameters->setFiberDimension(vertices, 1, 1); // slip time
- faultMesh->allocate(_parameters);
- assert(!_parameters.isNull());
+ delete _parameters; _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(faultMesh);
+ assert(0 != _parameters);
+ _parameters->add("slip rate", "slip rate");
+ topology::Field<topology::SubMesh>& slipRate = _parameters->get("slip rate");
+ slipRate.newSection(vertices, spaceDim);
+ slipRate.allocate();
+ slipRate.scale(velocityScale);
+ slipRate.vectorFieldType(topology::FieldBase::VECTOR);
+ const ALE::Obj<RealSection>& slipRateSection = slipRate.section();
+ assert(!slipRateSection.isNull());
+ _parameters->add("slip time", "slip time");
+ topology::Field<topology::SubMesh>& slipTime = _parameters->get("slip time");
+ slipTime.newSection(slipRateSection->getChart(), 1);
+ slipTime.allocate();
+ slipTime.scale(timeScale);
+ slipTime.vectorFieldType(topology::FieldBase::SCALAR);
+ const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
+ assert(!slipTimeSection.isNull());
+
// Open databases and set query values
_dbSlipRate->open();
switch (spaceDim)
@@ -114,18 +126,13 @@
_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);
+ _slipRateVertex.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,
@@ -133,7 +140,7 @@
normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(),
lengthScale);
- int err = _dbSlipRate->query(¶msVertex[indexSlipRate], spaceDim,
+ int err = _dbSlipRate->query(&_slipRateVertex[0], _slipRateVertex.size(),
&vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -143,10 +150,10 @@
msg << ") using spatial database " << _dbSlipRate->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(¶msVertex[indexSlipRate], spaceDim,
+ normalizer.nondimensionalize(&_slipRateVertex[0], _slipRateVertex.size(),
velocityScale);
- err = _dbSlipTime->query(¶msVertex[indexSlipTime], 1,
+ err = _dbSlipTime->query(&_slipTimeVertex, 1,
&vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -156,12 +163,12 @@
msg << ") using spatial database " << _dbSlipTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(¶msVertex[indexSlipTime], 1,
- timeScale);
+ normalizer.nondimensionalize(&_slipTimeVertex, 1, timeScale);
// add origin time to rupture time
- paramsVertex[indexSlipTime] += originTime;
+ _slipTimeVertex += originTime;
- _parameters->updatePoint(*v_iter, ¶msVertex[0]);
+ slipRateSection->updatePoint(*v_iter, &_slipRateVertex[0]);
+ slipTimeSection->updatePoint(*v_iter, &_slipTimeVertex);
} // for
// Close databases
@@ -172,114 +179,117 @@
// ----------------------------------------------------------------------
// Get slip on fault surface at time t.
void
-pylith::faults::ConstRateSlipFn::slip(const ALE::Obj<pylith::real_section_type>& slipField,
- const double t,
- const ALE::Obj<Mesh>& faultMesh)
+pylith::faults::ConstRateSlipFn::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 indexSlipRate = 0;
- const int indexSlipTime = spaceDim + _ConstRateSlipFn::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>& slipRate =
+ _parameters->get("final slip");
+ const ALE::Obj<RealSection>& slipRateSection = slipRate.section();
+ assert(!slipRateSection.isNull());
+ const topology::Field<topology::SubMesh>& slipTime =
+ _parameters->get("slip time");
+ const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
+ assert(!slipTimeSection.isNull());
+ const ALE::Obj<RealSection>& slipSection = slip->section();
+ assert(!slipSection.isNull());
+
+ 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);
+ slipRateSection->restrictPoint(*v_iter, &_slipRateVertex[0],
+ _slipRateVertex.size());
+ slipTimeSection->restrictPoint(*v_iter, &_slipTimeVertex, 1);
- const double* slipRate = ¶msVertex[indexSlipRate];
- const double slipTime = paramsVertex[indexSlipTime];
- slipValues = 0.0;
+ const double relTime = t - _slipTimeVertex;
+ const double elapsedTime = (relTime > 0) ? relTime : 0.0;
+ _slipRateVertex *= elapsedTime; // Convert slip rate to slip
- const double relTime = t - slipTime;
- if (relTime > 0)
- for (int i=0; i < spaceDim; ++i)
- slipValues[i] = slipRate[i] * relTime;
-
// Update field
- slipField->updateAddPoint(*v_iter, &slipValues[0]);
+ slipSection->updateAddPoint(*v_iter, &_slipRateVertex[0]);
} // for
- PetscLogFlops(numVertices * (1+1 + 4*spaceDim));
+ PetscLogFlops(vertices->size() * (1 + _slipRateVertex.size()));
} // slip
// ----------------------------------------------------------------------
// Get increment of slip on fault surface between time t0 and t1.
void
-pylith::faults::ConstRateSlipFn::slipIncr(const ALE::Obj<pylith::real_section_type>& slipField,
- const double t0,
- const double t1,
- const ALE::Obj<Mesh>& faultMesh)
+pylith::faults::ConstRateSlipFn::slipIncr(
+ topology::Field<topology::SubMesh>* slip,
+ const double t0,
+ 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 indexSlipRate = 0;
- const int indexSlipTime = spaceDim + _ConstRateSlipFn::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>& slipRate =
+ _parameters->get("final slip");
+ const ALE::Obj<RealSection>& slipRateSection = slipRate.section();
+ assert(!slipRateSection.isNull());
+ const topology::Field<topology::SubMesh>& slipTime =
+ _parameters->get("slip time");
+ const ALE::Obj<RealSection>& slipTimeSection = slipTime.section();
+ assert(!slipTimeSection.isNull());
+ const ALE::Obj<RealSection>& slipSection = slip->section();
+ assert(!slipSection.isNull());
+
+ 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);
+ slipRateSection->restrictPoint(*v_iter, &_slipRateVertex[0],
+ _slipRateVertex.size());
+ slipTimeSection->restrictPoint(*v_iter, &_slipTimeVertex, 1);
- const double* slipRate = ¶msVertex[indexSlipRate];
- const double slipTime = paramsVertex[indexSlipTime];
- slipValues = 0.0;
-
- const double relTime0 = t0 - slipTime;
- const double relTime1 = t1 - slipTime;
+ const double relTime0 = t0 - _slipTimeVertex;
+ const double relTime1 = t1 - _slipTimeVertex;
double elapsedTime = 0.0;
if (relTime0 > 0)
elapsedTime = t1 - t0;
else if (relTime1 > 0)
- elapsedTime = t1 - slipTime;
- for (int i=0; i < spaceDim; ++i)
- slipValues[i] = slipRate[i] * elapsedTime;
+ elapsedTime = t1 - _slipTimeVertex;
+ _slipRateVertex *= elapsedTime; // Convert slip rate to slip
// Update field
- slipField->updateAddPoint(*v_iter, &slipValues[0]);
+ slipSection->updateAddPoint(*v_iter, &_slipRateVertex[0]);
} // for
- PetscLogFlops(numVertices * (2 + spaceDim));
+ PetscLogFlops(vertices->size() * (4 + _slipRateVertex.size()));
} // slipIncr
// ----------------------------------------------------------------------
-// Get final slip (slip rate in this case).
-ALE::Obj<pylith::real_section_type>
+// Get final slip.
+const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::ConstRateSlipFn::finalSlip(void)
{ // finalSlip
- // This is actually slip rate.
- return _parameters->getFibration(0);
+ // Slip rate is parameter instead of final slip.
+ return _parameters->get("slip rate");
} // finalSlip
// ----------------------------------------------------------------------
// Get time when slip begins at each point.
-ALE::Obj<pylith::real_section_type>
+const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::ConstRateSlipFn::slipTime(void)
{ // slipTime
- return _parameters->getFibration(1);
+ return _parameters->get("slip time");
} // slipTime
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.hh 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.hh 2009-04-04 23:31:36 UTC (rev 14595)
@@ -23,27 +23,14 @@
#if !defined(pylith_faults_constrateslipfn_hh)
#define pylith_faults_constrateslipfn_hh
+// Include directives ---------------------------------------------------
#include "SlipTimeFn.hh"
-/// Namespace for pylith package
-namespace pylith {
- namespace faults {
- class ConstRateSlipFn;
- class TestConstRateSlipFn; // 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 ConstRate slip time function.
+// SlipTimeFn -----------------------------------------------------------
class pylith::faults::ConstRateSlipFn : public SlipTimeFn
{ // class ConstRateSlipFn
friend class TestConstRateSlipFn; // unit testing
@@ -73,12 +60,11 @@
/** Initialize slip time function.
*
* @param faultMesh Finite-element mesh of fault.
- * @param cs Coordinate system for mesh.
+ * @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);
@@ -86,64 +72,58 @@
*
* @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
- ConstRateSlipFn(const ConstRateSlipFn& m);
+ ConstRateSlipFn(const ConstRateSlipFn&); ///< Not implemented
+ const ConstRateSlipFn& operator=(const ConstRateSlipFn&); ///< Not implemented
- /// Not implemented
- const ConstRateSlipFn& operator=(const ConstRateSlipFn& f);
-
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
- /// Parameters for ConstRate slip time function.
- /// Slip rate (vector), slip time (scalar).
- ALE::Obj<real_section_type> _parameters;
+ double _slipTimeVertex; ///< Slip time at a vertex.
+ double_array _slipRateVertex; ///< Slip rate at a vertex.
+ /// Parameters for constant slip rate slip time function, slip rate
+ /// (vector) and slip time (scalar).
+ topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
+
/// Spatial database for slip rate.
spatialdata::spatialdb::SpatialDB* _dbSlipRate;
/// Spatial database for slip time.
spatialdata::spatialdb::SpatialDB* _dbSlipTime;
- int _spaceDim; ///< Spatial dimension for slip field.
-
}; // class ConstRateSlipFn
#include "ConstRateSlipFn.icc" // inline methods
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.icc 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/ConstRateSlipFn.icc 2009-04-04 23:31:36 UTC (rev 14595)
@@ -14,9 +14,6 @@
#error "ConstRateSlipFn.icc can only be included from ConstRateSlipFn.hh"
#endif
-#include <math.h> // USES exp()
-#include <assert.h> // USES assert()
-
// Set spatial database for slip rate.
inline
void
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.cc 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.cc 2009-04-04 23:31:36 UTC (rev 14595)
@@ -14,6 +14,8 @@
#include "Fault.hh" // implementation of object methods
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::Fault::Fault(void) :
@@ -27,14 +29,15 @@
// Destructor.
pylith::faults::Fault::~Fault(void)
{ // destructor
+ delete _faultMesh; _faultMesh = 0;
} // destructor
// ----------------------------------------------------------------------
// Get mesh associated with fault fields.
-const ALE::Obj<pylith::SubMesh>&
+const pylith::topology::SubMesh&
pylith::faults::Fault:: faultMesh(void) const
{ // faultMesh
- return _faultMesh;
+ return *_faultMesh;
} // faultMesh
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.hh 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.hh 2009-04-04 23:31:36 UTC (rev 14595)
@@ -24,36 +24,17 @@
#if !defined(pylith_faults_fault_hh)
#define pylith_faults_fault_hh
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh, real_section_type
+// Include directives ---------------------------------------------------
+#include "faultsfwd.hh" // forward declarations
+
+#include "pylith/topology/topologyfwd.hh" // USES Field<SubMesh>, SubMesh
#include "pylith/utils/arrayfwd.hh" // USES double_array
-#include "pylith/utils/vectorfields.hh" // USES VectorFieldEnum
+#include "spatialdata/spatialdb/spatialdbfwd.hh" // USES SpatialDB
+
#include <string> // HASA std::string
-/// Namespace for pylith package
-namespace pylith {
- namespace faults {
- class Fault;
- class TestFault; // unit testing
- } // faults
-
- namespace topology {
- class FieldsManager;
- } // topology
-} // pylith
-
-/// Namespace for spatialdata package
-namespace spatialdata {
- namespace geocoords {
- class CoordSys;
- } // geocoords
-
- namespace spatialdb {
- class SpatialDB; // USES SpatialDB
- } // spatialdb
-} // spatialdata
-
-/// C++ abstract base class for Fault object.
+// Fault ----------------------------------------------------------------
class pylith::faults::Fault
{ // class Fault
friend class TestFault; // unit testing
@@ -90,14 +71,15 @@
*
* @returns Label of fault
*/
- const std::string& label(void) const;
+ const char* label(void) const;
/** Adjust mesh topology for fault implementation.
*
* @param mesh PETSc mesh
*/
virtual
- void adjustTopology(const ALE::Obj<Mesh>& mesh, const bool flipFault = false) = 0;
+ void adjustTopology(const topology::Mesh& mesh,
+ const bool flipFault =false) = 0;
/** Initialize fault. Determine orientation and setup boundary
* condition parameters.
@@ -114,8 +96,7 @@
* (used to improve conditioning of Jacobian matrix)
*/
virtual
- void initialize(const ALE::Obj<Mesh>& mesh,
- const spatialdata::geocoords::CoordSys* cs,
+ void initialize(const topology::Mesh& mesh,
const double_array& upDir,
const double_array& normalDir,
spatialdata::spatialdb::SpatialDB* matDB) = 0;
@@ -124,51 +105,40 @@
*
* @returns PETSc mesh object
*/
- const ALE::Obj<SubMesh>& faultMesh(void) const;
+ const topology::SubMesh& faultMesh(void) const;
/** Get vertex field associated with integrator.
*
- * @param fieldType Type of field.
* @param name Name of vertex field.
- * @param mesh PETSc mesh for problem.
- * @param fields Fields manager.
+ * @param fields Solution fields.
* @returns Vertex field.
*/
virtual
- const ALE::Obj<real_section_type>&
- vertexField(VectorFieldEnum* fieldType,
- const char* name,
- const ALE::Obj<Mesh>& mesh,
- topology::FieldsManager* const fields) = 0;
+ const topology::Field<topology::SubMesh>&
+ vertexField(const char* name,
+ const topology::SolutionFields& fields) = 0;
/** Get cell field associated with integrator.
*
- * @param fieldType Type of field.
* @param name Name of cell field.
- * @param mesh PETSc mesh for problem.
- * @param fields Fields manager.
+ * @param fields Solution fields.
* @returns Cell field.
*/
virtual
- const ALE::Obj<real_section_type>&
- cellField(VectorFieldEnum* fieldType,
- const char* name,
- const ALE::Obj<Mesh>& mesh,
- topology::FieldsManager* const fields) = 0;
+ const topology::Field<topology::SubMesh>&
+ cellField(const char* name,
+ topology::SolutionFields& fields) = 0;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
- /// Not implemented
- Fault(const Fault& m);
+ Fault(const Fault&); ///< Not implemented
+ const Fault& operator=(const Fault&); ///< Not implemented
- /// Not implemented
- const Fault& operator=(const Fault& m);
-
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
- ALE::Obj<SubMesh> _faultMesh; ///< Mesh over fault surface
+ topology::SubMesh* _faultMesh; ///< Mesh over fault surface
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.icc 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.icc 2009-04-04 23:31:36 UTC (rev 14595)
@@ -37,9 +37,9 @@
// Get label of fault.
inline
-const std::string&
+const char*
pylith::faults::Fault::label(void) const {
- return _label;
+ return _label.c_str();
}
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc 2009-04-04 23:31:36 UTC (rev 14595)
@@ -16,11 +16,7 @@
#include "CohesiveTopology.hh" // USES CohesiveTopology::create()
-#if defined(NEWPYLITHMEHS)
#include "pylith/meshio/UCDFaultFile.hh" // USES UCDFaultFile::read()
-#endif
-
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
#include "pylith/utils/array.hh" // USES double_array
#include <cassert> // USES assert()
@@ -61,33 +57,33 @@
// ----------------------------------------------------------------------
// Adjust mesh topology for fault implementation.
void
-pylith::faults::FaultCohesive::adjustTopology(const ALE::Obj<Mesh>& mesh,
+pylith::faults::FaultCohesive::adjustTopology(const topology::Mesh& mesh,
const bool flipFault)
{ // adjustTopology
assert(std::string("") != label());
+ SubMesh faultMesh;
+ Mesh fault
Obj<SubMesh> faultMesh = NULL;
Obj<ALE::Mesh> faultBd = NULL;
if (_useFaultMesh) {
const int faultDim = 2;
+ assert(3 == mesh.dimension());
//MPI_Bcast(&faultDim, 1, MPI_INT, 0, comm);
faultMesh = new SubMesh(mesh->comm(), faultDim, mesh->debug());
-#if defined(NEWPYLITHMESH)
meshio::UCDFaultFile::readFault(_faultMeshFilename, mesh,
faultMesh, faultBd);
-#endif
// Get group of vertices associated with fault
- const ALE::Obj<int_section_type>& groupField =
- mesh->getIntSection(label());
+ const ALE::Obj<IntSection>& groupField = sieveMesh->getIntSection(label());
faultMesh->setRealSection("coordinates",
mesh->getRealSection("coordinates"));
CohesiveTopology::create(faultMesh, faultBd, mesh, groupField, id(),
_useLagrangeConstraints());
} else {
- if (!mesh->hasIntSection(label())) {
+ if (!sieveMesh->hasIntSection(label())) {
std::ostringstream msg;
msg << "Mesh missing group of vertices '" << label()
<< " for fault interface condition.";
@@ -95,8 +91,7 @@
} // if
// Get group of vertices associated with fault
- const ALE::Obj<int_section_type>& groupField =
- mesh->getIntSection(label());
+ const ALE::Obj<IntSection>& groupField = sieveMesh->getIntSection(label());
assert(!groupField.isNull());
faultMesh =
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.hh 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.hh 2009-04-04 23:31:36 UTC (rev 14595)
@@ -19,21 +19,10 @@
#if !defined(pylith_faults_faultcohesive_hh)
#define pylith_faults_faultcohesive_hh
+// Include directives ---------------------------------------------------
#include "Fault.hh" // ISA Fault
-#include "pylith/utils/sievefwd.hh" // HOLDSA PETSc Mesh
-#include "pylith/utils/petscfwd.h" // USES PetscMat
-
-/// Namespace for pylith package
-namespace pylith {
- namespace faults {
- class FaultCohesive;
- class TestFaultCohesive; // unit testing
- } // faults
-} // pylith
-
-/// C++ abstract base class for a fault surface implemented with
-/// cohesive elements.
+// FaultCohesive --------------------------------------------------------
class pylith::faults::FaultCohesive : public Fault
{ // class FaultCohesive
friend class TestFaultCohesive; // unit testing
@@ -64,9 +53,11 @@
/** Adjust mesh topology for fault implementation.
*
- * @param mesh PETSc mesh
+ * @param mesh PETSc mesh.
+ * @param flipFault Flip fault orientation.
*/
- void adjustTopology(const ALE::Obj<Mesh>& mesh, const bool flipFault = false);
+ void adjustTopology(const topology::Mesh& mesh,
+ const bool flipFault =false);
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/Makefile.am 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/Makefile.am 2009-04-04 23:31:36 UTC (rev 14595)
@@ -33,7 +33,10 @@
StepSlipFn.hh \
StepSlipFn.icc
-noinst_HEADERS =
+noinst_HEADERS = \
+ TopologyOps.hh \
+ TopologyVisitors.hh \
+ TopologyVisitors.cc
# export
clean-local: clean-subpkgincludeHEADERS
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/SlipTimeFn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/SlipTimeFn.hh 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/SlipTimeFn.hh 2009-04-04 23:31:36 UTC (rev 14595)
@@ -20,27 +20,12 @@
#if !defined(pylith_faults_sliptimefn_hh)
#define pylith_faults_sliptimefn_hh
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+// Include directives ---------------------------------------------------
+#include "Fault.hh" // ISA Fault
-/// Namespace for pylith package
-namespace pylith {
- namespace faults {
- class SlipTimeFn;
- class TestSlipTimeFn; // unit testing
- } // faults
-} // pylith
+#include "spatialdata/units/unitsfwd.hh" // USES Nondimensional
-/// Namespace for spatialdata package
-namespace spatialdata {
- namespace geocoords {
- class CoordSys;
- } // geocoords
- namespace units {
- class Nondimensional;
- } // units
-} // spatialdata
-
-/// C++ abstract base class for Fault object.
+// SlipTimeFn -----------------------------------------------------------
class pylith::faults::SlipTimeFn
{ // class SlipTimeFn
friend class TestSlipTimeFn; // unit testing
@@ -63,8 +48,7 @@
* @param originTime Origin time for earthquake source.
*/
virtual
- 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) = 0;
@@ -72,53 +56,46 @@
*
* @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.
*/
virtual
- void slip(const ALE::Obj<real_section_type>& slipField,
- const double t,
- const ALE::Obj<Mesh>& faultMesh) = 0;
+ void slip(topology::Field<topology::SubMesh>* const slipField,
+ const double t) = 0;
/** 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.
*/
virtual
- 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) = 0;
+ const double t1) = 0;
/** Get final slip.
*
* @returns Final slip.
*/
virtual
- ALE::Obj<real_section_type> finalSlip(void) = 0;
+ const topology::Field<topology::SubMesh>& finalSlip(void) = 0;
/** Get time when slip begins at each point.
*
* @returns Time when slip begins.
*/
virtual
- ALE::Obj<real_section_type> slipTime(void) = 0;
+ const topology::Field<topology::SubMesh>& slipTime(void) = 0;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
- /// Not implemented.
- SlipTimeFn(const SlipTimeFn& f);
+ SlipTimeFn(const SlipTimeFn&); ///< Not implemented
+ const SlipTimeFn& operator=(const SlipTimeFn&); ///< Not implemented
- /// Not implemented
- const SlipTimeFn& operator=(const SlipTimeFn& f);
-
}; // class SlipTimeFn
#endif // pylith_faults_sliptimefn_hh
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.cc 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.cc 2009-04-04 23:31:36 UTC (rev 14595)
@@ -14,8 +14,9 @@
#include "StepSlipFn.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,20 +26,18 @@
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
-namespace pylith {
- namespace faults {
- namespace _StepSlipFn {
- const int offsetSlipTime = 0;
- } // _StepSlipFn
- } // 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::StepSlipFn::StepSlipFn(void) :
+ _slipTimeVertex(0),
+ _parameters(0),
_dbFinalSlip(0),
- _dbSlipTime(0),
- _spaceDim(0)
+ _dbSlipTime(0)
{ // constructor
} // constructor
@@ -46,45 +45,56 @@
// Destructor.
pylith::faults::StepSlipFn::~StepSlipFn(void)
{ // destructor
- _dbFinalSlip = 0;
- _dbSlipTime = 0;
+ delete _parameters; _parameters = 0;
+ _dbFinalSlip = 0; // :TODO: Use shared pointer
+ _dbSlipTime = 0; // :TODO: Use shared pointer
} // destructor
// ----------------------------------------------------------------------
// Initialize slip time function.
void
pylith::faults::StepSlipFn::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);
- _spaceDim = cs->spaceDim();
- const int spaceDim = _spaceDim;
- const int indexFinalSlip = 0;
- const int indexSlipTime = spaceDim + _StepSlipFn::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 + 1;
- _parameters = new real_section_type(faultMesh->comm(), faultMesh->debug());
- _parameters->addSpace(); // final slip
- _parameters->addSpace(); // slip time
- assert(2 == _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); // 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());
+
// Open databases and set query values
_dbFinalSlip->open();
switch (spaceDim)
@@ -114,16 +124,13 @@
_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,
@@ -131,7 +138,7 @@
normalizer.dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(),
lengthScale);
- int err = _dbFinalSlip->query(¶msVertex[indexFinalSlip], spaceDim,
+ int err = _dbFinalSlip->query(&_slipVertex[0], _slipVertex.size(),
&vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -141,10 +148,10 @@
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 = _dbSlipTime->query(¶msVertex[indexSlipTime], 1,
+ err = _dbSlipTime->query(&_slipTimeVertex, 1,
&vCoordsGlobal[0], vCoordsGlobal.size(), cs);
if (err) {
std::ostringstream msg;
@@ -154,12 +161,12 @@
msg << ") using spatial database " << _dbSlipTime->label() << ".";
throw std::runtime_error(msg.str());
} // if
- normalizer.nondimensionalize(¶msVertex[indexSlipTime], 1,
- timeScale);
+ normalizer.nondimensionalize(&_slipTimeVertex, 1, timeScale);
// add origin time to rupture time
- paramsVertex[indexSlipTime] += originTime;
+ _slipTimeVertex += originTime;
- _parameters->updatePoint(*v_iter, ¶msVertex[0]);
+ finalSlipSection->updatePoint(*v_iter, &_slipVertex[0]);
+ slipTimeSection->updatePoint(*v_iter, &_slipTimeVertex);
} // for
// Close databases
@@ -170,109 +177,109 @@
// ----------------------------------------------------------------------
// Get slip on fault surface at time t.
void
-pylith::faults::StepSlipFn::slip(const ALE::Obj<pylith::real_section_type>& slipField,
- const double t,
- const ALE::Obj<Mesh>& faultMesh)
+pylith::faults::StepSlipFn::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 indexSlipTime = spaceDim + _StepSlipFn::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 ALE::Obj<RealSection>& slipSection = slip->section();
+ assert(!slipSection.isNull());
+
+ 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);
- const double* finalSlip = ¶msVertex[indexFinalSlip];
- const double slipTime = paramsVertex[indexSlipTime];
- slipValues = 0.0;
-
- const double relTime = t - slipTime;
- if (relTime >= 0.0)
- for (int i=0; i < spaceDim; ++i)
- slipValues[i] = finalSlip[i];
+ const double relTime = t - _slipTimeVertex;
+ if (relTime < 0.0)
+ _slipVertex = 0.0;
// Update field
- slipField->updateAddPoint(*v_iter, &slipValues[0]);
+ slipSection->updateAddPoint(*v_iter, &_slipVertex[0]);
} // for
- PetscLogFlops(numVertices * 1);
+ PetscLogFlops(vertices->size() * 1);
} // slip
// ----------------------------------------------------------------------
// Get increment of slip on fault surface between time t0 and t1.
void
-pylith::faults::StepSlipFn::slipIncr(const ALE::Obj<pylith::real_section_type>& slipField,
- const double t0,
- const double t1,
- const ALE::Obj<Mesh>& faultMesh)
+pylith::faults::StepSlipFn::slipIncr(topology::Field<topology::SubMesh>* slip,
+ const double t0,
+ 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 indexSlipTime = spaceDim + _StepSlipFn::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 ALE::Obj<RealSection>& slipSection = slip->section();
+ assert(!slipSection.isNull());
+
+ 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);
- const double* finalSlip = ¶msVertex[indexFinalSlip];
- const double slipTime = paramsVertex[indexSlipTime];
- slipValues = 0.0;
-
- const double relTime0 = t0 - slipTime;
- const double relTime1 = t1 - slipTime;
- if (relTime1 >= 0.0 && relTime0 < 0.0)
- for (int i=0; i < spaceDim; ++i)
- slipValues[i] = finalSlip[i];
+ const double relTime0 = t0 - _slipTimeVertex;
+ const double relTime1 = t1 - _slipTimeVertex;
+ if (relTime1 < 0.0 || relTime0 >= 0.0)
+ _slipVertex = 0.0;
// Update field
- slipField->updateAddPoint(*v_iter, &slipValues[0]);
+ slipSection->updateAddPoint(*v_iter, &_slipVertex[0]);
} // for
- PetscLogFlops(numVertices * 3);
+ PetscLogFlops(vertices->size() * 2);
} // slipIncr
// ----------------------------------------------------------------------
// Get final slip.
-ALE::Obj<pylith::real_section_type>
+const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::StepSlipFn::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::StepSlipFn::slipTime(void)
{ // slipTime
- return _parameters->getFibration(1);
+ return _parameters->get("slip time");
} // slipTime
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.hh 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.hh 2009-04-04 23:31:36 UTC (rev 14595)
@@ -22,27 +22,14 @@
#if !defined(pylith_faults_stepslipfn_hh)
#define pylith_faults_stepslipfn_hh
+// Include directives ---------------------------------------------------
#include "SlipTimeFn.hh"
-/// Namespace for pylith package
-namespace pylith {
- namespace faults {
- class StepSlipFn;
- class TestStepSlipFn; // 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 Step slip time function.
+// SlipTimeFn -----------------------------------------------------------
class pylith::faults::StepSlipFn : public SlipTimeFn
{ // class StepSlipFn
friend class TestStepSlipFn; // unit testing
@@ -71,12 +58,11 @@
/** Initialize slip time function.
*
* @param faultMesh Finite-element mesh of fault.
- * @param cs Coordinate system for mesh.
+ * @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);
@@ -84,64 +70,58 @@
*
* @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
- StepSlipFn(const StepSlipFn& m);
+ StepSlipFn(const StepSlipFn&); ///< Not implemented.
+ const StepSlipFn& operator=(const StepSlipFn&); ///< Not implemented
- /// Not implemented
- const StepSlipFn& operator=(const StepSlipFn& f);
-
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
- /// Parameters for Step slip time function.
- /// Final slip (vector), slip time (scalar).
- ALE::Obj<real_section_type> _parameters;
+ double _slipTimeVertex; ///< Slip time at a vertex.
+ double_array _slipVertex; ///< Final slip at a vertex.
+ /// Parameters for step slip time function, final slip (vector) and
+ /// slip time (scalar).
+ topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
+
/// Spatial database for final slip
spatialdata::spatialdb::SpatialDB* _dbFinalSlip;
/// Spatial database for slip time
spatialdata::spatialdb::SpatialDB* _dbSlipTime;
- int _spaceDim; ///< Spatial dimension for slip field.
-
}; // class StepSlipFn
#include "StepSlipFn.icc" // inline methods
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.icc 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/StepSlipFn.icc 2009-04-04 23:31:36 UTC (rev 14595)
@@ -14,9 +14,6 @@
#error "StepSlipFn.icc can only be included from StepSlipFn.hh"
#endif
-#include <math.h> // USES exp()
-#include <assert.h> // USES assert()
-
// Set spatial database for final slip.
inline
void
Added: short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.cc 2009-04-04 23:31:36 UTC (rev 14595)
@@ -0,0 +1,510 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TopologyOps.hh" // implementation of object methods
+//#include <Selection.hh> // Algorithms for submeshes
+
+//#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+
+#include <cassert> // USES assert()
+
+// ----------------------------------------------------------------------
+template<class InputPoints>
+bool
+pylith::faults::CohesiveTopology::compatibleOrientation(const ALE::Obj<Mesh>& mesh,
+ const Mesh::point_type& p,
+ const Mesh::point_type& q,
+ const int numFaultCorners,
+ const int faultFaceSize,
+ const int faultDepth,
+ const Obj<InputPoints>& points,
+ int indices[],
+ PointArray *origVertices,
+ PointArray *faceVertices,
+ PointArray *neighborVertices)
+{
+ typedef ALE::Selection<Mesh> selection;
+ const int debug = mesh->debug();
+ bool compatible;
+
+ bool eOrient = selection::getOrientedFace(mesh, p, points, numFaultCorners, indices, origVertices, faceVertices);
+ bool nOrient = selection::getOrientedFace(mesh, q, points, numFaultCorners, indices, origVertices, neighborVertices);
+
+ if (faultFaceSize > 1) {
+ if (debug) {
+ for(PointArray::iterator v_iter = faceVertices->begin(); v_iter != faceVertices->end(); ++v_iter) {
+ std::cout << " face vertex " << *v_iter << std::endl;
+ }
+ for(PointArray::iterator v_iter = neighborVertices->begin(); v_iter != neighborVertices->end(); ++v_iter) {
+ std::cout << " neighbor vertex " << *v_iter << std::endl;
+ }
+ }
+ compatible = !(*faceVertices->begin() == *neighborVertices->begin());
+ } else {
+ compatible = !(nOrient == eOrient);
+ }
+ return compatible;
+}
+
+// ----------------------------------------------------------------------
+void
+pylith::faults::CohesiveTopology::computeCensoredDepth(const ALE::Obj<Mesh::label_type>& depth,
+ const ALE::Obj<Mesh::sieve_type>& sieve,
+ const Mesh::point_type& firstCohesiveCell)
+{
+ Mesh::DepthVisitor d(*sieve, firstCohesiveCell, *depth);
+
+ sieve->roots(d);
+ while(d.isModified()) {
+ // FIX: Avoid the copy here somehow by fixing the traversal
+ std::vector<Mesh::point_type> modifiedPoints(d.getModifiedPoints().begin(), d.getModifiedPoints().end());
+
+ d.clear();
+ sieve->support(modifiedPoints, d);
+ }
+}
+
+// ----------------------------------------------------------------------
+void
+pylith::faults::CohesiveTopology::classifyCells(const ALE::Obj<Mesh::sieve_type>& sieve,
+ const Mesh::point_type& vertex,
+ const int depth,
+ const int faceSize,
+ const Mesh::point_type& firstCohesiveCell,
+ PointSet& replaceCells,
+ PointSet& noReplaceCells,
+ const int debug)
+{
+ // Replace all cells on a given side of the fault with a vertex on the fault
+ ClassifyVisitor<Mesh::sieve_type> cV(*sieve, replaceCells, noReplaceCells, firstCohesiveCell, faceSize, debug);
+ const PointSet& vReplaceCells = cV.getReplaceCells();
+ const PointSet& vNoReplaceCells = cV.getNoReplaceCells();
+
+ if (debug) {std::cout << "Checking fault vertex " << vertex << std::endl;}
+ sieve->support(vertex, cV);
+ cV.setMode(false);
+ const int classifyTotal = cV.getSize();
+ int classifySize = vReplaceCells.size() + vNoReplaceCells.size();
+
+ while(cV.getModified() && (classifySize < classifyTotal)) {
+ cV.reset();
+ sieve->support(vertex, cV);
+ if (debug) {
+ std::cout << "classifySize: " << classifySize << std::endl;
+ std::cout << "classifyTotal: " << classifyTotal << std::endl;
+ std::cout << "vReplaceCells.size: " << vReplaceCells.size() << std::endl;
+ std::cout << "vNoReplaceCells.size: " << vNoReplaceCells.size() << std::endl;
+ }
+ assert(classifySize < vReplaceCells.size() + vNoReplaceCells.size());
+ classifySize = vReplaceCells.size() + vNoReplaceCells.size();
+ assert(classifySize <= classifyTotal);
+ }
+ replaceCells.insert(vReplaceCells.begin(), vReplaceCells.end());
+ // More checking
+ noReplaceCells.insert(vNoReplaceCells.begin(), vNoReplaceCells.end());
+}
+
+// ----------------------------------------------------------------------
+void
+pylith::faults::CohesiveTopology::createFaultSieveFromVertices(const int dim,
+ const int firstCell,
+ const PointSet& faultVertices,
+ const Obj<Mesh>& mesh,
+ const Obj<ALE::Mesh::arrow_section_type>& orientation,
+ const Obj<ALE::Mesh::sieve_type>& faultSieve,
+ const bool flipFault)
+{
+ typedef ALE::Selection<ALE::Mesh> selection;
+ const Obj<sieve_type>& sieve = mesh->getSieve();
+ const PointSet::const_iterator fvBegin = faultVertices.begin();
+ const PointSet::const_iterator fvEnd = faultVertices.end();
+ int curCell = firstCell;
+ int curVertex = 0;
+ int newElement = curCell + dim*faultVertices.size();
+ int o = 1;
+ ALE::Mesh::point_type f = firstCell;
+ const int debug = mesh->debug();
+ Obj<PointSet> face = new PointSet();
+ int numCorners = 0; // The number of vertices in a mesh cell
+ int faceSize = 0; // The number of vertices in a mesh face
+ int *indices = NULL; // The indices of a face vertex set in a cell
+ std::map<int,int*> curElement;
+ std::map<int,PointArray> bdVertices;
+ std::map<int,PointArray> faultFaces;
+ std::map<int,oPointArray> oFaultFaces;
+ PointSet faultCells;
+ PointArray origVertices;
+ PointArray faceVertices;
+
+ if (!faultSieve->commRank()) {
+ numCorners = mesh->getNumCellCorners();
+ faceSize = selection::numFaceVertices(mesh);
+ indices = new int[faceSize];
+ }
+
+ curElement[0] = &curVertex;
+ curElement[dim] = &curCell;
+ for(int d = 1; d < dim; d++) {
+ curElement[d] = &newElement;
+ }
+
+ // This only works for uninterpolated meshes
+ assert((mesh->depth() == 1) || (mesh->depth() == -1));
+ ALE::ISieveVisitor::PointRetriever<sieve_type> sV(std::max(1, sieve->getMaxSupportSize()));
+ ALE::ISieveVisitor::PointRetriever<sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
+ for(PointSet::const_iterator fv_iter = fvBegin; fv_iter != fvEnd; ++fv_iter) {
+ sieve->support(*fv_iter, sV);
+ const Mesh::point_type *support = sV.getPoints();
+
+ if (debug) std::cout << "Checking fault vertex " << *fv_iter << std::endl;
+ const int sVsize = sV.getSize();
+ for (int i=0; i < sVsize; ++i) {
+ const int s = (!flipFault) ? i : sVsize - i - 1;
+ sieve->cone(support[s], cV);
+ const Mesh::point_type *cone = cV.getPoints();
+
+ if (debug) std::cout << " Checking cell " << support[s] << std::endl;
+ if (faultCells.find(support[s]) != faultCells.end()) {
+ cV.clear();
+ continue;
+ }
+ face->clear();
+ for(int c = 0; c < cV.getSize(); ++c) {
+ if (faultVertices.find(cone[c]) != fvEnd) {
+ if (debug) std::cout << " contains fault vertex " << cone[c] << std::endl;
+ face->insert(face->end(), cone[c]);
+ } // if
+ } // for
+ if (face->size() > faceSize)
+ throw ALE::Exception("Invalid fault mesh: Too many vertices of an "
+ "element on the fault");
+ if (face->size() == faceSize) {
+ if (debug) std::cout << " Contains a face on the fault" << std::endl;
+ ALE::Obj<sieve_type::supportSet> preFace;
+ if (dim < 2) {
+ preFace = faultSieve->nJoin1(face);
+ } else {
+ preFace = faultSieve->nJoin(face, dim);
+ }
+
+ if (preFace->size() > 1) {
+ throw ALE::Exception("Invalid fault sieve: Multiple faces from vertex set");
+ } else if (preFace->size() == 1) {
+ // Add the other cell neighbor for this face
+ if (dim == 0) {
+ faultSieve->addArrow(*faceVertices.begin(), support[s]);
+ } else {
+ faultSieve->addArrow(*preFace->begin(), support[s]);
+ }
+ } else if (preFace->size() == 0) {
+ if (debug) std::cout << " Orienting face " << f << std::endl;
+ selection::getOrientedFace(mesh, support[s], face, numCorners, indices, &origVertices, &faceVertices);
+ bdVertices[dim].clear();
+ for(PointArray::const_iterator v_iter = faceVertices.begin(); v_iter != faceVertices.end(); ++v_iter) {
+ bdVertices[dim].push_back(*v_iter);
+ if (debug) std::cout << " Boundary vertex " << *v_iter << std::endl;
+ }
+ if (dim == 0) {
+ f = *faceVertices.begin();
+ }
+ if (faceSize != dim+1) {
+ if (debug) std::cout << " Adding hex face " << f << std::endl;
+ ALE::SieveBuilder<ALE::Mesh>::buildHexFaces(faultSieve, orientation, dim, curElement, bdVertices, oFaultFaces, f, o);
+ } else {
+ if (debug) std::cout << " Adding simplicial face " << f << std::endl;
+ ALE::SieveBuilder<ALE::Mesh>::buildFaces(faultSieve, orientation, dim, curElement, bdVertices, oFaultFaces, f, o);
+ }
+ faultSieve->addArrow(f, support[s]);
+ //faultSieve->view("");
+ f++;
+ } // if/else
+ faultCells.insert(support[s]);
+ } // if
+ cV.clear();
+ } // for
+ sV.clear();
+ } // for
+ if (!faultSieve->commRank()) delete [] indices;
+}
+
+// ----------------------------------------------------------------------
+void
+pylith::faults::CohesiveTopology::createFaultSieveFromFaces(const int dim,
+ const int firstCell,
+ const int numFaces,
+ const int faultVertices[],
+ const int faultCells[],
+ const Obj<Mesh>& mesh,
+ const Obj<ALE::Mesh::arrow_section_type>& orientation,
+ const Obj<ALE::Mesh::sieve_type>& faultSieve)
+{
+ typedef ALE::Selection<ALE::Mesh> selection;
+ int faceSize = 0; // The number of vertices in a mesh face
+ int curCell = firstCell;
+ int curVertex = 0;
+ int newElement = curCell + dim*numFaces;
+ int o = 1;
+ int f = firstCell;
+ const int debug = mesh->debug();
+ std::map<int,int*> curElement;
+ std::map<int,PointArray> bdVertices;
+ std::map<int,oPointArray> oFaultFaces;
+
+ if (!faultSieve->commRank()) {
+ faceSize = selection::numFaceVertices(mesh);
+ }
+
+ curElement[0] = &curVertex;
+ curElement[dim] = &curCell;
+ for(int d = 1; d < dim; d++) {
+ curElement[d] = &newElement;
+ }
+
+ // Loop over fault faces
+ for(int face = 0; face < numFaces; ++face) {
+ // Push oriented vertices of face
+ bdVertices[dim].clear();
+ for(int i = 0; i < faceSize; ++i) {
+ bdVertices[dim].push_back(faultVertices[face*faceSize+i]);
+ if (debug) std::cout << " Boundary vertex " << faultVertices[face*faceSize+i] << std::endl;
+ }
+ // Create face
+ if (faceSize != dim+1) {
+ if (debug) std::cout << " Adding hex face " << f << std::endl;
+ ALE::SieveBuilder<ALE::Mesh>::buildHexFaces(faultSieve, orientation, dim, curElement, bdVertices, oFaultFaces, f, o);
+ } else {
+ if (debug) std::cout << " Adding simplicial face " << f << std::endl;
+ ALE::SieveBuilder<ALE::Mesh>::buildFaces(faultSieve, orientation, dim, curElement, bdVertices, oFaultFaces, f, o);
+ }
+ // Add arrow to cells
+ faultSieve->addArrow(face, faultCells[face*2+0]);
+ faultSieve->addArrow(face, faultCells[face*2+1]);
+ }
+}
+
+// ----------------------------------------------------------------------
+void
+pylith::faults::CohesiveTopology::orientFaultSieve(const int dim,
+ const Obj<Mesh>& mesh,
+ const Obj<ALE::Mesh::arrow_section_type>& orientation,
+ const Obj<ALE::Mesh>& fault)
+{
+ // Must check the orientation here
+ typedef ALE::Selection<ALE::Mesh> selection;
+ const Obj<ALE::Mesh::sieve_type>& faultSieve = fault->getSieve();
+ const Mesh::point_type firstFaultCell = *fault->heightStratum(1)->begin();
+ const Obj<ALE::Mesh::label_sequence>& fFaces = fault->heightStratum(2);
+ const int numFaultFaces = fFaces->size();
+ const int faultDepth = fault->depth()-1; // Depth of fault cells
+ int numFaultCorners = 0; // The number of vertices in a fault cell
+ int faultFaceSize = 0; // The number of vertices in a face between fault cells
+ int faceSize = 0; // The number of vertices in a mesh face
+ const int debug = fault->debug();
+ Obj<PointSet> newCells = new PointSet();
+ Obj<PointSet> loopCells = new PointSet();
+ PointSet flippedCells; // Incorrectly oriented fault cells
+ PointSet facesSeen; // Fault faces already considered
+ PointSet cellsSeen; // Fault cells already matched
+ PointArray faceVertices;
+
+ if (!fault->commRank()) {
+ faceSize = selection::numFaceVertices(mesh);
+ numFaultCorners = faultSieve->nCone(firstFaultCell, faultDepth)->size();
+ if (debug) std::cout << " Fault corners " << numFaultCorners << std::endl;
+ if (dim == 0) {
+ assert(numFaultCorners == faceSize-1);
+ } else {
+ assert(numFaultCorners == faceSize);
+ }
+ if (faultDepth == 1) {
+ faultFaceSize = 1;
+ } else {
+ faultFaceSize = faultSieve->nCone(*fFaces->begin(), faultDepth-1)->size();
+ }
+ }
+ if (debug) std::cout << " Fault face size " << faultFaceSize << std::endl;
+
+ newCells->insert(firstFaultCell);
+ while(facesSeen.size() != numFaultFaces) {
+ Obj<PointSet> tmp = newCells; newCells = loopCells; loopCells = tmp;
+
+ newCells->clear();
+ if (!loopCells->size()) {throw ALE::Exception("Fault surface not a single connected component.");}
+ // Loop over new cells
+ for(PointSet::iterator c_iter = loopCells->begin(); c_iter != loopCells->end(); ++c_iter) {
+ // Loop over edges of this cell
+ const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& cone = faultSieve->cone(*c_iter);
+ const ALE::Mesh::sieve_type::traits::coneSequence::iterator eBegin = cone->begin();
+ const ALE::Mesh::sieve_type::traits::coneSequence::iterator eEnd = cone->end();
+
+ for(ALE::Mesh::sieve_type::traits::coneSequence::iterator e_iter = eBegin; e_iter != eEnd; ++e_iter) {
+ if (facesSeen.find(*e_iter) != facesSeen.end()) continue;
+ facesSeen.insert(*e_iter);
+ if (debug) std::cout << " Checking orientation of fault face " << *e_iter << std::endl;
+ const Obj<ALE::Mesh::sieve_type::traits::supportSequence>& support = faultSieve->support(*e_iter);
+ ALE::Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin();
+
+ // Throw out boundary fault faces
+ if (support->size() < 2) continue;
+ ALE::Mesh::point_type cellA = *s_iter; ++s_iter;
+ ALE::Mesh::point_type cellB = *s_iter;
+ bool flippedA = (flippedCells.find(cellA) != flippedCells.end());
+ bool flippedB = (flippedCells.find(cellB) != flippedCells.end());
+ bool seenA = (cellsSeen.find(cellA) != cellsSeen.end());
+ bool seenB = (cellsSeen.find(cellB) != cellsSeen.end());
+
+ if (!seenA) newCells->insert(cellA);
+ if (!seenB) newCells->insert(cellB);
+ if (debug) std::cout << " neighboring cells " << cellA << " and " << cellB << std::endl;
+ // In 1D, just check that vertices match
+ if (dim == 1) {
+ const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneA = faultSieve->cone(cellA);
+ ALE::Mesh::sieve_type::traits::coneSequence::iterator iterA = coneA->begin();
+ const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneB = faultSieve->cone(cellB);
+ ALE::Mesh::sieve_type::traits::coneSequence::iterator iterB = coneB->begin();
+ int posA, posB;
+
+ for(posA = 0; posA < 2; ++posA, ++iterA) if (*iterA == *e_iter) break;
+ for(posB = 0; posB < 2; ++posB, ++iterB) if (*iterB == *e_iter) break;
+ if (debug) std::cout << " with face positions " << posA << " and " << posB << std::endl;
+ if ((posA == 2) || (posB == 2)) {throw ALE::Exception("Could not find fault face in cone");}
+ if ((posA == posB) ^ (flippedA || flippedB)) {
+ if (debug) {
+ std::cout << "Invalid orientation in fault mesh" << std::endl;
+ std::cout << " fault face: " << *e_iter << " cellA: " << cellA << " cellB: " << cellB << std::endl;
+ }
+ if (flippedA && flippedB) {throw ALE::Exception("Attempt to flip already flipped cell: Fault mesh is non-orientable");}
+ if (seenA && seenB) {throw ALE::Exception("Previous cells do not match: Fault mesh is non-orientable");}
+ if (!seenA && !flippedA) {
+ flippedCells.insert(cellA);
+ } else if (!seenB && !flippedB) {
+ flippedCells.insert(cellB);
+ } else {
+ throw ALE::Exception("Inconsistent mesh orientation: Fault mesh is non-orientable");
+ }
+ }
+ } else if (dim == 2) {
+ // Check orientation
+ ALE::MinimalArrow<ALE::Mesh::sieve_type::point_type,ALE::Mesh::sieve_type::point_type> arrowA(*e_iter, cellA);
+ const int oA = orientation->restrictPoint(arrowA)[0];
+ ALE::MinimalArrow<ALE::Mesh::sieve_type::point_type,ALE::Mesh::sieve_type::point_type> arrowB(*e_iter, cellB);
+ const int oB = orientation->restrictPoint(arrowB)[0];
+ const bool mismatch = (oA == oB);
+
+ // Truth Table
+ // mismatch flips action mismatch flipA ^ flipB action
+ // F 0 flips no F F F
+ // F 1 flip yes F T T
+ // F 2 flips no T F T
+ // T 0 flips yes T T F
+ // T 1 flip no
+ // T 2 flips yes
+ if (mismatch ^ (flippedA ^ flippedB)) {
+ if (debug) {
+ std::cout << "Invalid orientation in fault mesh" << std::endl;
+ std::cout << " fault face: " << *e_iter << " cellA: " << cellA << " cellB: " << cellB << std::endl;
+ }
+ if (flippedA && flippedB) {throw ALE::Exception("Attempt to flip already flipped cell: Fault mesh is non-orientable");}
+ if (seenA && seenB) {throw ALE::Exception("Previous cells do not match: Fault mesh is non-orientable");}
+ if (!seenA && !flippedA) {
+ flippedCells.insert(cellA);
+ if (debug) {std::cout << " Scheduling cell " << cellA << " for flipping" << std::endl;}
+ } else if (!seenB && !flippedB) {
+ flippedCells.insert(cellB);
+ if (debug) {std::cout << " Scheduling cell " << cellB << " for flipping" << std::endl;}
+ } else {
+ throw ALE::Exception("Inconsistent mesh orientation: Fault mesh is non-orientable");
+ }
+ }
+ }
+ cellsSeen.insert(cellA);
+ cellsSeen.insert(cellB);
+ }
+ }
+ }
+ for(PointSet::const_iterator f_iter = flippedCells.begin(); f_iter != flippedCells.end(); ++f_iter) {
+ if (debug) std::cout << " Reversing fault face " << *f_iter << std::endl;
+ faceVertices.clear();
+ const ALE::Obj<ALE::Mesh::sieve_type::traits::coneSequence>& cone = faultSieve->cone(*f_iter);
+ for(ALE::Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
+ faceVertices.insert(faceVertices.begin(), *v_iter);
+ }
+ faultSieve->clearCone(*f_iter);
+ int color = 0;
+ for(PointArray::const_iterator v_iter = faceVertices.begin(); v_iter != faceVertices.end(); ++v_iter) {
+ faultSieve->addArrow(*v_iter, *f_iter, color++);
+ }
+
+ if (dim > 1) {
+ // Here, they are edges, not vertices
+ for(PointArray::const_iterator e_iter = faceVertices.begin(); e_iter != faceVertices.end(); ++e_iter) {
+ ALE::MinimalArrow<ALE::Mesh::sieve_type::point_type,ALE::Mesh::sieve_type::point_type> arrow(*e_iter, *f_iter);
+ int o = orientation->restrictPoint(arrow)[0];
+
+ if (debug) std::cout << " Reversing orientation of " << *e_iter <<"-->"<<*f_iter << " from " << o << " to " << -(o+1) << std::endl;
+ o = -(o+1);
+ orientation->updatePoint(arrow, &o);
+ }
+ }
+ }
+ flippedCells.clear();
+ for(ALE::Mesh::label_sequence::iterator e_iter = fFaces->begin(); e_iter != fFaces->end(); ++e_iter) {
+ if (debug) std::cout << " Checking orientation of fault face " << *e_iter << std::endl;
+ // for each face get the support (2 fault cells)
+ const Obj<ALE::Mesh::sieve_type::traits::supportSequence>& support = faultSieve->support(*e_iter);
+ ALE::Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin();
+
+ // Throw out boundary fault faces
+ if (support->size() > 1) {
+ ALE::Mesh::point_type cellA = *s_iter; ++s_iter;
+ ALE::Mesh::point_type cellB = *s_iter;
+
+ if (debug) std::cout << " neighboring cells " << cellA << " and " << cellB << std::endl;
+ // In 1D, just check that vertices match
+ if (dim == 1) {
+ const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneA = faultSieve->cone(cellA);
+ ALE::Mesh::sieve_type::traits::coneSequence::iterator iterA = coneA->begin();
+ const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneB = faultSieve->cone(cellB);
+ ALE::Mesh::sieve_type::traits::coneSequence::iterator iterB = coneB->begin();
+ int posA, posB;
+
+ for(posA = 0; posA < 2; ++posA, ++iterA) if (*iterA == *e_iter) break;
+ for(posB = 0; posB < 2; ++posB, ++iterB) if (*iterB == *e_iter) break;
+ if (debug) std::cout << " with face positions " << posA << " and " << posB << std::endl;
+ if ((posA == 2) || (posB == 2)) {throw ALE::Exception("Could not find fault face in cone");}
+ if (posA == posB) {
+ std::cout << "Invalid orientation in fault mesh" << std::endl;
+ std::cout << " fault face: " << *e_iter << " cellA: " << cellA << " cellB: " << cellB << std::endl;
+ throw ALE::Exception("Invalid orientation in fault mesh");
+ }
+ } else {
+ // Check orientation
+ ALE::MinimalArrow<ALE::Mesh::sieve_type::point_type,ALE::Mesh::sieve_type::point_type> arrowA(*e_iter, cellA);
+ const int oA = orientation->restrictPoint(arrowA)[0];
+ ALE::MinimalArrow<ALE::Mesh::sieve_type::point_type,ALE::Mesh::sieve_type::point_type> arrowB(*e_iter, cellB);
+ const int oB = orientation->restrictPoint(arrowB)[0];
+
+ if (oA == oB) {
+ std::cout << "Invalid orientation in fault mesh" << std::endl;
+ std::cout << " fault face: " << *e_iter << " cellA: " << cellA << " cellB: " << cellB << std::endl;
+ throw ALE::Exception("Invalid orientation in fault mesh");
+ }
+ }
+ }
+ }
+ if (debug) fault->view("Oriented Fault mesh");
+}
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.hh 2009-04-04 23:31:36 UTC (rev 14595)
@@ -0,0 +1,94 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/faults/TopologyOps.hh
+ *
+ * @brief C++ object to manage creation of cohesive cells.
+ */
+
+#if !defined(pylith_faults_topologyops_hh)
+#define pylith_faults_topologyops_hh
+
+// Include directives ---------------------------------------------------
+#include "faultsfwd.hh" // forward declarations
+
+// TopologyOps ----------------------------------------------------------
+class pylith::faults::TopologyOps
+{ // class TopologyOps
+public :
+ typedef std::set<SieveMesh::point_type> PointSet;
+ typedef std::vector<sieve_type::point_type> PointArray;
+ typedef std::pair<sieve_type::point_type, int> oPoint_type;
+ typedef std::vector<oPoint_type> oPointArray;
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ template<class InputPoints>
+ static
+ bool compatibleOrientation(const topology::& mesh,
+ const SieveMesh::point_type& p,
+ const SieveMesh::point_type& q,
+ const int numFaultCorners,
+ const int faultFaceSize,
+ const int faultDepth,
+ const ALE::Obj<InputPoints>& points,
+ int indices[],
+ PointArray *origVertices,
+ PointArray *faceVertices,
+ PointArray *neighborVertices);
+
+ static
+ void computeCensoredDepth(const ALE::Obj<Mesh::label_type>& depth,
+ const ALE::Obj<Mesh::sieve_type>& sieve,
+ const Mesh::point_type& firstCohesiveCell);
+
+ static
+ void classifyCells(const ALE::Obj<Mesh::sieve_type>& sieve,
+ const Mesh::point_type& vertex,
+ const int depth,
+ const int faceSize,
+ const Mesh::point_type& firstCohesiveCell,
+ PointSet& replaceCells,
+ PointSet& noReplaceCells,
+ const int debug);
+
+ static
+ void createFaultSieveFromVertices(const int dim,
+ const int firstCell,
+ const PointSet& faultVertices,
+ const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<ALE::Mesh::arrow_section_type>& orientation,
+ const ALE::Obj<ALE::Mesh::sieve_type>& faultSieve,
+ const bool flipFault);
+
+ static
+ void createFaultSieveFromFaces(const int dim,
+ const int firstCell,
+ const int numFaces,
+ const int faultVertices[],
+ const int faultCells[],
+ const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<ALE::Mesh::arrow_section_type>& orientation,
+ const ALE::Obj<ALE::Mesh::sieve_type>& faultSieve);
+
+ static
+ void orientFaultSieve(const int dim,
+ const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<ALE::Mesh::arrow_section_type>& orientation,
+ const ALE::Obj<ALE::Mesh>& fault);
+}; // class CohesiveTopology
+
+#endif // pylith_faults_cohesivetopology_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.cc 2009-04-04 23:31:36 UTC (rev 14595)
@@ -0,0 +1,257 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+// ----------------------------------------------------------------------
+template<typename Sieve, typename Renumbering>
+pylith::faults::ReplaceVisitor<Sieve,Renumbering>::ReplaceVisitor(
+ Renumbering& r,
+ const int size,
+ const int debug = 0) :
+ renumbering(r),
+ size(size),
+ i(0),
+ debug(debug)
+{ // constructor
+ this->points = new point_type[this->size];
+ this->mapped = false;
+} // constructor
+
+// ----------------------------------------------------------------------
+template<typename Sieve, typename Renumbering>
+pylith::faults::ReplaceVisitor<Sieve,Renumbering>::~ReplaceVisitor()
+{ // destructor
+ delete[] this->points;
+} // destructor
+
+// ----------------------------------------------------------------------
+template<typename Sieve, typename Renumbering>
+void
+pylith::faults::ReplaceVisitor<Sieve,Renumbering>::visitPoint(
+ const point_type& point)
+{ // visitPoint
+ if (i >= this->size)
+ throw ALE::Exception("Too many points for ReplaceVisitor");
+ if (this->renumbering.find(point) != this->renumbering.end()) {
+ if (debug)
+ std::cout << " point " << this->renumbering[point] << std::endl;
+ points[i] = this->renumbering[point];
+ this->mapped = true;
+ } else {
+ if (debug) std::cout << " point " << point << std::endl;
+ points[i] = point;
+ } // if/else
+ ++i;
+} // visitPoint
+
+// ----------------------------------------------------------------------
+template<typename Sieve, typename Renumbering>
+void
+pylith::faults::ReplaceVisitor<Sieve,Renumbering>::visitArrow(
+ const typename Sieve::arrow_type&)
+{ // visitArrow
+} // visitArrow
+
+// ----------------------------------------------------------------------
+template<typename Sieve, typename Renumbering>
+inline
+const point_type*
+pylith::faults::ReplaceVisitor<Sieve,Renumbering>::getPoints(void)
+{ // getPoints
+ return this->points;
+} // getPoints
+
+// ----------------------------------------------------------------------
+template<typename Sieve, typename Renumbering>
+inline
+bool
+pylith::faults::ReplaceVisitor<Sieve,Renumbering>::mappedPoint(void)
+{ // mappedPoint
+ return this->mapped;
+} // mappedPoint
+
+// ----------------------------------------------------------------------
+template<typename Sieve, typename Renumbering>
+inline
+void
+pylith::faults::ReplaceVisitor<Sieve,Renumbering>::clear(void)
+{ // clear
+ this->i = 0; this->mapped = false;
+} // clear
+
+
+// ClassifyVisitor ------------------------------------------------------
+template<typename Sieve>
+pylith::faults::ClassifyVisitor<Sieve>::ClassifyVisitor(const Sieve& s,
+ const PointSet& rC,
+ const PointSet& nrC,
+ const point_type& fC,
+ const int fS,
+ const int debug =0) :
+ sieve(s),
+ replaceCells(rC),
+ noReplaceCells(nrC),
+ firstCohesiveCell(fC),
+ faceSize(fS),
+ debug(debug),
+ modified(false),
+ setupMode(true),
+ size(0)
+{ // constructor
+ pR.setSize(s.getMaxConeSize());
+} // constructor
+
+// ----------------------------------------------------------------------
+template<typename Sieve>
+pylith::faults::ClassifyVisitor<Sieve>::~ClassifyVisitor(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+template<typename Sieve>
+void
+pylith::faults::ClassifyVisitor<Sieve>::visitPoint(const point_type& point)
+{ // visitPoint
+ if (this->setupMode) {
+ if (replaceCells.find(point) != replaceCells.end())
+ vReplaceCells.insert(point);
+ if (noReplaceCells.find(point) != noReplaceCells.end())
+ vNoReplaceCells.insert(point);
+ if (point >= firstCohesiveCell) return;
+ this->modified = true;
+ this->size++;
+ return;
+ } // if
+ bool classified = false;
+
+ if (debug)
+ std::cout << "Checking neighbor " << point << std::endl;
+ if (vReplaceCells.find(point) != vReplaceCells.end()) {
+ if (debug)
+ std::cout << " already in replaceCells" << std::endl;
+ return;
+ } // if
+ if (vNoReplaceCells.find(point) != vNoReplaceCells.end()) {
+ if (debug)
+ std::cout << " already in noReplaceCells" << std::endl;
+ return;
+ } // if
+ if (point >= firstCohesiveCell) {
+ if (debug)
+ std::cout << " already a cohesive cell" << std::endl;
+ return;
+ } // if
+ // If neighbor shares a face with anyone in replaceCells, then add
+ for(PointSet::const_iterator c_iter = vReplaceCells.begin();
+ c_iter != vReplaceCells.end();
+ ++c_iter) {
+ sieve.meet(*c_iter, point, pR);
+
+ if (pR.getSize() == faceSize) {
+ if (debug)
+ std::cout << " Scheduling " << point << " for replacement"
+ << std::endl;
+ vReplaceCells.insert(point);
+ modified = true;
+ classified = true;
+ pR.clear();
+ break;
+ } // if
+ pR.clear();
+ } // for
+ if (classified)
+ return;
+ // It is unclear whether taking out the noReplace cells will speed this up
+ for(PointSet::const_iterator c_iter = vNoReplaceCells.begin();
+ c_iter != vNoReplaceCells.end();
+ ++c_iter) {
+ sieve.meet(*c_iter, point, pR);
+
+ if (pR.getSize() == faceSize) {
+ if (debug)
+ std::cout << " Scheduling " << point << " for no replacement"
+ << std::endl;
+ vNoReplaceCells.insert(point);
+ modified = true;
+ classified = true;
+ pR.clear();
+ break;
+ } // for
+ pR.clear();
+ } // for
+} // visitPoint
+
+// ----------------------------------------------------------------------
+template<typename Sieve>
+inline
+void
+pylith::faults::ClassifyVisitor<Sieve>::visitArrow(const typename Sieve::arrow_type&)
+{ // visitArrow
+} // visitArrow
+
+// ----------------------------------------------------------------------
+template<typename Sieve>
+inline
+const PointSet&
+pylith::faults::ClassifyVisitor<Sieve>::getReplaceCells(void) const
+{ // getReplaceCells
+ return this->vReplaceCells;
+} // getReplaceCells
+
+// ----------------------------------------------------------------------
+template<typename Sieve>
+inline
+const PointSet&
+pylith::faults::ClassifyVisitor<Sieve>::getNoReplaceCells() const
+{ // getNoReplaceCells
+ return this->vNoReplaceCells;
+} // getNoReplaceCells
+
+// ----------------------------------------------------------------------
+template<typename Sieve>
+inline
+bool
+pylith::faults::ClassifyVisitor<Sieve>::getModified() const
+{ // getModified
+ return this->modified;
+} // getModified
+
+// ----------------------------------------------------------------------
+template<typename Sieve>
+inline
+int
+pylith::faults::ClassifyVisitor<Sieve>::getSize() const
+{ // getSize
+ return this->size;
+} // getSize
+
+// ----------------------------------------------------------------------
+template<typename Sieve>
+inline
+void
+pylith::faults::ClassifyVisitor<Sieve>::setMode(const bool isSetup)
+{ // setMode
+ this->setupMode = isSetup;
+} // setMode
+
+// ----------------------------------------------------------------------
+template<typename Sieve>
+inline
+void
+pylith::faults::ClassifyVisitor<Sieve>::reset(void)
+{ // reset
+ this->modified = false;
+} // reset
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.hh 2009-04-04 23:31:36 UTC (rev 14595)
@@ -0,0 +1,88 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/faults/TopologyVisitors.hh
+ *
+ * @brief C++ objects for visitors to topology data structures.
+ */
+
+#if !defined(pylith_faults_topologyvisitors_hh)
+#define pylith_faults_topologyvisitors_hh
+
+// Include directives ---------------------------------------------------
+#include "faultsfwd.hh" // forward declarations
+
+// ReplaceVisitor -------------------------------------------------------
+template<typename Sieve, typename Renumbering>
+class pylith::faults::ReplaceVisitor {
+public:
+ typedef typename Sieve::point_type point_type;
+protected:
+ Renumbering& renumbering;
+ const int size;
+ int i;
+ const int debug;
+ point_type *points;
+ bool mapped;
+public:
+ ReplaceVisitor(Renumbering& r,
+ const int size,
+ const int debug =0);
+ ~ReplaceVisitor(void);
+ void visitPoint(const point_type& point);
+ void visitArrow(const typename Sieve::arrow_type&);
+ const point_type *getPoints(void);
+ bool mappedPoint(void);
+ void clear(void);
+};
+
+// ClassifyVisitor ------------------------------------------------------
+template<typename Sieve>
+class pylith::faults::ClassifyVisitor {
+public:
+ typedef typename Sieve::point_type point_type;
+protected:
+ const Sieve& sieve;
+ const PointSet& replaceCells;
+ const PointSet& noReplaceCells;
+ const point_type firstCohesiveCell;
+ const int faceSize;
+ const int debug;
+ PointSet vReplaceCells;
+ PointSet vNoReplaceCells;
+ bool modified;
+ bool setupMode;
+ int size;
+ ALE::ISieveVisitor::PointRetriever<Sieve> pR;
+public:
+ ClassifyVisitor(const Sieve& s,
+ const PointSet& rC,
+ const PointSet& nrC,
+ const point_type& fC,
+ const int fS,
+ const int debug =0);
+ ~ClassifyVisitor(void);
+ void visitPoint(const point_type& point);
+ void visitArrow(const typename Sieve::arrow_type&);
+ const PointSet& getReplaceCells() const;
+ const PointSet& getNoReplaceCells() const;
+ const bool getModified() const;
+ const int getSize() const;
+ void setMode(const bool isSetup);
+ void reset();
+};
+
+
+#endif // pylith_faults_topologyvisitors_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/faults/faultsfwd.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/faultsfwd.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/faultsfwd.hh 2009-04-04 23:31:36 UTC (rev 14595)
@@ -0,0 +1,51 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/** @file libsrc/faults/faultsfwd.hh
+ *
+ * @brief Forward declarations for PyLith faults objects.
+ *
+ * Including this header file eliminates the need to use separate
+ * forward declarations.
+ */
+
+#if !defined(pylith_faults_faultsfwd_hh)
+#define pylith_faults_faultsfwd_hh
+
+namespace pylith {
+ namespace faults {
+
+ class CohesiveTopology;
+
+ class Fault;
+ class FaultCohesiveDyn;
+ class FaultCohesiveKin;
+
+ class EqSinSrc;
+ class SlipTimeFn;
+ class BrundSlipFn;
+ class ConstRateSlipFn;
+ class LiuCosSlipFn;
+ class StepSlipFn;
+
+ class TopologyOps;
+ template<typename Sieve, typename Renumbering> class ReplaceVisitor;
+ template<typename Sieve> class ClassifyVisitor;
+
+ } // faults
+} // pylith
+
+
+#endif // pylith_faults_bcfwd_hh
+
+
+// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Fields.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Fields.hh 2009-04-04 20:23:37 UTC (rev 14594)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Fields.hh 2009-04-04 23:31:36 UTC (rev 14595)
@@ -22,6 +22,8 @@
// Include directives ---------------------------------------------------
#include "topologyfwd.hh" // forward declarations
+#include "pylith/topology/FieldBase.hh" // USES FieldBase::DomainEnum
+
#include <string> // USES std::string
// Fields ---------------------------------------------------------------
More information about the CIG-COMMITS
mailing list