[cig-commits] r6701 - in short/3D/PyLith/trunk: . libsrc/faults
brad at geodynamics.org
brad at geodynamics.org
Thu Apr 26 16:47:56 PDT 2007
Author: brad
Date: 2007-04-26 16:47:56 -0700 (Thu, 26 Apr 2007)
New Revision: 6701
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc
short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh
short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc
short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh
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/SlipTimeFn.hh
Log:
Worked on implementing cohesive cells for kinematic earthquake source.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-04-26 22:06:24 UTC (rev 6700)
+++ short/3D/PyLith/trunk/TODO 2007-04-26 23:47:56 UTC (rev 6701)
@@ -19,31 +19,43 @@
Use MeshIOAscii to get mesh information into Python tests?
-0. Start implementing faults
+0. Implement faults for kinematic source
a. Add creation of cohesive cells
i. Add tests for interpolated meshes and constrained cohesive cells.
+ Double check consistency in ordering of vertices (positive/negative).
b. Implement slip time function
- i. Add Python
- ii. Add unit tests
+ i. Allow 1-D, 2-D, or 3-D specification of slip.
+ 1-D uses "slip"
+ 2-D uses "slip", "fault-opening"
+ 3-D uses "left-lateral", "reverse", "fault-opening"
+ ii. Add Python
+ iii. Add unit tests
+
c. Start implementing integrator for faults
+
i. Update Reference cell
(1) Compute basis and basisDeriv at vertices
(2) Add unit tests
+
ii. Update Quadrature
(1) Add basis and basisDeriv at vertices to Quadrature
(2) Add computeGeometryQuad/computeGeometryVert to Quadrature
(3) Add unit tests
+
iii. FaultCohesive
- (1) orientXD
- (2) Unit tests
+ (1) Unit tests
+
iv. FaultCohesiveKin
- (1) initialize(), setField(), integrateResidual(), integrateJacobian()
+ (1) initialize(), integrateResidual(), integrateJacobian()
(2) Add Python unit tests
(3) Add C++ unit tests
1. Finish implementing ExplicitElasticity and Explicit
a. Replace integrateConstant() with integrateResidual()
+
+ {f}-[A]{u} {u} is "guess" (zero for explicit)
+
a. Double check loops for integrateResidual() and integrateJacobian()
b. Create unit test (construction of residual term and Jacobian)
Modified: short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc 2007-04-26 22:06:24 UTC (rev 6700)
+++ short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc 2007-04-26 23:47:56 UTC (rev 6701)
@@ -57,8 +57,12 @@
void
pylith::faults::BruneSlipFn::initialize(const ALE::Obj<Mesh>& mesh,
const ALE::Obj<Mesh>& faultMesh,
+ const std::vector<Mesh::point_type>& vertices,
+
const spatialdata::geocoords::CoordSys* cs)
{ // initialize
+ typedef std::vector<Mesh::point_type>::const_iterator vert_iterator;
+
assert(!mesh.isNull());
assert(!faultMesh.isNull());
assert(0 != cs);
@@ -66,9 +70,6 @@
assert(0 != _dbSlipTime);
assert(0 != _dbPeakRate);
- // Get fault vertices
- const ALE::Obj<Mesh::label_sequence>& vertices = faultMesh->depthStratum(0);
-
// Create and allocate sections for parameters
delete _parameters;
_parameters = new feassemble::ParameterManager(faultMesh);
@@ -124,11 +125,9 @@
double slipData[3];
double slipTimeData;
double peakRateData;
- const Mesh::label_sequence::iterator vBegin = vertices->begin();
- const Mesh::label_sequence::iterator vEnd = vertices->end();
- for (Mesh::label_sequence::iterator v_iter=vBegin;
- v_iter != vEnd;
- ++v_iter) {
+ const vert_iterator vBegin = vertices.begin();
+ const vert_iterator vEnd = vertices.end();
+ for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter) {
// Get coordinates of vertex
const real_section_type::value_type* vCoords =
coordinates->restrictPoint(*v_iter);
@@ -182,14 +181,13 @@
// Get slip on fault surface at time t.
const ALE::Obj<pylith::real_section_type>&
pylith::faults::BruneSlipFn::slip(const double t,
- const ALE::Obj<Mesh>& faultMesh)
+ const std::vector<Mesh::point_type>& vertices)
{ // slip
+ typedef std::vector<Mesh::point_type>::const_iterator vert_iterator;
+
assert(0 != _parameters);
assert(!_slipField.isNull());
- // Get fault vertices
- const ALE::Obj<Mesh::label_sequence>& vertices = faultMesh->depthStratum(0);
-
// Get parameters
const ALE::Obj<real_section_type>& finalSlip =
_parameters->getReal("final slip");
@@ -204,11 +202,9 @@
assert(!peakRate.isNull());
double slipValues[3];
- const Mesh::label_sequence::iterator vBegin = vertices->begin();
- const Mesh::label_sequence::iterator vEnd = vertices->end();
- for (Mesh::label_sequence::iterator v_iter=vBegin;
- v_iter != vEnd;
- ++v_iter) {
+ const vert_iterator vBegin = vertices.begin();
+ const vert_iterator vEnd = vertices.end();
+ for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter) {
// Get values of parameters at vertex
const real_section_type::value_type* vFinalSlip =
finalSlip->restrictPoint(*v_iter);
Modified: short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh 2007-04-26 22:06:24 UTC (rev 6700)
+++ short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh 2007-04-26 23:47:56 UTC (rev 6701)
@@ -83,20 +83,22 @@
/** Initialize slip time function.
*
* @param mesh Finite-element mesh.
- * @param faultMesh Fault finite-element mesh.
+ * @param faultMesh Finite-element mesh of fault.
+ * @param vertices Vertices where slip will be prescribed.
* @param cs Coordinate system for mesh
*/
void initialize(const ALE::Obj<Mesh>& mesh,
const ALE::Obj<Mesh>& faultMesh,
+ const std::vector<Mesh::point_type>& vertices,
const spatialdata::geocoords::CoordSys* cs);
/** Get slip on fault surface at time t.
*
* @param t Time t.
- * @param mesh Fault finite-element mesh.
+ * @param vertices Vertices where slip will be prescribed.
*/
const ALE::Obj<real_section_type>& slip(const double t,
- const ALE::Obj<Mesh>& faultMesh);
+ const std::vector<Mesh::point_type>& vertices);
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc 2007-04-26 22:06:24 UTC (rev 6700)
+++ short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc 2007-04-26 23:47:56 UTC (rev 6701)
@@ -52,23 +52,24 @@
// ----------------------------------------------------------------------
// Initialize slip time function.
void
-pylith::faults::EqKinSrc::initialize(const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<Mesh>& faultMesh,
-
- const spatialdata::geocoords::CoordSys* cs)
+pylith::faults::EqKinSrc::initialize(
+ const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<Mesh>& faultMesh,
+ const std::vector<Mesh::point_type>& vertices,
+ const spatialdata::geocoords::CoordSys* cs)
{ // initialize
assert(0 != _slipfn);
- _slipfn->initialize(mesh, faultMesh, cs);
+ _slipfn->initialize(mesh, faultMesh, vertices, cs);
} // initialize
// ----------------------------------------------------------------------
// Get slip on fault surface at time t.
const ALE::Obj<pylith::real_section_type>&
pylith::faults::EqKinSrc::slip(const double t,
- const ALE::Obj<Mesh>& faultMesh)
+ const std::vector<Mesh::point_type>& vertices)
{ // slip
assert(0 != _slipfn);
- return _slipfn->slip(t, faultMesh);
+ return _slipfn->slip(t, vertices);
} // slip
Modified: short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh 2007-04-26 22:06:24 UTC (rev 6700)
+++ short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh 2007-04-26 23:47:56 UTC (rev 6701)
@@ -76,22 +76,24 @@
/** Initialize slip time function.
*
* @param mesh Finite-element mesh.
- * @param faultMesh Fault finite-element mesh.
+ * @param faultMesh Finite-element mesh of fault.
+ * @param vertices Vertices where slip will be prescribed.
* @param cs Coordinate system for mesh
*/
virtual
void initialize(const ALE::Obj<Mesh>& mesh,
const ALE::Obj<Mesh>& faultMesh,
+ const std::vector<Mesh::point_type>& vertices,
const spatialdata::geocoords::CoordSys* cs);
/** Get slip on fault surface at time t.
*
* @param t Time t.
- * @param mesh Fault finite-element mesh.
+ * @param vertices Vertices where slip will be prescribed.
*/
virtual
const ALE::Obj<real_section_type>& slip(const double t,
- const ALE::Obj<Mesh>& faultMesh);
+ const std::vector<Mesh::point_type>& vertices);
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2007-04-26 22:06:24 UTC (rev 6700)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2007-04-26 23:47:56 UTC (rev 6701)
@@ -17,6 +17,7 @@
#include "CohesiveTopology.hh" // USES CohesiveTopology::create()
#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+#include "pylith/utils/array.hh" // USES double_array
#include <assert.h> // USES assert()
@@ -59,5 +60,113 @@
CohesiveTopology::create(_faultMesh, mesh, groupField, id());
} // adjustTopology
+// ----------------------------------------------------------------------
+// Compute weighted orientation of fault for cohesive cell between
+// 1-D elements.
+void
+pylith::faults::FaultCohesive::_orient1D(double_array* orientation,
+ const double_array& jacobian,
+ const double_array& jacobianDet,
+ const double_array& upDir,
+ const int numLocs)
+{ // _orient1D
+ assert(0 != orientation);
+ assert(numLocs == orientation->size());
+ (*orientation) = 1.0;
+} // _orient1D
+
+// ----------------------------------------------------------------------
+// Compute weighted orientation of fault for cohesive cell between
+// 2-D elements.
+void
+pylith::faults::FaultCohesive::_orient2D(double_array* orientation,
+ const double_array& jacobian,
+ const double_array& jacobianDet,
+ const double_array& upDir,
+ const int numLocs)
+{ // _orient2D
+ const int orientSize = 4;
+ assert(0 != orientation);
+ assert(orientSize*numLocs == orientation->size());
+ const int jacobianSize = 2;
+ assert(numLocs*jacobianSize == jacobian.size());
+ // cellDim is 1
+ const int spaceDim = 2;
+ for (int iLoc=0; iLoc < numLocs; ++iLoc) {
+ const double j1 = jacobian[iLoc*spaceDim ];
+ const double j2 = jacobian[iLoc*spaceDim+1];
+ (*orientation)[iLoc*orientSize ] = j1;
+ (*orientation)[iLoc*orientSize+1] = j2;
+ (*orientation)[iLoc*orientSize+2] = j2;
+ (*orientation)[iLoc*orientSize+3] = -j1;
+ } // for
+} // _orient2D
+
+// ----------------------------------------------------------------------
+// Compute weighted orientation of fault for cohesive cell between
+// 3-D elements.
+void
+pylith::faults::FaultCohesive::_orient3D(double_array* orientation,
+ const double_array& jacobian,
+ const double_array& jacobianDet,
+ const double_array& upDir,
+ const int numLocs)
+{ // _orient3D
+ const int orientSize = 9;
+ assert(0 != orientation);
+ assert(orientSize*numLocs == orientation->size());
+ const int jacobianSize = 6;
+ assert(numLocs*jacobianSize == jacobian.size());
+ assert(numLocs == jacobianDet.size());
+
+ const int cellDim = 2;
+ const int spaceDim = 3;
+ for (int iLoc=0; iLoc < numLocs; ++iLoc) {
+ const double j00 = jacobian[iLoc*spaceDim ];
+ const double j10 = jacobian[iLoc*spaceDim+1];
+ const double j20 = jacobian[iLoc*spaceDim+2];
+ const double j01 = jacobian[iLoc*spaceDim+3];
+ const double j11 = jacobian[iLoc*spaceDim+4];
+ const double j21 = jacobian[iLoc*spaceDim+5];
+
+ // Compute normal using Jacobian
+ double r0 = j10*j21 - j20*j11;
+ double r1 = -j00*j21 + j20*j01;
+ double r2 = j00*j11 - j10*j01;
+ // Make unit vector
+ double mag = sqrt(r0*r0 + r1*r1 + r2*r2);
+ r0 /= mag;
+ r1 /= mag;
+ r2 /= mag;
+
+ // Compute along-strike direction; cross product of "up" and normal
+ double p0 = upDir[1]*r2 - upDir[2]*r1;
+ double p1 = -upDir[0]*r2 + upDir[2]*r0;
+ double p2 = upDir[0]*r1 - upDir[1]*r0;
+ // Make unit vector
+ mag = sqrt(p0*p0 + p1*p1 + p2*p2);
+ p0 /= mag;
+ p1 /= mag;
+ p2 /= mag;
+
+ // Compute up-dip direction; cross product of normal and along-strike
+ const double q0 = r1*p2 - r2*p1;
+ const double q1 = -r0*p2 + r2*p0;
+ const double q2 = r0*p1 - r1*p0;
+
+ const double wt = jacobianDet[iLoc];
+ (*orientation)[iLoc*orientSize ] = p0*wt;
+ (*orientation)[iLoc*orientSize+1] = q0*wt;
+ (*orientation)[iLoc*orientSize+2] = r0*wt;
+ (*orientation)[iLoc*orientSize+3] = p1*wt;
+ (*orientation)[iLoc*orientSize+4] = q1*wt;
+ (*orientation)[iLoc*orientSize+5] = r1*wt;
+ (*orientation)[iLoc*orientSize+6] = p2*wt;
+ (*orientation)[iLoc*orientSize+7] = q2*wt;
+ (*orientation)[iLoc*orientSize+8] = r2*wt;
+ } // for
+} // _orient3D
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2007-04-26 22:06:24 UTC (rev 6700)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2007-04-26 23:47:56 UTC (rev 6701)
@@ -92,11 +92,13 @@
/** Set field.
*
+ * @param t Current time
* @param disp Displacement field
* @param mesh Finite-element mesh
*/
virtual
- void setField(const ALE::Obj<real_section_type>& disp,
+ void setField(const double t,
+ const ALE::Obj<real_section_type>& disp,
const ALE::Obj<Mesh>& mesh) = 0;
// PROTECTED METHODS //////////////////////////////////////////////////
@@ -124,11 +126,13 @@
* @param upDir Direction perpendicular to along-strike direction that is
* not collinear with fault normal (usually "up" direction but could
* be up-dip direction).
+ * @param numLocs Number of locations where values are given.
*/
- void orient1D(double_array* orientation,
- const double_array& jacobian,
- const double_array& jacobianDet,
- const double_array& upDir);
+ void _orient1D(double_array* orientation,
+ const double_array& jacobian,
+ const double_array& jacobianDet,
+ const double_array& upDir,
+ const int numLocs);
/** Compute weighted orientation of fault for cohesive cell between
* 2-D elements. Orientation is either at vertices or quadrature
@@ -146,11 +150,13 @@
* @param upDir Direction perpendicular to along-strike direction that is
* not collinear with fault normal (usually "up" direction but could
* be up-dip direction).
+ * @param numLocs Number of locations where values are given.
*/
- void orient2D(double_array* orientation,
- const double_array& jacobian,
- const double_array& jacobianDet,
- const double_array& upDir);
+ void _orient2D(double_array* orientation,
+ const double_array& jacobian,
+ const double_array& jacobianDet,
+ const double_array& upDir,
+ const int numLocs);
/** Compute weighted orientation of fault for cohesive cell between
* 3-D elements. Orientation is either at vertices or quadrature
@@ -168,11 +174,13 @@
* @param upDir Direction perpendicular to along-strike direction that is
* not collinear with fault normal (usually "up" direction but could
* be up-dip direction).
+ * @param numLocs Number of locations where values are given.
*/
- void orient3D(double_array* orientation,
- const double_array& jacobian,
- const double_array& jacobianDet,
- const double_array& upDir);
+ void _orient3D(double_array* orientation,
+ const double_array& jacobian,
+ const double_array& jacobianDet,
+ const double_array& upDir,
+ const int numLocs);
// NOT IMPLEMENTED ////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-04-26 22:06:24 UTC (rev 6700)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-04-26 23:47:56 UTC (rev 6701)
@@ -66,8 +66,21 @@
throw std::runtime_error("Up direction for fault orientation must be "
"a vector with 3 components.");
-
-
+ // Allocate section for orientation information
+
+ // Loop over cells
+ // Compute cell geometry at vertices
+ // Compute weighted orientation of face at vertices (using geometry info)
+ // Update weighted orientations (wt is |J|)
+
+ // Assemble orientation information
+
+ // Loop over vertices
+ // Make orientation information unit magnitude
+
+ // Create list of constraint vertices
+
+ // Only store orientation information at constraint vertices
} // initialize
// ----------------------------------------------------------------------
@@ -78,6 +91,10 @@
const ALE::Obj<real_section_type>& disp,
const ALE::Obj<Mesh>& mesh)
{ // integrateResidual
+
+ // Subtract constraint forces (which are in disp at the constraint
+ // DOF) to residual; contributions are at DOF of normal vertices (i and j)
+
} // integrateResidual
// ----------------------------------------------------------------------
@@ -88,15 +105,31 @@
const ALE::Obj<real_section_type>& dispT,
const ALE::Obj<Mesh>& mesh)
{ // integrateJacobian
+
+ // Add constraint information to Jacobian matrix; these are the
+ // direction cosines. Entries are associated with vertices ik, jk,
+ // ki, and kj.
+
} // integrateJacobian
// ----------------------------------------------------------------------
// Set field.
void
pylith::faults::FaultCohesiveKin::setField(
+ const double t,
const ALE::Obj<real_section_type>& disp,
const ALE::Obj<Mesh>& mesh)
{ // setField
+ typedef std::vector<Mesh::point_type>::const_iterator vert_iterator;
+
+ assert(0 != _eqsrc);
+
+ const ALE::Obj<real_section_type>& slip = _eqsrc->slip(t, _constraintVert);
+ assert(!slip.isNull());
+ const vert_iterator vBegin = _constraintVert.begin();
+ const vert_iterator vEnd = _constraintVert.end();
+ for (vert_iterator v_iter=vBegin; v_iter != vEnd; ++v_iter)
+ disp->updatePoint(*v_iter, slip->restrictPoint(*v_iter));
} // setField
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-04-26 22:06:24 UTC (rev 6700)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-04-26 23:47:56 UTC (rev 6701)
@@ -12,8 +12,18 @@
/** @file libsrc/faults/FaultCohesiveKin.hh
*
- * @brief C++ abstract base class for a fault surface implemented with
- * cohesive elements.
+ * @brief C++ object implementing a fault surface with a kinematic
+ * earthquake source.
+ *
+ * Fault boundary condition is specified using Lagrange
+ * multipliers. The constraints are associated with "constraint"
+ * vertices which sit between the pair of vertices on each side of the
+ * fault.
+ *
+ * The ordering of vertices in a cohesive cell is the vertices on the
+ * POSITIVE/NEGATIVE (CHECK WHICH IT IS) side of the fault, the
+ * corresponding entries on the other side of the fault, and then the
+ * corresponding constraint vertices.
*/
#if !defined(pylith_faults_faultcohesivekin_hh)
@@ -92,10 +102,12 @@
/** Set field.
*
+ * @param t Current time
* @param disp Displacement field
* @param mesh Finite-element mesh
*/
- void setField(const ALE::Obj<real_section_type>& disp,
+ void setField(const double t,
+ const ALE::Obj<real_section_type>& disp,
const ALE::Obj<Mesh>& mesh);
// PROTECTED METHODS //////////////////////////////////////////////////
@@ -118,6 +130,9 @@
EqKinSrc* _eqsrc; ///< Kinematic earthquake source information
+ /// Fault vertices associated with constraints
+ std::vector<Mesh::point_type> _constraintVert;
+
}; // class FaultCohesiveKin
#include "FaultCohesiveKin.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh 2007-04-26 22:06:24 UTC (rev 6700)
+++ short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh 2007-04-26 23:47:56 UTC (rev 6701)
@@ -66,22 +66,24 @@
/** Initialize slip time function.
*
* @param mesh Finite-element mesh.
- * @param faultMesh Fault finite-element mesh.
+ * @param faultMesh Finite-element mesh of fault.
+ * @param vertices Vertices where slip will be prescribed.
* @param cs Coordinate system for mesh
*/
virtual
void initialize(const ALE::Obj<Mesh>& mesh,
const ALE::Obj<Mesh>& faultMesh,
+ const std::vector<Mesh::point_type>& vertices,
const spatialdata::geocoords::CoordSys* cs) = 0;
/** Get slip on fault surface at time t.
*
* @param t Time t.
- * @param mesh Fault finite-element mesh.
+ * @param vertices Vertices where slip will be prescribed.
*/
virtual
const ALE::Obj<real_section_type>& slip(const double t,
- const ALE::Obj<Mesh>& faultMesh) = 0;
+ const std::vector<Mesh::point_type>& vertices) = 0;
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
More information about the cig-commits
mailing list