[cig-commits] r6835 - in short/3D/PyLith/trunk: . libsrc
libsrc/faults unittests/libtests/faults
brad at geodynamics.org
brad at geodynamics.org
Wed May 9 18:34:20 PDT 2007
Author: brad
Date: 2007-05-09 18:34:20 -0700 (Wed, 09 May 2007)
New Revision: 6835
Added:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/trunk/libsrc/faults/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/faults/TestFault.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFault.hh
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.hh
Log:
Worked on unit tests for FaultCohesive.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/TODO 2007-05-10 01:34:20 UTC (rev 6835)
@@ -20,26 +20,63 @@
Use MeshIOAscii to get mesh information into Python tests?
0. Implement faults for kinematic source
- a. Add creation of cohesive cells
- i. Add tests for interpolated meshes and constrained cohesive cells.
+ a. Creation of cohesive cells
+ i. Add tests for interpolated meshes.
Double check consistency in ordering of vertices (positive/negative).
- b. Implement slip time function
- ii. Add Python
- iii. Add unit tests
+ b. Implement integrator for faults
- c. Start implementing integrator for faults
+ i. FaultCohesive
+ (1) C++ Unit tests
+ _orient1D()
+ _orient2D()
+ _orient3D()
+ (2) Python unit tests
+ constructor
+ initialize()
- i. Update Reference cell
- (1) Improve how we get Jacobian at vertices of cell
+ ii. EqKinSrc
+ (1) C++ unit tests
+ constructor
+ clone()
+ slipfn()
+ initialize()
+ slip()
+ (2) Python unit tests
+ constructor
+ initialize()
- iii. FaultCohesive
- (1) Unit tests
+ iii. BruneSlipFn
+ (1) C++ unit tests
+ constructor
+ clone()
+ dbFinalSlip()
+ dbSlipTime()
+ dbPeakRate()
+ initialize()
+ slip()
+ _slip()
+ (2) Python unit tests
+ constructor
+ initialize()
iv. FaultCohesiveKin
- (2) Add Python unit tests
- (3) Add C++ unit tests
+ (1) C++ unit tests
+ constructor
+ clone()
+ eqsrc()
+ initialize()
+ integrateResidual()
+ integrateJacobian()
+ setField()
+ _useLagrangeConstraint()
+ (2) Python unit tests
+ constructor
+ initialize()
+ v. Update Reference cell
+ (1) Improve how we get Jacobian at vertices of cell
+
1. Finish implementing ExplicitElasticity and Explicit
a. Replace integrateConstant() with integrateResidual()
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2007-05-10 01:34:20 UTC (rev 6835)
@@ -27,6 +27,7 @@
faults/Fault.cc \
faults/FaultCohesive.cc \
faults/FaultCohesiveKin.cc \
+ faults/FaultCohesiveDyn.cc \
faults/SlipTimeFn.cc \
feassemble/ExplicitElasticity.cc \
feassemble/ImplicitElasticity.cc \
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2007-05-10 01:34:20 UTC (rev 6835)
@@ -16,6 +16,8 @@
#include "CohesiveTopology.hh" // USES CohesiveTopology::create()
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+
#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
#include "pylith/utils/array.hh" // USES double_array
@@ -170,5 +172,16 @@
} // for
} // _orient3D
+// ----------------------------------------------------------------------
+// Get size (fiber dimension) of orientation information.
+int
+pylith::faults::FaultCohesive::_orientationSize(void) const
+{ // _orientationSize
+ assert(0 != _quadrature);
+ const int size = _quadrature->cellDim()*_quadrature->spaceDim();
+ return size;
+} // _orientationSize
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2007-05-10 01:34:20 UTC (rev 6835)
@@ -191,6 +191,11 @@
const double_array& upDir,
const int numLocs);
+ /** Get size (fiber dimension) of orientation information.
+ *
+ * @returns Size of orientation information.
+ */
+ int _orientationSize(void) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Added: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2007-05-10 01:34:20 UTC (rev 6835)
@@ -0,0 +1,274 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "FaultCohesiveDyn.hh" // implementation of object methods
+
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/utils/array.hh" // USES double_array
+
+#include <assert.h> // USES assert()
+#include <sstream> // USES std::ostringstream
+#include <stdexcept> // USES std::runtime_error
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::faults::FaultCohesiveDyn::FaultCohesiveDyn(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::faults::FaultCohesiveDyn::~FaultCohesiveDyn(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Copy constructor.
+pylith::faults::FaultCohesiveDyn::FaultCohesiveDyn(const FaultCohesiveDyn& f) :
+ FaultCohesive(f)
+{ // copy constructor
+} // copy constructor
+
+// ----------------------------------------------------------------------
+// Initialize fault. Determine orientation and setup boundary
+void
+pylith::faults::FaultCohesiveDyn::initialize(const ALE::Obj<ALE::Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ const double_array& upDir)
+{ // initialize
+ assert(0 != _quadrature);
+ assert(0 != _faultMesh);
+ assert(!_faultMesh->isNull());
+
+ if (3 != upDir.size())
+ throw std::runtime_error("Up direction for fault orientation must be "
+ "a vector with 3 components.");
+
+#if 0
+ // Allocate section for orientation at quadrature points
+ ALE::Obj<real_section_type> orientation =
+ new real_section_type((*_faultMesh)->comm(), (*_faultMesh)->debug());
+ assert(!orientation.isNull());
+ const int orientationSize = _orientationSize();
+ orientation->setFiberDimension((*_faultMesh)->depthStratum(0),
+ orientationSize);
+ (*_faultMesh)->allocate(orientation);
+
+ // Get section containing coordinates of vertices
+ const ALE::Obj<real_section_type>& coordinates =
+ mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
+ // Set orientation method
+ const int cellDim = _quadrature->cellDim();
+ const int spaceDim = _quadrature->spaceDim();
+ orient_fn_type orientFn;
+ switch (cellDim)
+ { // switch
+ case 1 :
+ orientFn = _orient1D;
+ break;
+ case 2 :
+ orientFn = _orient2D;
+ break;
+ case 3 :
+ orientFn = _orient3D;
+ break;
+ default :
+ assert(0);
+ } // switch
+
+ // Loop over cells, computing orientation at each vertex in cell
+ const ALE::Obj<sieve_type>& sieve = (*_faultMesh)->getSieve();
+ assert(!sieve.isNull());
+ const ALE::Obj<Mesh::label_sequence>& cells =
+ (*_faultMesh)->heightStratum(0);
+ const Mesh::label_sequence::iterator cBegin = cells->begin();
+ const Mesh::label_sequence::iterator cEnd = cells->end();
+ double_array cellOrientation(_quadrature->numBasis()*orientationSize);
+ const int numVertices = _quadrature->numBasis();
+ for (Mesh::label_sequence::iterator c_iter=cBegin;
+ c_iter != cEnd;
+ ++c_iter) {
+ // Compute cell geometry at vertices
+ _quadrature->computeGeometryVert(*_faultMesh, coordinates, *c_iter);
+
+ const double_array& jacobian = _quadrature->jacobianVert();
+ const double_array& jacobianDet = _quadrature->jacobianDetVert();
+
+ // Compute weighted orientation of face at vertices (using geometry info)
+ orientFn(&cellOrientation, jacobian, jacobianDet, upDir, numVertices);
+
+ // Update orientation section for vertices in cell
+ const ALE::Obj<sieve_type::traits::coneSequence>& cone =
+ sieve->cone(*c_iter);
+ assert(!cone.isNull());
+ const sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
+ const sieve_type::traits::coneSequence::iterator vEnd = cone->end();
+ int index = 0;
+ for(sieve_type::traits::coneSequence::iterator v_iter=vBegin;
+ v_iter != vEnd;
+ ++v_iter)
+ orientation->updatePoint(*v_iter, &cellOrientation[index]);
+ index += orientationSize;
+ } // for
+
+ // Assemble orientation information
+ //orientation->complete();
+
+ // Loop over vertices, make orientation information unit magnitude
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ (*_faultMesh)->depthStratum(0);
+ const Mesh::label_sequence::iterator vBegin = vertices->begin();
+ const Mesh::label_sequence::iterator vEnd = vertices->end();
+ double_array vertexDir(orientationSize);
+ for (Mesh::label_sequence::iterator v_iter=vBegin;
+ v_iter != vEnd;
+ ++v_iter) {
+ const real_section_type::value_type* vertexOrient =
+ orientation->restrictPoint(*v_iter);
+
+ assert(cellDim*spaceDim == orientationSize);
+ for (int iDim=0, index=0; iDim < cellDim; ++iDim, index+=cellDim) {
+ double mag = 0;
+ for (int jDim=0; jDim < spaceDim; ++jDim)
+ mag *= vertexOrient[index*cellDim+jDim];
+ for (int jDim=0; jDim < cellDim; ++jDim)
+ vertexDir[index*cellDim+jDim] = vertexOrient[index*cellDim+jDim] / mag;
+ } // for
+ orientation->updatePoint(*v_iter, &vertexDir[0]);
+ } // for
+
+ // Create set of constraint vertices
+ std::set<Mesh::point_type> setVert;
+ for (Mesh::label_sequence::iterator c_iter=cBegin;
+ c_iter != cEnd;
+ ++c_iter) {
+ // Vertices for each cohesive cell are in groups of N.
+ // 0 to N-1: vertices on negative side of the fault
+ // N-1 to 2N-1: vertices on positive side of the fault
+ // 2N to 3N-1: vertices associated with constraint forces
+ const ALE::Obj<sieve_type::traits::coneSequence>& cone =
+ sieve->cone(*c_iter);
+ assert(!cone.isNull());
+ const sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
+ const sieve_type::traits::coneSequence::iterator vEnd = cone->end();
+ const int coneSize = cone->size();
+ assert(coneSize % 3 == 0);
+ sieve_type::traits::coneSequence::iterator v_iter = vBegin;
+ // Skip over non-constraint vertices
+ for (int i=0, numSkip=2*coneSize/3; i < numSkip; ++i)
+ ++v_iter;
+ // Add constraint vertices to set
+ for(int i=0, numConstraintVert = coneSize/3;
+ i < numConstraintVert;
+ ++i, ++v_iter)
+ setVert.insert(*v_iter);
+ } // for
+
+ // Only store orientation information at constraint vertices
+ _orientation =
+ new real_section_type((*_faultMesh)->comm(), (*_faultMesh)->debug());
+ assert(!_orientation.isNull());
+ const std::set<Mesh::point_type>::const_iterator cvBegin =
+ _constraintVert.begin();
+ const std::set<Mesh::point_type>::const_iterator cvEnd =
+ _constraintVert.end();
+ for (std::set<Mesh::point_type>::const_iterator v_iter=cvBegin;
+ v_iter != cvEnd;
+ ++v_iter)
+ _orientation->setFiberDimension(*v_iter, orientationSize);
+ (*_faultMesh)->allocate(_orientation);
+ for (std::set<Mesh::point_type>::const_iterator v_iter=cvBegin;
+ v_iter != cvEnd;
+ ++v_iter) {
+ const real_section_type::value_type* vertexOrient =
+ orientation->restrictPoint(*v_iter);
+ _orientation->updatePoint(*v_iter, vertexOrient);
+ } // for
+#endif
+} // initialize
+
+// ----------------------------------------------------------------------
+// Integrate contribution of cohesive cells to residual term.
+void
+pylith::faults::FaultCohesiveDyn::integrateResidual(
+ const ALE::Obj<real_section_type>& residual,
+ const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateResidual
+#if 0
+ // Subtract constraint forces (which are disp at the constraint
+ // DOF) to residual; contributions are at DOF of normal vertices (i and j)
+
+ const ALE::Obj<Mesh::label_sequence>& cells =
+ (*_faultMesh)->heightStratum(0);
+ const Mesh::label_sequence::iterator cBegin = cells->begin();
+ const Mesh::label_sequence::iterator cEnd = cells->end();
+
+ // Allocate vector for cell values (if necessary)
+ _initCellVector();
+
+ // Loop over cohesive cells
+ const int numVertices = _quadrature->numBasis();
+ const int numConstraintVert = numVertices / 3;
+ assert(numVertices == numConstraintVert * 3);
+ for (Mesh::label_sequence::iterator c_iter=cBegin;
+ c_iter != cEnd;
+ ++c_iter) {
+ _resetCellVector();
+
+ // Get values at vertices (want constraint forces in disp vector)
+ const real_section_type::value_type* cellDisp =
+ mesh->restrict(disp, *c_iter);
+
+ // Transfer constraint forces to cell's constribution to residual vector
+ for (int i=0; i < numConstraintVert; ++i) {
+ const double constraintForce = cellDisp[2*numConstraintVert+i];
+ _cellVector[ i] = -constraintForce;
+ _cellVector[numConstraintVert+i] = -constraintForce;
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numConstraintVert*2);
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // Update residual
+ mesh->updateAdd(residual, *c_iter, _cellVector);
+ } // for
+#endif
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Compute Jacobian matrix (A) associated with operator.
+void
+pylith::faults::FaultCohesiveDyn::integrateJacobian(
+ PetscMat* mat,
+ const ALE::Obj<real_section_type>& dispT,
+ const ALE::Obj<Mesh>& mesh)
+{ // integrateJacobian
+} // integrateJacobian
+
+// ----------------------------------------------------------------------
+// Set field.
+void
+pylith::faults::FaultCohesiveDyn::setField(
+ const double t,
+ const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh)
+{ // setField
+} // setField
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2007-05-10 01:34:20 UTC (rev 6835)
@@ -0,0 +1,137 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/faults/FaultCohesiveDyn.hh
+ *
+ * @brief C++ implementation for a fault surface with spontaneous
+ * (dynamic) slip implemented with cohesive elements.
+ *
+ * The ordering of vertices in a cohesive cell is the vertices on the
+ * POSITIVE/NEGATIVE (CHECK WHICH IT IS) side of the fault and then the
+ * corresponding entries on the other side of the fault.
+ */
+
+#if !defined(pylith_faults_faultcohesivedyn_hh)
+#define pylith_faults_faultcohesivedyn_hh
+
+#include "FaultCohesive.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace faults {
+ class FaultCohesiveDyn;
+ class TestFaultCohesiveDyn; // unit testing
+ } // faults
+} // pylith
+
+/// @brief C++ implementation for a fault surface with spontaneous
+/// (dynamic) slip implemented with cohesive elements.
+class pylith::faults::FaultCohesiveDyn : public FaultCohesive
+{ // class FaultCohesiveDyn
+ friend class TestFaultCohesiveDyn; // unit testing
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor.
+ FaultCohesiveDyn(void);
+
+ /// Destructor.
+ virtual
+ ~FaultCohesiveDyn(void);
+
+ /** Create copy of fault.
+ *
+ * @returns Copy of fault.
+ */
+ Fault* clone(void) const;
+
+ /** Initialize fault. Determine orientation and setup boundary
+ * condition parameters.
+ *
+ * @param mesh PETSc mesh
+ * @param cs Coordinate system for mesh
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction).
+ */
+ void initialize(const ALE::Obj<ALE::Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs,
+ const double_array& upDir);
+
+ /** Integrate contribution of cohesive cells to residual term.
+ *
+ * @param residual Residual field (output)
+ * @param disp Displacement field at time t
+ * @param mesh Finite-element mesh
+ */
+ void integrateResidual(const ALE::Obj<real_section_type>& residual,
+ const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Compute Jacobian matrix (A) associated with operator.
+ *
+ * @param mat Sparse matrix
+ * @param disp Displacement field
+ * @param mesh Finite-element mesh
+ */
+ void integrateJacobian(PetscMat* mat,
+ const ALE::Obj<real_section_type>& dispT,
+ const ALE::Obj<Mesh>& mesh);
+
+ /** Set field.
+ *
+ * @param t Current time
+ * @param disp Displacement field
+ * @param mesh Finite-element mesh
+ */
+ void setField(const double t,
+ const ALE::Obj<real_section_type>& disp,
+ const ALE::Obj<Mesh>& mesh);
+
+ // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+ /** Copy constructor.
+ *
+ * @param m Fault to copy
+ */
+ FaultCohesiveDyn(const FaultCohesiveDyn& m);
+
+ /** Cohesive cells use Lagrange multiplier constraints?
+ *
+ * @returns True if implementation using Lagrange multiplier
+ * constraints, false otherwise.
+ */
+ bool _useLagrangeConstraints(void) const;
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ const FaultCohesiveDyn& operator=(const FaultCohesiveDyn& m);
+
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+ /// Orientation of fault surface at vertices (fiber dimension is
+ /// nonzero only at constraint vertices)
+ ALE::Obj<real_section_type> _orientation;
+
+}; // class FaultCohesiveDyn
+
+#include "FaultCohesiveDyn.icc" // inline methods
+
+#endif // pylith_faults_faultcohesivedyn_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc 2007-05-10 01:34:20 UTC (rev 6835)
@@ -0,0 +1,32 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_faults_faultcohesivedyn_hh)
+#error "FaultCohesiveDyn.icc can only be included from FaultCohesiveDyn.hh"
+#endif
+
+// Create copy of fault.
+inline
+pylith::faults::Fault*
+pylith::faults::FaultCohesiveDyn::clone(void) const {
+ return new FaultCohesiveDyn(*this);
+} // clone
+
+// Cohesive cells use Lagrange multiplier constraints?
+inline
+bool
+pylith::faults::FaultCohesiveDyn::_useLagrangeConstraints(void) const {
+ return false;
+} // useLagrangeConstraints
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-05-10 01:34:20 UTC (rev 6835)
@@ -367,16 +367,5 @@
disp->updatePoint(*v_iter, slip->restrictPoint(*v_iter));
} // setField
-// ----------------------------------------------------------------------
-// Get size (fiber dimension) of orientation information.
-int
-pylith::faults::FaultCohesiveKin::_orientationSize(void) const
-{ // _orientationSize
- assert(0 != _quadrature);
- const int size = _quadrature->cellDim()*_quadrature->spaceDim();
- return size;
-} // _orientationSize
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-05-10 01:34:20 UTC (rev 6835)
@@ -12,8 +12,8 @@
/** @file libsrc/faults/FaultCohesiveKin.hh
*
- * @brief C++ object implementing a fault surface with a kinematic
- * earthquake source.
+ * @brief C++ implementation for a fault surface with kinematic
+ * (prescribed) slip implemented with cohesive elements.
*
* Fault boundary condition is specified using Lagrange
* multipliers. The constraints are associated with "constraint"
@@ -134,15 +134,6 @@
/// Not implemented
const FaultCohesiveKin& operator=(const FaultCohesiveKin& m);
- // PRIVATE METHODS ////////////////////////////////////////////////////
-private :
-
- /** Get size (fiber dimension) of orientation information.
- *
- * @returns Size of orientation information.
- */
- int _orientationSize(void) const;
-
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/faults/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/Makefile.am 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/libsrc/faults/Makefile.am 2007-05-10 01:34:20 UTC (rev 6835)
@@ -21,6 +21,8 @@
Fault.hh \
Fault.icc \
FaultCohesive.hh \
+ FaultCohesiveDyn.hh \
+ FaultCohesiveDyn.icc \
FaultCohesiveKin.hh \
FaultCohesiveKin.icc \
SlipTimeFn.hh
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFault.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFault.cc 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFault.cc 2007-05-10 01:34:20 UTC (rev 6835)
@@ -17,7 +17,6 @@
#include "pylith/faults/FaultCohesiveKin.hh" // USES FaultCohesiveKin
#include <string> // USES std::string
-#include <stdexcept> // USES std::logic_error TEMPORARY
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestFault );
@@ -46,15 +45,5 @@
CPPUNIT_ASSERT_EQUAL(label, fault.label());
} // testLabel
-#if 0
-// ----------------------------------------------------------------------
-// Test initialize()
-void
-pylith::faults::TestFault::testInitialize(void)
-{ // testInitialize
- throw std::logic_error("Need to implement test.");
-} // testInitialize
-#endif
-
// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFault.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFault.hh 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFault.hh 2007-05-10 01:34:20 UTC (rev 6835)
@@ -39,9 +39,6 @@
CPPUNIT_TEST_SUITE( TestFault );
CPPUNIT_TEST( testID );
CPPUNIT_TEST( testLabel );
-#if 0
- CPPUNIT_TEST( testInitialize );
-#endif
CPPUNIT_TEST_SUITE_END();
// PUBLIC METHODS /////////////////////////////////////////////////////
@@ -53,11 +50,6 @@
/// Test label()
void testLabel(void);
-#if 0
- /// Test initialize()
- void testInitialize(void);
-#endif
-
}; // class TestFault
#endif // pylith_faults_testfault_hh
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc 2007-05-10 01:34:20 UTC (rev 6835)
@@ -15,10 +15,12 @@
#include "TestFaultCohesive.hh" // Implementation of class methods
#include "pylith/faults/FaultCohesiveKin.hh" // USES FaultsCohesiveKin
+#include "pylith/faults/FaultCohesiveDyn.hh" // USES FaultsCohesiveDyn
#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
-#include "pylith/utils/array.hh" // USES int_array
+#include "pylith/utils/array.hh" // USES int_array, double_array
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
+#include "pylith/feassemble/Quadrature2Din3D.hh" // USES Quadrature2Din3D
#include "data/CohesiveDataLine2.hh" // USES CohesiveDataLine2
#include "data/CohesiveDataTri3.hh" // USES CohesiveDataTri3
@@ -32,6 +34,8 @@
#include "data/CohesiveDataTet4Lagrange.hh" // USES CohesiveDataTet4Lagrange
#include "data/CohesiveDataHex8Lagrange.hh" // USES CohesiveDataHex8Lagrange
+#include <stdexcept> // TEMPORARY
+
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestFaultCohesive );
@@ -41,7 +45,8 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyLine2(void)
{ // testAdjustTopologyLine2
CohesiveDataLine2 data;
- _testAdjustTopology(data);
+ FaultCohesiveDyn fault;
+ _testAdjustTopology(&fault, data);
} // testAdjustTopologyLine2
// ----------------------------------------------------------------------
@@ -50,7 +55,8 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTri3(void)
{ // testAdjustTopologyTri3
CohesiveDataTri3 data;
- _testAdjustTopology(data);
+ FaultCohesiveDyn fault;
+ _testAdjustTopology(&fault, data);
} // testAdjustTopologyTri3
// ----------------------------------------------------------------------
@@ -59,7 +65,8 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyQuad4(void)
{ // testAdjustTopologyQuad4
CohesiveDataQuad4 data;
- _testAdjustTopology(data);
+ FaultCohesiveDyn fault;
+ _testAdjustTopology(&fault, data);
} // testAdjustTopologyQuad4
// ----------------------------------------------------------------------
@@ -68,7 +75,8 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTet4(void)
{ // testAdjustTopologyTet4
CohesiveDataTet4 data;
- _testAdjustTopology(data);
+ FaultCohesiveDyn fault;
+ _testAdjustTopology(&fault, data);
} // testAdjustTopologyTet4
// ----------------------------------------------------------------------
@@ -77,7 +85,8 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyHex8(void)
{ // testAdjustTopologyHex8
CohesiveDataHex8 data;
- _testAdjustTopology(data);
+ FaultCohesiveDyn fault;
+ _testAdjustTopology(&fault, data);
} // testAdjustTopologyHex8
// ----------------------------------------------------------------------
@@ -87,7 +96,8 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyLine2Lagrange(void)
{ // testAdjustTopologyLine2Lagrange
CohesiveDataLine2Lagrange data;
- _testAdjustTopology(data);
+ FaultCohesiveKin fault;
+ _testAdjustTopology(&fault, data);
} // testAdjustTopologyLine2Lagrange
// ----------------------------------------------------------------------
@@ -97,7 +107,8 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTri3Lagrange(void)
{ // testAdjustTopologyTri3Lagrange
CohesiveDataTri3Lagrange data;
- _testAdjustTopology(data);
+ FaultCohesiveKin fault;
+ _testAdjustTopology(&fault, data);
} // testAdjustTopologyTri3Lagrange
// ----------------------------------------------------------------------
@@ -107,7 +118,8 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyQuad4Lagrange(void)
{ // testAdjustTopologyQuad4Lagrange
CohesiveDataQuad4Lagrange data;
- _testAdjustTopology(data);
+ FaultCohesiveKin fault;
+ _testAdjustTopology(&fault, data);
} // testAdjustTopologyQuad4Lagrange
// ----------------------------------------------------------------------
@@ -117,7 +129,8 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTet4Lagrange(void)
{ // testAdjustTopologyTet4Lagrange
CohesiveDataTet4Lagrange data;
- _testAdjustTopology(data);
+ FaultCohesiveKin fault;
+ _testAdjustTopology(&fault, data);
} // testAdjustTopologyTet4Lagrange
// ----------------------------------------------------------------------
@@ -127,13 +140,137 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyHex8Lagrange(void)
{ // testAdjustTopologyHex8Lagrange
CohesiveDataHex8Lagrange data;
- _testAdjustTopology(data);
+ FaultCohesiveKin fault;
+ _testAdjustTopology(&fault, data);
} // testAdjustTopologyHex8Lagrange
// ----------------------------------------------------------------------
+// Test _orientationSize().
+void
+pylith::faults::TestFaultCohesive::testOrientationSize(void)
+{ // testOrientationSize
+ const int cellDim = 2;
+ const int spaceDim = 3;
+ const int numBasis = 1;
+ const int numQuadPts = 1;
+ const double basisVert[] = { 0.5, 0.3, 0.7 };
+ const double basisDerivVert[] = { -0.5, 0.5, -0.4 };
+ const double basisQuad[] = { 0.5, 0.5, 0.4 };
+ const double basisDerivQuad[] = { 0.5, 0.3, -0.4 };
+ const double quadPtsRef[] = { 0.0, 3.0 };
+ const double quadWts[] = { 2.0 };
+ const double minJacobian = 1.0;
+
+ feassemble::Quadrature2Din3D q;
+ q.initialize(basisVert, basisDerivVert,
+ basisQuad, basisDerivQuad, quadPtsRef, quadWts,
+ cellDim, numBasis, numQuadPts, spaceDim);
+
+ FaultCohesiveKin fault;
+ fault.quadrature(&q);
+ CPPUNIT_ASSERT_EQUAL(cellDim*spaceDim, fault._orientationSize());
+} // testOrientationSize
+
+// ----------------------------------------------------------------------
+// Test _orient1D().
+void
+pylith::faults::TestFaultCohesive::testOrient1D(void)
+{ // testOrient1D
+ const int numLocs = 3;
+ double_array jacobian;
+ double_array jacobianDet;
+ double_array upDir;
+ double_array orientation(numLocs);
+
+ FaultCohesive::_orient1D(&orientation,
+ jacobian, jacobianDet, upDir, numLocs);
+
+ const int size = orientation.size();
+ CPPUNIT_ASSERT_EQUAL(numLocs, size);
+ const double tolerance = 1.0e-6;
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, orientation[i], tolerance);
+} // testOrient1D
+
+// ----------------------------------------------------------------------
+// Test _orient2D().
+void
+pylith::faults::TestFaultCohesive::testOrient2D(void)
+{ // testOrient2D
+ const int numLocs = 2;
+ const int spaceDim = 2;
+ const int orientSize = 4;
+
+ const double jacobianVals[] = {
+ -1.0, 2.0,
+ -0.5, 1.0
+ };
+ double_array jacobian(jacobianVals, numLocs*spaceDim*(spaceDim-1));
+ double_array jacobianDet;
+ double_array upDir;
+ double_array orientation(numLocs*orientSize);
+
+ const double orientationE[] = {
+ -1.0, 2.0, 2.0, 1.0,
+ -0.5, 1.0, 1.0, 0.5
+ };
+
+ FaultCohesive::_orient2D(&orientation,
+ jacobian, jacobianDet, upDir, numLocs);
+
+ const int size = orientation.size();
+ CPPUNIT_ASSERT_EQUAL(numLocs*orientSize, size);
+ const double tolerance = 1.0e-6;
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(orientationE[i], orientation[i], tolerance);
+} // testOrient2D
+
+// ----------------------------------------------------------------------
+// Test _orient3D().
+void
+pylith::faults::TestFaultCohesive::testOrient3D(void)
+{ // testOrient3D
+ const int numLocs = 2;
+ const int spaceDim = 3;
+ const int orientSize = 9;
+
+ const double jacobianVals[] = {
+ 2.0, 1.0, 0.5, -0.5, -0.2, 2.0,
+ };
+ double_array jacobian(jacobianVals, numLocs*spaceDim*(spaceDim-1));
+ const double jacobianDetVals[] = {
+ 1.3, 0.7
+ };
+ double_array jacobianDet(jacobianDetVals, numLocs);
+ const double upDirVals[] = { 0.0, 0.0, 1.0 };
+ double_array upDir(upDirVals, 3);
+ double_array orientation(numLocs*orientSize);
+
+ const double orientationE[] = {
+ 1.1654847299258313, -0.012145479112634533, 0.57575848378190342,
+ 0.57588657243394026, 0.024580136299379406, -1.1652255028919474,
+ 0.0, 1.2997108540889502, 0.027417070656281111,
+
+ -0.063065020894967128, 0.68889496382717197, -0.10698846644884991,
+ -0.6957991174916025, -0.0688894963827172, -0.033433895765265592,
+ -0.043432605694620839, 0.10333424457407581, 0.69096717914882233
+ };
+
+ FaultCohesive::_orient3D(&orientation,
+ jacobian, jacobianDet, upDir, numLocs);
+
+ const int size = orientation.size();
+ CPPUNIT_ASSERT_EQUAL(numLocs*orientSize, size);
+ const double tolerance = 1.0e-6;
+ for (int i=0; i < size; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(orientationE[i], orientation[i], tolerance);
+} // testOrient3D
+
+// ----------------------------------------------------------------------
// Test adjustTopology().
void
-pylith::faults::TestFaultCohesive::_testAdjustTopology(const CohesiveData& data)
+pylith::faults::TestFaultCohesive::_testAdjustTopology(Fault* fault,
+ const CohesiveData& data)
{ // _testAdjustTopology
ALE::Obj<ALE::Mesh> mesh;
meshio::MeshIOAscii iohandler;
@@ -142,10 +279,10 @@
iohandler.interpolate(false);
iohandler.read(&mesh);
- FaultCohesiveKin fault;
- fault.id(1);
- fault.label("fault");
- fault.adjustTopology(mesh);
+ CPPUNIT_ASSERT(0 != fault);
+ fault->id(1);
+ fault->label("fault");
+ fault->adjustTopology(mesh);
CPPUNIT_ASSERT_EQUAL(data.cellDim, mesh->getDimension());
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.hh 2007-05-10 00:42:52 UTC (rev 6834)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.hh 2007-05-10 01:34:20 UTC (rev 6835)
@@ -41,14 +41,11 @@
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
CPPUNIT_TEST_SUITE( TestFaultCohesive );
-#if 0
- // awaiting FaultCohesiveDyn implementation
CPPUNIT_TEST( testAdjustTopologyLine2 );
CPPUNIT_TEST( testAdjustTopologyTri3 );
CPPUNIT_TEST( testAdjustTopologyQuad4 );
CPPUNIT_TEST( testAdjustTopologyTet4 );
CPPUNIT_TEST( testAdjustTopologyHex8 );
-#endif
CPPUNIT_TEST( testAdjustTopologyLine2Lagrange );
CPPUNIT_TEST( testAdjustTopologyTri3Lagrange );
@@ -56,6 +53,12 @@
CPPUNIT_TEST( testAdjustTopologyTet4Lagrange );
CPPUNIT_TEST( testAdjustTopologyHex8Lagrange );
+ CPPUNIT_TEST( testOrientationSize );
+
+ CPPUNIT_TEST( testOrient1D );
+ CPPUNIT_TEST( testOrient2D );
+ CPPUNIT_TEST( testOrient3D );
+
CPPUNIT_TEST_SUITE_END();
// PUBLIC METHODS /////////////////////////////////////////////////////
@@ -96,11 +99,28 @@
/// multipliers.
void testAdjustTopologyHex8Lagrange(void);
+ /// Test _orientationSize().
+ void testOrientationSize(void);
+
+ /// Test _orient1D().
+ void testOrient1D(void);
+
+ /// Test _orient2D().
+ void testOrient2D(void);
+
+ /// Test _orient3D().
+ void testOrient3D(void);
+
// PROTECTED METHODS //////////////////////////////////////////////////
public :
- /// Test adjustTopology().
- void _testAdjustTopology(const CohesiveData& data);
+ /** Test adjustTopology().
+ *
+ * @param fault Fault for cohesive elements.
+ * @param data Cohesive element data.
+ */
+ void _testAdjustTopology(Fault* fault,
+ const CohesiveData& data);
}; // class TestFaultCohesive
More information about the cig-commits
mailing list