[cig-commits] r6742 - in short/3D/PyLith/trunk: . libsrc/faults
libsrc/feassemble
brad at geodynamics.org
brad at geodynamics.org
Tue May 1 14:34:27 PDT 2007
Author: brad
Date: 2007-05-01 14:34:23 -0700 (Tue, 01 May 2007)
New Revision: 6742
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.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
short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
Log:
Initial C++ implementation of kinematic earthquake source. Not tested.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-05-01 20:27:04 UTC (rev 6741)
+++ short/3D/PyLith/trunk/TODO 2007-05-01 21:34:23 UTC (rev 6742)
@@ -31,7 +31,7 @@
c. Start implementing integrator for faults
i. Update Reference cell
- (1) Get vertices from FIAT [requested Matt fix]
+ (1) Improve how we get Jacobian at vertices of cell
iii. FaultCohesive
(1) Unit tests
Modified: short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc 2007-05-01 20:27:04 UTC (rev 6741)
+++ short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc 2007-05-01 21:34:23 UTC (rev 6742)
@@ -58,11 +58,11 @@
void
pylith::faults::BruneSlipFn::initialize(const ALE::Obj<Mesh>& mesh,
const ALE::Obj<Mesh>& faultMesh,
- const std::vector<Mesh::point_type>& vertices,
+ const std::set<Mesh::point_type>& vertices,
const spatialdata::geocoords::CoordSys* cs)
{ // initialize
- typedef std::vector<Mesh::point_type>::const_iterator vert_iterator;
+ typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
assert(!mesh.isNull());
assert(!faultMesh.isNull());
@@ -198,9 +198,9 @@
// Get slip on fault surface at time t.
const ALE::Obj<pylith::real_section_type>&
pylith::faults::BruneSlipFn::slip(const double t,
- const std::vector<Mesh::point_type>& vertices)
+ const std::set<Mesh::point_type>& vertices)
{ // slip
- typedef std::vector<Mesh::point_type>::const_iterator vert_iterator;
+ typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
assert(0 != _parameters);
assert(!_slipField.isNull());
Modified: short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh 2007-05-01 20:27:04 UTC (rev 6741)
+++ short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh 2007-05-01 21:34:23 UTC (rev 6742)
@@ -89,7 +89,7 @@
*/
void initialize(const ALE::Obj<Mesh>& mesh,
const ALE::Obj<Mesh>& faultMesh,
- const std::vector<Mesh::point_type>& vertices,
+ const std::set<Mesh::point_type>& vertices,
const spatialdata::geocoords::CoordSys* cs);
/** Get slip on fault surface at time t.
@@ -98,7 +98,7 @@
* @param vertices Vertices where slip will be prescribed.
*/
const ALE::Obj<real_section_type>& slip(const double t,
- const std::vector<Mesh::point_type>& vertices);
+ const std::set<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-05-01 20:27:04 UTC (rev 6741)
+++ short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.cc 2007-05-01 21:34:23 UTC (rev 6742)
@@ -55,7 +55,7 @@
pylith::faults::EqKinSrc::initialize(
const ALE::Obj<Mesh>& mesh,
const ALE::Obj<Mesh>& faultMesh,
- const std::vector<Mesh::point_type>& vertices,
+ const std::set<Mesh::point_type>& vertices,
const spatialdata::geocoords::CoordSys* cs)
{ // initialize
assert(0 != _slipfn);
@@ -66,7 +66,7 @@
// Get slip on fault surface at time t.
const ALE::Obj<pylith::real_section_type>&
pylith::faults::EqKinSrc::slip(const double t,
- const std::vector<Mesh::point_type>& vertices)
+ const std::set<Mesh::point_type>& vertices)
{ // slip
assert(0 != _slipfn);
return _slipfn->slip(t, vertices);
Modified: short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh 2007-05-01 20:27:04 UTC (rev 6741)
+++ short/3D/PyLith/trunk/libsrc/faults/EqKinSrc.hh 2007-05-01 21:34:23 UTC (rev 6742)
@@ -83,7 +83,7 @@
virtual
void initialize(const ALE::Obj<Mesh>& mesh,
const ALE::Obj<Mesh>& faultMesh,
- const std::vector<Mesh::point_type>& vertices,
+ const std::set<Mesh::point_type>& vertices,
const spatialdata::geocoords::CoordSys* cs);
/** Get slip on fault surface at time t.
@@ -93,7 +93,7 @@
*/
virtual
const ALE::Obj<real_section_type>& slip(const double t,
- const std::vector<Mesh::point_type>& vertices);
+ const std::set<Mesh::point_type>& vertices);
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2007-05-01 20:27:04 UTC (rev 6741)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2007-05-01 21:34:23 UTC (rev 6742)
@@ -101,6 +101,16 @@
const ALE::Obj<real_section_type>& disp,
const ALE::Obj<Mesh>& mesh) = 0;
+ // PROTECTED TYPEDEFS /////////////////////////////////////////////////
+protected :
+
+ /// Function type for orientation methods.
+ typedef void (*orient_fn_type)(double_array*,
+ const double_array&,
+ const double_array&,
+ const double_array&,
+ const int);
+
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -136,6 +146,7 @@
* be up-dip direction).
* @param numLocs Number of locations where values are given.
*/
+ static
void _orient1D(double_array* orientation,
const double_array& jacobian,
const double_array& jacobianDet,
@@ -160,6 +171,7 @@
* be up-dip direction).
* @param numLocs Number of locations where values are given.
*/
+ static
void _orient2D(double_array* orientation,
const double_array& jacobian,
const double_array& jacobianDet,
@@ -184,6 +196,7 @@
* be up-dip direction).
* @param numLocs Number of locations where values are given.
*/
+ static
void _orient3D(double_array* orientation,
const double_array& jacobian,
const double_array& jacobianDet,
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-05-01 20:27:04 UTC (rev 6741)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2007-05-01 21:34:23 UTC (rev 6742)
@@ -83,12 +83,34 @@
mesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- // Loop over cells
+ // 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) {
@@ -99,22 +121,95 @@
const double_array& jacobianDet = _quadrature->jacobianDetVert();
// Compute weighted orientation of face at vertices (using geometry info)
- // NEED TO CALL APPROPRIATE ORIENT FOR CELL DIM
- //_orient1D(&cellOrientation, jacobian, jacobianDet, upDir, numVertices);
+ orientFn(&cellOrientation, jacobian, jacobianDet, upDir, numVertices);
- // Update weighted orientations
- // Loop over vertices in cell
- // Update section with orientation for each vertex
+ // 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
+ // 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 list of constraint vertices
+ // 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
} // initialize
// ----------------------------------------------------------------------
@@ -125,10 +220,44 @@
const ALE::Obj<real_section_type>& disp,
const ALE::Obj<Mesh>& mesh)
{ // integrateResidual
-
- // Subtract constraint forces (which are in disp at the constraint
+ // 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
} // integrateResidual
// ----------------------------------------------------------------------
@@ -144,6 +273,74 @@
// direction cosines. Entries are associated with vertices ik, jk,
// ki, and kj.
+ 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 matrix for cell values (if necessary)
+ _initCellMatrix();
+
+ const int numVertices = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ for (Mesh::label_sequence::iterator c_iter=cBegin;
+ c_iter != cEnd;
+ ++c_iter) {
+ _resetCellMatrix();
+
+ // Get orientations at cells vertices (only valid at constraint vertices)
+ const real_section_type::value_type* cellOrientation =
+ mesh->restrict(_orientation, *c_iter);
+
+ const int numConstraintVert = numVertices / 3;
+ assert(numVertices == numConstraintVert * 3);
+ const int orientationSize = _orientationSize();
+ for (int iConstraint=0; iConstraint < numConstraintVert; ++iConstraint) {
+ // Get orientations at constraint vertex
+ const real_section_type::value_type* constraintOrient =
+ &cellOrientation[iConstraint*orientationSize];
+
+ // Blocks in cell matrix associated with normal cohesive
+ // vertices i and j and constraint vertex k
+ const int iBasis = 3*iConstraint;
+ const int jBasis = 3*iConstraint + 1;
+ const int kBasis = 3*iConstraint + 2;
+
+ // Entries associated with constraint forces applied at node i
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int kDim=0; kDim < spaceDim; ++kDim) {
+ const int row = iBasis*spaceDim+iDim;
+ const int col = kBasis*spaceDim+kDim;
+ _cellMatrix[row*numVertices*spaceDim+col] =
+ constraintOrient[iDim*spaceDim+kDim];
+ _cellMatrix[col*numVertices*spaceDim+row] =
+ _cellMatrix[row*numVertices*spaceDim+col]; // symmetric
+ } // for
+
+ // Entries associated with constraint forces applied at node j
+ for (int jDim=0; jDim < spaceDim; ++jDim)
+ for (int kDim=0; kDim < spaceDim; ++kDim) {
+ const int row = jBasis*spaceDim+jDim;
+ const int col = kBasis*spaceDim+kDim;
+ _cellMatrix[row*numVertices*spaceDim+col] =
+ -constraintOrient[jDim*spaceDim+kDim];
+ _cellMatrix[col*numVertices*spaceDim+row] =
+ _cellMatrix[row*numVertices*spaceDim+col]; // symmetric
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numConstraintVert*spaceDim*spaceDim*4);
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // Assemble cell contribution into PETSc Matrix
+ const ALE::Obj<Mesh::order_type>& globalOrder =
+ mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
+ err = updateOperator(*mat, mesh, dispT, globalOrder,
+ *c_iter, _cellMatrix, ADD_VALUES);
+ if (err)
+ throw std::runtime_error("Update to PETSc Mat failed.");
+ } // for
} // integrateJacobian
// ----------------------------------------------------------------------
@@ -154,7 +351,7 @@
const ALE::Obj<real_section_type>& disp,
const ALE::Obj<Mesh>& mesh)
{ // setField
- typedef std::vector<Mesh::point_type>::const_iterator vert_iterator;
+ typedef std::set<Mesh::point_type>::const_iterator vert_iterator;
assert(0 != _eqsrc);
@@ -173,20 +370,7 @@
{ // _orientationSize
assert(0 != _quadrature);
- int size = 0;
- switch (_quadrature->cellDim()) {
- case 1 :
- size = 1;
- break;
- case 2 :
- size = 4;
- break;
- case 3 :
- size = 9;
- break;
- default :
- assert(0);
- } // switch
+ const int size = _quadrature->cellDim()*_quadrature->spaceDim();
return size;
} // _orientationSize
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-05-01 20:27:04 UTC (rev 6741)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2007-05-01 21:34:23 UTC (rev 6742)
@@ -151,7 +151,7 @@
ALE::Obj<real_section_type> _orientation;
/// Fault vertices associated with constraints
- std::vector<Mesh::point_type> _constraintVert;
+ std::set<Mesh::point_type> _constraintVert;
}; // class FaultCohesiveKin
Modified: short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh 2007-05-01 20:27:04 UTC (rev 6741)
+++ short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh 2007-05-01 21:34:23 UTC (rev 6742)
@@ -73,7 +73,7 @@
virtual
void initialize(const ALE::Obj<Mesh>& mesh,
const ALE::Obj<Mesh>& faultMesh,
- const std::vector<Mesh::point_type>& vertices,
+ const std::set<Mesh::point_type>& vertices,
const spatialdata::geocoords::CoordSys* cs) = 0;
/** Get slip on fault surface at time t.
@@ -83,7 +83,7 @@
*/
virtual
const ALE::Obj<real_section_type>& slip(const double t,
- const std::vector<Mesh::point_type>& vertices) = 0;
+ const std::set<Mesh::point_type>& vertices) = 0;
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc 2007-05-01 20:27:04 UTC (rev 6741)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc 2007-05-01 21:34:23 UTC (rev 6742)
@@ -336,7 +336,7 @@
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
- // Assemble cell contribution into field
+ // Assemble cell contribution into PETSc Matrix
const ALE::Obj<Mesh::order_type>& globalOrder =
mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
More information about the cig-commits
mailing list