[cig-commits] r16208 - short/3D/PyLith/trunk/libsrc/faults
brad at geodynamics.org
brad at geodynamics.org
Mon Feb 1 18:40:38 PST 2010
Author: brad
Date: 2010-02-01 18:40:38 -0800 (Mon, 01 Feb 2010)
New Revision: 16208
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
Log:
Worked on better data structure for pairing coincident vertices.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2010-02-01 23:06:25 UTC (rev 16207)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2010-02-02 02:40:38 UTC (rev 16208)
@@ -13,23 +13,18 @@
#include <portinfo>
#include "FaultCohesiveKin.hh" // implementation of object methods
-
#include "EqKinSrc.hh" // USES EqKinSrc
#include "CohesiveTopology.hh" // USES CohesiveTopology
-
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
#include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
-
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/topology/Jacobian.hh" // USES Jacobian
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
-
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
-
#include <cmath> // USES pow(), sqrt()
#include <strings.h> // USES strcasecmp()
#include <cstring> // USES strlen()
@@ -37,7 +32,6 @@
#include <cassert> // USES assert()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
-
//#define PRECOMPUTE_GEOMETRY
// ----------------------------------------------------------------------
@@ -47,40 +41,34 @@
// ----------------------------------------------------------------------
// Default constructor.
-pylith::faults::FaultCohesiveKin::FaultCohesiveKin(void)
-{ // constructor
+pylith::faults::FaultCohesiveKin::FaultCohesiveKin(void) { // constructor
} // constructor
// ----------------------------------------------------------------------
// Destructor.
-pylith::faults::FaultCohesiveKin::~FaultCohesiveKin(void)
-{ // destructor
+pylith::faults::FaultCohesiveKin::~FaultCohesiveKin(void) { // destructor
deallocate();
} // destructor
// ----------------------------------------------------------------------
// Deallocate PETSc and local data structures.
-void
-pylith::faults::FaultCohesiveKin::deallocate(void)
-{ // deallocate
+void pylith::faults::FaultCohesiveKin::deallocate(void) { // deallocate
FaultCohesive::deallocate();
// :TODO: Use shared pointers for earthquake sources
} // deallocate
-
+
// ----------------------------------------------------------------------
// Set kinematic earthquake source.
-void
-pylith::faults::FaultCohesiveKin::eqsrcs(const char* const* names,
- const int numNames,
- EqKinSrc** sources,
- const int numSources)
-{ // eqsrcs
+void pylith::faults::FaultCohesiveKin::eqsrcs(const char* const * names,
+ const int numNames,
+ EqKinSrc** sources,
+ const int numSources) { // eqsrcs
assert(numNames == numSources);
// :TODO: Use shared pointers for earthquake sources
_eqSrcs.clear();
- for (int i=0; i < numSources; ++i) {
+ for (int i = 0; i < numSources; ++i) {
if (0 == sources[i])
throw std::runtime_error("Null earthquake source.");
_eqSrcs[std::string(names[i])] = sources[i];
@@ -89,11 +77,9 @@
// ----------------------------------------------------------------------
// Initialize fault. Determine orientation and setup boundary
-void
-pylith::faults::FaultCohesiveKin::initialize(const topology::Mesh& mesh,
- const double upDir[3],
- const double normalDir[3])
-{ // initialize
+void pylith::faults::FaultCohesiveKin::initialize(const topology::Mesh& mesh,
+ const double upDir[3],
+ const double normalDir[3]) { // initialize
assert(0 != upDir);
assert(0 != normalDir);
assert(0 != _quadrature);
@@ -103,18 +89,19 @@
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
assert(0 != cs);
-
- delete _faultMesh; _faultMesh = new topology::SubMesh();
- CohesiveTopology::createFaultParallel(_faultMesh, &_cohesiveToFault,
- mesh, id(), useLagrangeConstraints());
- delete _fields;
- _fields = new topology::Fields<topology::Field<topology::SubMesh> >(*_faultMesh);
+ delete _faultMesh;
+ _faultMesh = new topology::SubMesh();
+ CohesiveTopology::createFaultParallel(_faultMesh, &_cohesiveToFault, mesh, id(), useLagrangeConstraints());
+ //_initializeCohesiveInfo(mesh);
+
+ delete _fields;
+ _fields = new topology::Fields<topology::Field<topology::SubMesh> >(
+ *_faultMesh);
+
const srcs_type::const_iterator srcsEnd = _eqSrcs.end();
- for (srcs_type::iterator s_iter=_eqSrcs.begin();
- s_iter != srcsEnd;
- ++s_iter) {
+ for (srcs_type::iterator s_iter = _eqSrcs.begin(); s_iter != srcsEnd; ++s_iter) {
EqKinSrc* src = s_iter->second;
assert(0 != src);
src->initialize(*_faultMesh, *_normalizer);
@@ -127,7 +114,7 @@
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
+ faultSieveMesh->depthStratum(0);
assert(!vertices.isNull());
_fields->add("slip", "slip");
topology::Field<topology::SubMesh>& slip = _fields->get("slip");
@@ -136,8 +123,8 @@
slip.vectorFieldType(topology::FieldBase::VECTOR);
slip.scale(_normalizer->lengthScale());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ faultSieveMesh->heightStratum(0);
assert(!cells.isNull());
const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
@@ -156,9 +143,8 @@
} // initialize
// ----------------------------------------------------------------------
-void
-pylith::faults::FaultCohesiveKin::splitField(topology::Field<topology::Mesh>* field)
-{ // splitFields
+void pylith::faults::FaultCohesiveKin::splitField(topology::Field<
+ topology::Mesh>* field) { // splitFields
assert(0 != field);
const ALE::Obj<RealSection>& section = field->section();
@@ -177,19 +163,18 @@
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
assert(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin =
- vertices->begin();
+ const SieveSubMesh::label_sequence::iterator verticesBegin =
+ vertices->begin();
const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
- SieveSubMesh::renumbering_type& renumbering =
- faultSieveMesh->getRenumbering();
+ SieveSubMesh::renumbering_type& renumbering =
+ faultSieveMesh->getRenumbering();
const SieveSubMesh::renumbering_type::const_iterator renumberingEnd =
- renumbering.end();
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter)
+ renumbering.end();
+ for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
+ != verticesEnd; ++v_iter)
if (renumbering.find(*v_iter) != renumberingEnd) {
const int vertexFault = renumbering[*v_iter];
const int vertexMesh = *v_iter;
@@ -205,12 +190,10 @@
// ----------------------------------------------------------------------
// Integrate contribution of cohesive cells to residual term that
// require assembly across processors.
-void
-pylith::faults::FaultCohesiveKin::integrateResidual(
- const topology::Field<topology::Mesh>& residual,
- const double t,
- topology::SolutionFields* const fields)
-{ // integrateResidual
+void pylith::faults::FaultCohesiveKin::integrateResidual(const topology::Field<
+ topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields) { // integrateResidual
assert(0 != fields);
assert(0 != _quadrature);
assert(0 != _fields);
@@ -233,20 +216,20 @@
// Get cell information and setup storage for cell data
const int spaceDim = _quadrature->spaceDim();
- const int orientationSize = spaceDim*spaceDim;
+ const int orientationSize = spaceDim * spaceDim;
const int numBasis = _quadrature->numBasis();
const int numConstraintVert = numBasis;
- const int numCorners = 3*numConstraintVert; // cohesive cell
+ const int numCorners = 3 * numConstraintVert; // cohesive cell
const int numQuadPts = _quadrature->numQuadPts();
const double_array& quadWts = _quadrature->quadWts();
assert(quadWts.size() == numQuadPts);
// Allocate vectors for cell values
- double_array orientationCell(numConstraintVert*orientationSize);
- double_array dispTCell(numCorners*spaceDim);
- double_array dispTIncrCell(numCorners*spaceDim);
- double_array dispTpdtCell(numCorners*spaceDim);
- double_array residualCell(numCorners*spaceDim);
+ double_array orientationCell(numConstraintVert * orientationSize);
+ double_array dispTCell(numCorners * spaceDim);
+ double_array dispTIncrCell(numCorners * spaceDim);
+ double_array dispTpdtCell(numCorners * spaceDim);
+ double_array residualCell(numCorners * spaceDim);
// Tributary area for the current for each vertex.
double_array areaVertexCell(numConstraintVert);
@@ -257,13 +240,13 @@
// Get cohesive cells
const ALE::Obj<SieveMesh>& sieveMesh = residual.mesh().sieveMesh();
assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cellsCohesive =
- sieveMesh->getLabelStratum("material-id", id());
+ const ALE::Obj<SieveMesh::label_sequence>& cellsCohesive =
+ sieveMesh->getLabelStratum("material-id", id());
assert(!cellsCohesive.isNull());
const SieveMesh::label_sequence::iterator cellsCohesiveBegin =
- cellsCohesive->begin();
+ cellsCohesive->begin();
const SieveMesh::label_sequence::iterator cellsCohesiveEnd =
- cellsCohesive->end();
+ cellsCohesive->end();
const int cellsCohesiveSize = cellsCohesive->size();
// Get fault Sieve mesh
@@ -271,50 +254,44 @@
assert(!faultSieveMesh.isNull());
// Get section information
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
+ const ALE::Obj<RealSection>& orientationSection =
+ _fields->get("orientation").section();
assert(!orientationSection.isNull());
topology::Mesh::RestrictVisitor orientationVisitor(*orientationSection,
- orientationCell.size(),
- &orientationCell[0]);
+ orientationCell.size(), &orientationCell[0]);
- const ALE::Obj<RealSection>& areaSection =
- _fields->get("area").section();
+ const ALE::Obj<RealSection>& areaSection = _fields->get("area").section();
assert(!areaSection.isNull());
- topology::Mesh::RestrictVisitor areaVisitor(*areaSection,
- areaCell.size(), &areaCell[0]);
+ topology::Mesh::RestrictVisitor areaVisitor(*areaSection, areaCell.size(),
+ &areaCell[0]);
topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
const ALE::Obj<RealSection>& dispTSection = dispT.section();
- assert(!dispTSection.isNull());
- topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
- dispTCell.size(),
- &dispTCell[0]);
+ assert(!dispTSection.isNull());
+ topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection, dispTCell.size(),
+ &dispTCell[0]);
topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
const ALE::Obj<RealSection>& dispTIncrSection = dispTIncr.section();
- assert(!dispTIncrSection.isNull());
+ assert(!dispTIncrSection.isNull());
topology::Mesh::RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
- dispTIncrCell.size(),
- &dispTIncrCell[0]);
+ dispTIncrCell.size(), &dispTIncrCell[0]);
const ALE::Obj<RealSection>& residualSection = residual.section();
topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
- &residualCell[0]);
+ &residualCell[0]);
- double_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- faultSieveMesh->getRealSection("coordinates");
+ double_array coordinatesCell(numBasis * spaceDim);
+ const ALE::Obj<RealSection>& coordinates = faultSieveMesh->getRealSection(
+ "coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
_logger->eventEnd(setupEvent);
- for (SieveMesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
- c_iter != cellsCohesiveEnd;
- ++c_iter) {
+ for (SieveMesh::label_sequence::iterator c_iter = cellsCohesiveBegin; c_iter
+ != cellsCohesiveEnd; ++c_iter) {
const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
areaVertexCell = 0.0;
residualCell = 0.0;
@@ -339,19 +316,19 @@
// Get orientations at fault cell's vertices.
orientationVisitor.clear();
faultSieveMesh->restrictClosure(c_fault, orientationVisitor);
-
+
// Get area at fault cell's vertices.
areaVisitor.clear();
faultSieveMesh->restrictClosure(c_fault, areaVisitor);
-
+
// Get disp(t) at cohesive cell's vertices.
dispTVisitor.clear();
sieveMesh->restrictClosure(*c_iter, dispTVisitor);
-
+
// Get dispIncr(t) at cohesive cell's vertices.
dispTIncrVisitor.clear();
sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
-
+
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -439,12 +416,10 @@
// ----------------------------------------------------------------------
// Integrate contribution of cohesive cells to residual term that do
// not require assembly across cells, vertices, or processors.
-void
-pylith::faults::FaultCohesiveKin::integrateResidualAssembled(
- const topology::Field<topology::Mesh>& residual,
- const double t,
- topology::SolutionFields* const fields)
-{ // integrateResidualAssembled
+void pylith::faults::FaultCohesiveKin::integrateResidualAssembled(const topology::Field<
+ topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields) { // integrateResidualAssembled
assert(0 != fields);
assert(0 != _fields);
assert(0 != _logger);
@@ -466,9 +441,7 @@
slip.zero();
// Compute slip field at current time step
const srcs_type::const_iterator srcsEnd = _eqSrcs.end();
- for (srcs_type::iterator s_iter=_eqSrcs.begin();
- s_iter != srcsEnd;
- ++s_iter) {
+ for (srcs_type::iterator s_iter = _eqSrcs.begin(); s_iter != srcsEnd; ++s_iter) {
EqKinSrc* src = s_iter->second;
assert(0 != src);
if (t >= src->originTime())
@@ -487,22 +460,21 @@
const ALE::Obj<RealSection>& residualSection = residual.section();
assert(!residualSection.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
assert(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin =
- vertices->begin();
+ const SieveSubMesh::label_sequence::iterator verticesBegin =
+ vertices->begin();
const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
- SieveSubMesh::renumbering_type& renumbering =
- faultSieveMesh->getRenumbering();
+ SieveSubMesh::renumbering_type& renumbering =
+ faultSieveMesh->getRenumbering();
const SieveSubMesh::renumbering_type::const_iterator renumberingEnd =
- renumbering.end();
+ renumbering.end();
_logger->eventEnd(setupEvent);
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter)
+ for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
+ != verticesEnd; ++v_iter)
if (renumbering.find(*v_iter) != renumberingEnd) {
const int vertexFault = renumbering[*v_iter];
const int vertexMesh = *v_iter;
@@ -524,12 +496,9 @@
// ----------------------------------------------------------------------
// Compute Jacobian matrix (A) associated with operator that do not
// require assembly across cells, vertices, or processors.
-void
-pylith::faults::FaultCohesiveKin::integrateJacobianAssembled(
- topology::Jacobian* jacobian,
- const double t,
- topology::SolutionFields* const fields)
-{ // integrateJacobianAssembled
+void pylith::faults::FaultCohesiveKin::integrateJacobianAssembled(topology::Jacobian* jacobian,
+ const double t,
+ topology::SolutionFields* const fields) { // integrateJacobianAssembled
assert(0 != jacobian);
assert(0 != fields);
assert(0 != _fields);
@@ -550,34 +519,33 @@
// Get cohesive cells
const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cellsCohesive =
- sieveMesh->getLabelStratum("material-id", id());
+ const ALE::Obj<SieveMesh::label_sequence>& cellsCohesive =
+ sieveMesh->getLabelStratum("material-id", id());
assert(!cellsCohesive.isNull());
const SieveMesh::label_sequence::iterator cellsCohesiveBegin =
- cellsCohesive->begin();
+ cellsCohesive->begin();
const SieveMesh::label_sequence::iterator cellsCohesiveEnd =
- cellsCohesive->end();
+ cellsCohesive->end();
const int cellsCohesiveSize = cellsCohesive->size();
const int spaceDim = _quadrature->spaceDim();
- const int orientationSize = spaceDim*spaceDim;
+ const int orientationSize = spaceDim * spaceDim;
const int numConstraintVert = _quadrature->numBasis();
- const int numCorners = 3*numConstraintVert; // cohesive cell
- double_array matrixCell(numCorners*spaceDim * numCorners*spaceDim);
- double_array orientationCell(numConstraintVert*orientationSize);
+ const int numCorners = 3 * numConstraintVert; // cohesive cell
+ double_array matrixCell(numCorners * spaceDim * numCorners * spaceDim);
+ double_array orientationCell(numConstraintVert * orientationSize);
// Get section information
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
const ALE::Obj<RealSection>& solutionSection = fields->solution().section();
- assert(!solutionSection.isNull());
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
+ assert(!solutionSection.isNull());
+ const ALE::Obj<RealSection>& orientationSection =
+ _fields->get("orientation").section();
assert(!orientationSection.isNull());
topology::Mesh::RestrictVisitor orientationVisitor(*orientationSection,
- orientationCell.size(),
- &orientationCell[0]);
+ orientationCell.size(), &orientationCell[0]);
#if 0 // DEBUGGING
// Check that fault cells match cohesive cells
@@ -587,16 +555,16 @@
const int rank = mesh->commRank();
for (Mesh::label_sequence::iterator c_iter = cellsCohesiveBegin;
- c_iter != cellsCohesiveEnd;
- ++c_iter) {
+ c_iter != cellsCohesiveEnd;
+ ++c_iter) {
mesh->getSieve()->cone(*c_iter, cV);
- const int coneSize = cV.getSize();
- const Mesh::point_type *cone = cV.getPoints();
- const int faceSize = coneSize / 3;
- const Mesh::point_type face = _cohesiveToFault[*c_iter];
+ const int coneSize = cV.getSize();
+ const Mesh::point_type *cone = cV.getPoints();
+ const int faceSize = coneSize / 3;
+ const Mesh::point_type face = _cohesiveToFault[*c_iter];
_faultMesh->getSieve()->cone(face, cV2);
- const int fConeSize = cV2.getSize();
- const Mesh::point_type *fCone = cV2.getPoints();
+ const int fConeSize = cV2.getSize();
+ const Mesh::point_type *fCone = cV2.getPoints();
assert(0 == coneSize % faceSize);
assert(faceSize == fConeSize);
@@ -611,20 +579,19 @@
const PetscMat jacobianMatrix = jacobian->matrix();
assert(0 != jacobianMatrix);
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solutionSection);
+ const ALE::Obj<SieveMesh::order_type>& globalOrder =
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
+ solutionSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
topology::Mesh::IndicesVisitor jacobianVisitor(*solutionSection,
- *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ *globalOrder, (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth()) * spaceDim);
_logger->eventEnd(setupEvent);
- for (SieveMesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
- c_iter != cellsCohesiveEnd;
- ++c_iter) {
+ for (SieveMesh::label_sequence::iterator c_iter = cellsCohesiveBegin; c_iter
+ != cellsCohesiveEnd; ++c_iter) {
const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
_logger->eventBegin(restrictEvent);
@@ -637,39 +604,39 @@
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
- for (int iConstraint=0; iConstraint < numConstraintVert; ++iConstraint) {
+ for (int iConstraint = 0; iConstraint < numConstraintVert; ++iConstraint) {
// Blocks in cell matrix associated with normal cohesive
// vertices i and j and constraint vertex k
const int indexI = iConstraint;
- const int indexJ = iConstraint + numConstraintVert;
- const int indexK = iConstraint + 2*numConstraintVert;
+ const int indexJ = iConstraint + numConstraintVert;
+ const int indexK = iConstraint + 2 * numConstraintVert;
// Get orientation at constraint vertex
- const double* orientationVertex =
- &orientationCell[iConstraint*orientationSize];
+ const double* orientationVertex = &orientationCell[iConstraint
+ * orientationSize];
assert(0 != orientationVertex);
// 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 = indexI*spaceDim+iDim;
- const int col = indexK*spaceDim+kDim;
- matrixCell[row*numCorners*spaceDim+col] =
- -orientationVertex[kDim*spaceDim+iDim];
- matrixCell[col*numCorners*spaceDim+row] =
- -orientationVertex[kDim*spaceDim+iDim];
- } // for
+ for (int iDim = 0; iDim < spaceDim; ++iDim)
+ for (int kDim = 0; kDim < spaceDim; ++kDim) {
+ const int row = indexI * spaceDim + iDim;
+ const int col = indexK * spaceDim + kDim;
+ matrixCell[row * numCorners * spaceDim + col]
+ = -orientationVertex[kDim * spaceDim + iDim];
+ matrixCell[col * numCorners * spaceDim + row]
+ = -orientationVertex[kDim * spaceDim + iDim];
+ } // 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 = indexJ*spaceDim+jDim;
- const int col = indexK*spaceDim+kDim;
- matrixCell[row*numCorners*spaceDim+col] =
- orientationVertex[kDim*spaceDim+jDim];
- matrixCell[col*numCorners*spaceDim+row] =
- orientationVertex[kDim*spaceDim+jDim];
- } // for
+ for (int jDim = 0; jDim < spaceDim; ++jDim)
+ for (int kDim = 0; kDim < spaceDim; ++kDim) {
+ const int row = indexJ * spaceDim + jDim;
+ const int col = indexK * spaceDim + kDim;
+ matrixCell[row * numCorners * spaceDim + col]
+ = orientationVertex[kDim * spaceDim + jDim];
+ matrixCell[col * numCorners * spaceDim + row]
+ = orientationVertex[kDim * spaceDim + jDim];
+ } // for
} // for
_logger->eventEnd(computeEvent);
@@ -678,8 +645,7 @@
_logger->eventBegin(updateEvent);
jacobianVisitor.clear();
PetscErrorCode err = updateOperator(jacobianMatrix, *sieveMesh->getSieve(),
- jacobianVisitor, *c_iter,
- &matrixCell[0], INSERT_VALUES);
+ jacobianVisitor, *c_iter, &matrixCell[0], INSERT_VALUES);
CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
_logger->eventEnd(updateEvent);
@@ -687,16 +653,14 @@
PetscLogFlops(cellsCohesiveSize*numConstraintVert*spaceDim*spaceDim*4);
_needNewJacobian = false;
} // integrateJacobianAssembled
-
+
// ----------------------------------------------------------------------
// Compute Jacobian matrix (A) associated with operator that do not
// require assembly across cells, vertices, or processors.
-void
-pylith::faults::FaultCohesiveKin::integrateJacobianAssembled(
- topology::Field<topology::Mesh>* jacobian,
- const double t,
- topology::SolutionFields* const fields)
-{ // integrateJacobianAssembled
+void pylith::faults::FaultCohesiveKin::integrateJacobianAssembled(topology::Field<
+ topology::Mesh>* jacobian,
+ const double t,
+ topology::SolutionFields* const fields) { // integrateJacobianAssembled
assert(0 != jacobian);
assert(0 != fields);
assert(0 != _fields);
@@ -723,16 +687,16 @@
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
assert(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin =
- vertices->begin();
+ const SieveSubMesh::label_sequence::iterator verticesBegin =
+ vertices->begin();
const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
- SieveSubMesh::renumbering_type& renumbering =
- faultSieveMesh->getRenumbering();
+ SieveSubMesh::renumbering_type& renumbering =
+ faultSieveMesh->getRenumbering();
const SieveSubMesh::renumbering_type::const_iterator renumberingEnd =
- renumbering.end();
+ renumbering.end();
// Get section information
const int spaceDim = _quadrature->spaceDim();
@@ -743,9 +707,8 @@
_logger->eventEnd(setupEvent);
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter)
+ for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
+ != verticesEnd; ++v_iter)
if (renumbering.find(*v_iter) != renumberingEnd) {
assert(jacobianSection->getFiberDimension(*v_iter) == spaceDim);
@@ -758,14 +721,11 @@
_needNewJacobian = false;
} // integrateJacobianAssembled
-
+
// ----------------------------------------------------------------------
// Update state variables as needed.
-void
-pylith::faults::FaultCohesiveKin::updateStateVars(
- const double t,
- topology::SolutionFields* const fields)
-{ // updateStateVars
+void pylith::faults::FaultCohesiveKin::updateStateVars(const double t,
+ topology::SolutionFields* const fields) { // updateStateVars
assert(0 != fields);
assert(0 != _fields);
@@ -774,10 +734,9 @@
// ----------------------------------------------------------------------
// Adjust solution from solver with lumped Jacobian to match Lagrange
// multiplier constraints.
-void
-pylith::faults::FaultCohesiveKin::adjustSolnLumped(topology::SolutionFields* const fields,
- const topology::Field<topology::Mesh>& jacobian)
-{ // adjustSolnLumped
+void pylith::faults::FaultCohesiveKin::adjustSolnLumped(topology::SolutionFields* const fields,
+ const topology::Field<
+ topology::Mesh>& jacobian) { // adjustSolnLumped
assert(0 != fields);
assert(0 != _quadrature);
@@ -802,25 +761,25 @@
// Get cell information and setup storage for cell data
const int spaceDim = _quadrature->spaceDim();
- const int orientationSize = spaceDim*spaceDim;
+ const int orientationSize = spaceDim * spaceDim;
const int numBasis = _quadrature->numBasis();
const int numConstraintVert = numBasis;
- const int numCorners = 3*numConstraintVert; // cohesive cell
+ const int numCorners = 3 * numConstraintVert; // cohesive cell
const int numQuadPts = _quadrature->numQuadPts();
const double_array& quadWts = _quadrature->quadWts();
assert(quadWts.size() == numQuadPts);
// Get cohesive cells
- const ALE::Obj<SieveMesh>& sieveMesh =
- fields->get("dispIncr(t->t+dt)").mesh().sieveMesh();
+ const ALE::Obj<SieveMesh>& sieveMesh =
+ fields->get("dispIncr(t->t+dt)").mesh().sieveMesh();
assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cellsCohesive =
- sieveMesh->getLabelStratum("material-id", id());
+ const ALE::Obj<SieveMesh::label_sequence>& cellsCohesive =
+ sieveMesh->getLabelStratum("material-id", id());
assert(!cellsCohesive.isNull());
const SieveMesh::label_sequence::iterator cellsCohesiveBegin =
- cellsCohesive->begin();
+ cellsCohesive->begin();
const SieveMesh::label_sequence::iterator cellsCohesiveEnd =
- cellsCohesive->end();
+ cellsCohesive->end();
const int cellsCohesiveSize = cellsCohesive->size();
// Get fault Sieve mesh
@@ -828,71 +787,62 @@
assert(!faultSieveMesh.isNull());
// Get section information
- double_array orientationCell(numConstraintVert*orientationSize);
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
+ double_array orientationCell(numConstraintVert * orientationSize);
+ const ALE::Obj<RealSection>& orientationSection =
+ _fields->get("orientation").section();
assert(!orientationSection.isNull());
topology::Mesh::RestrictVisitor orientationVisitor(*orientationSection,
- orientationCell.size(),
- &orientationCell[0]);
+ orientationCell.size(), &orientationCell[0]);
- double_array slipCell(numConstraintVert*spaceDim);
+ double_array slipCell(numConstraintVert * spaceDim);
const ALE::Obj<RealSection>& slipSection = _fields->get("slip").section();
assert(!slipSection.isNull());
- topology::Mesh::RestrictVisitor slipVisitor(*slipSection,
- slipCell.size(),
- &slipCell[0]);
+ topology::Mesh::RestrictVisitor slipVisitor(*slipSection, slipCell.size(),
+ &slipCell[0]);
// Tributary area for the current for each vertex.
double_array areaVertexCell(numConstraintVert);
// Total fault area associated with each vertex (assembled over all cells).
double_array areaCell(numConstraintVert);
- const ALE::Obj<RealSection>& areaSection =
- _fields->get("area").section();
+ const ALE::Obj<RealSection>& areaSection = _fields->get("area").section();
assert(!areaSection.isNull());
- topology::Mesh::RestrictVisitor areaVisitor(*areaSection,
- areaCell.size(), &areaCell[0]);
+ topology::Mesh::RestrictVisitor areaVisitor(*areaSection, areaCell.size(),
+ &areaCell[0]);
- double_array jacobianCell(numCorners*spaceDim);
+ double_array jacobianCell(numCorners * spaceDim);
const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
- assert(!jacobianSection.isNull());
+ assert(!jacobianSection.isNull());
topology::Mesh::RestrictVisitor jacobianVisitor(*jacobianSection,
- jacobianCell.size(),
- &jacobianCell[0]);
+ jacobianCell.size(), &jacobianCell[0]);
- double_array residualCell(numCorners*spaceDim);
- const ALE::Obj<RealSection>& residualSection =
- fields->get("residual").section();
+ double_array residualCell(numCorners * spaceDim);
+ const ALE::Obj<RealSection>& residualSection =
+ fields->get("residual").section();
topology::Mesh::RestrictVisitor residualVisitor(*residualSection,
- residualCell.size(),
- &residualCell[0]);
+ residualCell.size(), &residualCell[0]);
- double_array dispTCell(numCorners*spaceDim);
- const ALE::Obj<RealSection>& dispTSection =
- fields->get("disp(t)").section();
- topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
- dispTCell.size(),
- &dispTCell[0]);
+ double_array dispTCell(numCorners * spaceDim);
+ const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
+ topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection, dispTCell.size(),
+ &dispTCell[0]);
- double_array solutionCell(numCorners*spaceDim);
- const ALE::Obj<RealSection>& solutionSection =
- fields->get("dispIncr(t->t+dt)").section();
+ double_array solutionCell(numCorners * spaceDim);
+ const ALE::Obj<RealSection>& solutionSection = fields->get(
+ "dispIncr(t->t+dt)").section();
topology::Mesh::UpdateAddVisitor solutionVisitor(*solutionSection,
- &solutionCell[0]);
+ &solutionCell[0]);
- double_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- faultSieveMesh->getRealSection("coordinates");
+ double_array coordinatesCell(numBasis * spaceDim);
+ const ALE::Obj<RealSection>& coordinates = faultSieveMesh->getRealSection(
+ "coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
_logger->eventEnd(setupEvent);
- for (SieveMesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
- c_iter != cellsCohesiveEnd;
- ++c_iter) {
+ for (SieveMesh::label_sequence::iterator c_iter = cellsCohesiveBegin; c_iter
+ != cellsCohesiveEnd; ++c_iter) {
const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
areaVertexCell = 0.0;
solutionCell = 0.0;
@@ -917,74 +867,74 @@
// Get orientations at fault cell's vertices.
orientationVisitor.clear();
faultSieveMesh->restrictClosure(c_fault, orientationVisitor);
-
+
// Get area at fault cell's vertices.
areaVisitor.clear();
faultSieveMesh->restrictClosure(c_fault, areaVisitor);
-
+
// Get slip at fault cell's vertices.
slipVisitor.clear();
faultSieveMesh->restrictClosure(c_fault, slipVisitor);
-
+
// Get residual at cohesive cell's vertices.
residualVisitor.clear();
sieveMesh->restrictClosure(*c_iter, residualVisitor);
-
+
// Get jacobian at cohesive cell's vertices.
jacobianVisitor.clear();
sieveMesh->restrictClosure(*c_iter, jacobianVisitor);
-
+
// Get disp(t) at cohesive cell's vertices.
dispTVisitor.clear();
sieveMesh->restrictClosure(*c_iter, dispTVisitor);
-
+
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
// Compute contributory area for cell (to weight contributions)
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ for (int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- const double dArea = wt*basis[iQuad*numBasis+iBasis];
+ for (int iBasis = 0; iBasis < numBasis; ++iBasis) {
+ const double dArea = wt * basis[iQuad * numBasis + iBasis];
areaVertexCell[iBasis] += dArea;
} // for
} // for
- for (int iConstraint=0; iConstraint < numConstraintVert; ++iConstraint) {
+ for (int iConstraint = 0; iConstraint < numConstraintVert; ++iConstraint) {
// Blocks in cell matrix associated with normal cohesive
// vertices i and j and constraint vertex k
const int indexI = iConstraint;
- const int indexJ = iConstraint + numConstraintVert;
- const int indexK = iConstraint + 2*numConstraintVert;
+ const int indexJ = iConstraint + numConstraintVert;
+ const int indexK = iConstraint + 2 * numConstraintVert;
assert(areaCell[iConstraint] > 0);
const double wt = areaVertexCell[iConstraint] / areaCell[iConstraint];
-
+
// Get orientation at constraint vertex
- const double* orientationVertex =
- &orientationCell[iConstraint*orientationSize];
+ const double* orientationVertex = &orientationCell[iConstraint
+ * orientationSize];
assert(0 != orientationVertex);
// Get slip at constraint vertex
- const double* slipVertex = &slipCell[iConstraint*spaceDim];
+ const double* slipVertex = &slipCell[iConstraint * spaceDim];
assert(0 != slipVertex);
// Get Jacobian at conventional vertices i and j
- const double* Ai = &jacobianCell[indexI*spaceDim];
+ const double* Ai = &jacobianCell[indexI * spaceDim];
assert(0 != Ai);
- const double* Aj = &jacobianCell[indexJ*spaceDim];
+ const double* Aj = &jacobianCell[indexJ * spaceDim];
assert(0 != Aj);
// Get residual at conventional vertices i and j
- const double* ri = &residualCell[indexI*spaceDim];
+ const double* ri = &residualCell[indexI * spaceDim];
assert(0 != ri);
- const double* rj = &residualCell[indexJ*spaceDim];
+ const double* rj = &residualCell[indexJ * spaceDim];
assert(0 != rj);
-
+
// Get disp(t) at conventional vertices i and j
- const double* ui = &dispTCell[indexI*spaceDim];
+ const double* ui = &dispTCell[indexI * spaceDim];
assert(0 != ui);
- const double* uj = &dispTCell[indexJ*spaceDim];
+ const double* uj = &dispTCell[indexJ * spaceDim];
assert(0 != uj);
switch (spaceDim) { // switch
@@ -1129,7 +1079,7 @@
throw std::logic_error("Unknown spatial dimension.");
} // switch
} // for
-
+
_logger->eventEnd(computeEvent);
#if 0 // DEBUGGING
@@ -1160,10 +1110,7 @@
// ----------------------------------------------------------------------
// Verify configuration is acceptable.
-void
-pylith::faults::FaultCohesiveKin::verifyConfiguration(
- const topology::Mesh& mesh) const
-{ // verifyConfiguration
+void pylith::faults::FaultCohesiveKin::verifyConfiguration(const topology::Mesh& mesh) const { // verifyConfiguration
assert(0 != _quadrature);
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
@@ -1172,38 +1119,36 @@
if (!sieveMesh->hasIntSection(label())) {
std::ostringstream msg;
msg << "Mesh missing group of vertices '" << label()
- << " for boundary condition.";
+ << " for boundary condition.";
throw std::runtime_error(msg.str());
} // if
// check compatibility of mesh and quadrature scheme
- const int dimension = mesh.dimension()-1;
+ const int dimension = mesh.dimension() - 1;
if (_quadrature->cellDim() != dimension) {
std::ostringstream msg;
- msg << "Dimension of reference cell in quadrature scheme ("
- << _quadrature->cellDim()
- << ") does not match dimension of cells in mesh ("
- << dimension << ") for fault '" << label()
- << "'.";
+ msg << "Dimension of reference cell in quadrature scheme ("
+ << _quadrature->cellDim()
+ << ") does not match dimension of cells in mesh (" << dimension
+ << ") for fault '" << label() << "'.";
throw std::runtime_error(msg.str());
} // if
const int numCorners = _quadrature->refGeometry().numCorners();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", id());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", id());
assert(!cells.isNull());
const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for (SieveMesh::label_sequence::iterator c_iter = cellsBegin; c_iter
+ != cellsEnd; ++c_iter) {
const int cellNumCorners = sieveMesh->getNumCellCorners(*c_iter);
- if (3*numCorners != cellNumCorners) {
+ if (3 * numCorners != cellNumCorners) {
std::ostringstream msg;
- msg << "Number of vertices in reference cell (" << numCorners
- << ") is not compatible with number of vertices (" << cellNumCorners
- << ") in cohesive cell " << *c_iter << " for fault '"
- << label() << "'.";
+ msg << "Number of vertices in reference cell (" << numCorners
+ << ") is not compatible with number of vertices (" << cellNumCorners
+ << ") in cohesive cell " << *c_iter << " for fault '" << label()
+ << "'.";
throw std::runtime_error(msg.str());
} // if
} // for
@@ -1212,10 +1157,8 @@
// ----------------------------------------------------------------------
// Get vertex field associated with integrator.
const pylith::topology::Field<pylith::topology::SubMesh>&
-pylith::faults::FaultCohesiveKin::vertexField(
- const char* name,
- const topology::SolutionFields* fields)
-{ // vertexField
+pylith::faults::FaultCohesiveKin::vertexField(const char* name,
+ const topology::SolutionFields* fields) { // vertexField
assert(0 != _faultMesh);
assert(0 != _quadrature);
assert(0 != _normalizer);
@@ -1230,65 +1173,63 @@
double scale = 0.0;
int fiberDim = 0;
if (0 == strcasecmp("slip", name)) {
- const topology::Field<topology::SubMesh>& slip =
- _fields->get("slip");
+ const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
return slip;
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
+ const ALE::Obj<RealSection>& orientationSection = _fields->get(
+ "orientation").section();
assert(!orientationSection.isNull());
- const ALE::Obj<RealSection>& dirSection =
- orientationSection->getFibration(0);
+ const ALE::Obj<RealSection>& dirSection = orientationSection->getFibration(
+ 0);
assert(!dirSection.isNull());
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
+ _fields->get("buffer (vector)");
buffer.copy(dirSection);
buffer.label("strike_dir");
buffer.scale(1.0);
return buffer;
} else if (2 == cohesiveDim && 0 == strcasecmp("dip_dir", name)) {
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
+ const ALE::Obj<RealSection>& orientationSection = _fields->get(
+ "orientation").section();
assert(!orientationSection.isNull());
- const ALE::Obj<RealSection>& dirSection =
- orientationSection->getFibration(1);
+ const ALE::Obj<RealSection>& dirSection = orientationSection->getFibration(
+ 1);
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
+ _fields->get("buffer (vector)");
buffer.copy(dirSection);
buffer.label("dip_dir");
buffer.scale(1.0);
return buffer;
} else if (0 == strcasecmp("normal_dir", name)) {
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
+ const ALE::Obj<RealSection>& orientationSection = _fields->get(
+ "orientation").section();
assert(!orientationSection.isNull());
- const int space =
- (0 == cohesiveDim) ? 0 : (1 == cohesiveDim) ? 1 : 2;
- const ALE::Obj<RealSection>& dirSection =
- orientationSection->getFibration(space);
+ const int space = (0 == cohesiveDim) ? 0 : (1 == cohesiveDim) ? 1 : 2;
+ const ALE::Obj<RealSection>& dirSection = orientationSection->getFibration(
+ space);
assert(!dirSection.isNull());
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
+ _fields->get("buffer (vector)");
buffer.copy(dirSection);
buffer.label("normal_dir");
buffer.scale(1.0);
return buffer;
} else if (0 == strncasecmp("final_slip_X", name, slipStrLen)) {
- const std::string value = std::string(name).substr(slipStrLen+1);
+ const std::string value = std::string(name).substr(slipStrLen + 1);
const srcs_type::const_iterator s_iter = _eqSrcs.find(value);
assert(s_iter != _eqSrcs.end());
return s_iter->second->finalSlip();
} else if (0 == strncasecmp("slip_time_X", name, timeStrLen)) {
- const std::string value = std::string(name).substr(timeStrLen+1);
+ const std::string value = std::string(name).substr(timeStrLen + 1);
const srcs_type::const_iterator s_iter = _eqSrcs.find(value);
assert(s_iter != _eqSrcs.end());
return s_iter->second->slipTime();
@@ -1298,51 +1239,117 @@
const topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
+ _fields->get("buffer (vector)");
_calcTractionsChange(&buffer, dispT);
return buffer;
} else {
std::ostringstream msg;
- msg << "Request for unknown vertex field '" << name
- << "' for fault '" << label() << "'.";
+ msg << "Request for unknown vertex field '" << name << "' for fault '"
+ << label() << "'.";
throw std::runtime_error(msg.str());
} // else
-
+
// Satisfy return values
assert(0 != _fields);
- const topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
+ const topology::Field<topology::SubMesh>& buffer = _fields->get(
+ "buffer (vector)");
return buffer;
} // vertexField
// ----------------------------------------------------------------------
// Get cell field associated with integrator.
const pylith::topology::Field<pylith::topology::SubMesh>&
-pylith::faults::FaultCohesiveKin::cellField(
- const char* name,
- const topology::SolutionFields* fields)
-{ // cellField
+pylith::faults::FaultCohesiveKin::cellField(const char* name,
+ const topology::SolutionFields* fields) { // cellField
// Should not reach this point if requested field was found
std::ostringstream msg;
- msg << "Request for unknown cell field '" << name
- << "' for fault '" << label() << ".";
+ msg << "Request for unknown cell field '" << name << "' for fault '"
+ << label() << ".";
throw std::runtime_error(msg.str());
// Satisfy return values
assert(0 != _fields);
- const topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
+ const topology::Field<topology::SubMesh>& buffer = _fields->get(
+ "buffer (vector)");
return buffer;
} // cellField
// ----------------------------------------------------------------------
+// Initialize auxiliary cohesive cell information.
+void pylith::faults::FaultCohesiveKin::_initializeCohesiveInfo(const topology::Mesh& mesh) { // _initializeCohesiveInfo
+ assert(0 != _quadrature);
+
+ // Get cohesive cells
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", id());
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ const int numConstraintVert = _quadrature->numBasis();
+ const int numCorners = 3 * numConstraintVert; // cohesive cell
+
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+ SubMesh::renumbering_type& renumbering = faultSieveMesh->getRenumbering();
+ const SieveSubMesh::renumbering_type::const_iterator renumberingEnd =
+ renumbering.end();
+
+ typedef std::map<int, int> indexmap_type;
+ indexmap_type indexMap;
+ _cohesiveVertices.resize(renumbering.size());
+ int index = 0;
+
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = mesh.sieveMesh()->getSieve();
+ assert(!sieve.isNull());
+ typedef ALE::SieveAlg<SieveMesh> SieveAlg;
+
+ ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type> ncV(*sieve,
+ (size_t) pow(sieve->getMaxConeSize(), std::max(0, sieveMesh->depth())));
+
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Get oriented closure
+ ncV.clear();
+ ALE::ISieveTraversal<SieveMesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
+ const int coneSize = ncV.getSize();
+ assert(coneSize == numCorners);
+ const Mesh::point_type *cone = ncV.getPoints();
+ assert(0 != cone);
+
+ for (int iConstraint = 0; iConstraint < numConstraintVert; ++iConstraint) {
+ // normal cohesive vertices i and j and constraint vertex k
+ const int indexI = iConstraint;
+ const int indexJ = iConstraint + numConstraintVert;
+ const int indexK = iConstraint + 2 * numConstraintVert;
+
+ const int v_lagrange = cone[indexK];
+ const int v_negative = cone[indexI];
+ const int v_positive = cone[indexJ];
+ const int v_fault = renumbering[v_lagrange];
+
+ if (indexMap.end() == indexMap.find(v_lagrange)) {
+ _cohesiveVertices[index].lagrange = v_lagrange;
+ _cohesiveVertices[index].positive = v_positive;
+ _cohesiveVertices[index].negative = v_negative;
+ _cohesiveVertices[index].fault = v_fault;
+ indexMap[v_lagrange] = index; // add index to map
+ ++index;
+ } // if
+ } // for
+ } // for
+} // _initializeCohesiveInfo
+
+// ----------------------------------------------------------------------
// Initialize logger.
-void
-pylith::faults::FaultCohesiveKin::_initializeLogger(void)
-{ // initializeLogger
- delete _logger; _logger = new utils::EventLogger;
+void pylith::faults::FaultCohesiveKin::_initializeLogger(void) { // initializeLogger
+ delete _logger;
+ _logger = new utils::EventLogger;
assert(0 != _logger);
_logger->className("FaultCohesiveKin");
_logger->initialize();
@@ -1358,7 +1365,7 @@
_logger->registerEvent("FkIR compute");
_logger->registerEvent("FkIR restrict");
_logger->registerEvent("FkIR update");
-
+
_logger->registerEvent("FkIJ setup");
_logger->registerEvent("FkIJ geometry");
_logger->registerEvent("FkIJ compute");
@@ -1368,10 +1375,8 @@
// ----------------------------------------------------------------------
// Calculate orientation at fault vertices.
-void
-pylith::faults::FaultCohesiveKin::_calcOrientation(const double upDir[3],
- const double normalDir[3])
-{ // _calcOrientation
+void pylith::faults::FaultCohesiveKin::_calcOrientation(const double upDir[3],
+ const double normalDir[3]) { // _calcOrientation
assert(0 != upDir);
assert(0 != normalDir);
assert(0 != _faultMesh);
@@ -1382,16 +1387,17 @@
// Get vertices in fault mesh.
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
+ faultSieveMesh->depthStratum(0);
+ const SieveSubMesh::label_sequence::iterator verticesBegin =
+ vertices->begin();
const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
-
+
// Containers for orientation information.
const int cohesiveDim = _faultMesh->dimension();
const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
- const int orientationSize = spaceDim*spaceDim;
+ const int orientationSize = spaceDim * spaceDim;
const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
const double_array& verticesRef = cellGeometry.vertices();
const int jacobianSize = (cohesiveDim > 0) ? spaceDim * cohesiveDim : 1;
@@ -1399,7 +1405,7 @@
double_array jacobian(jacobianSize);
double jacobianDet = 0;
double_array orientationVertex(orientationSize);
- double_array coordinatesCell(numBasis*spaceDim);
+ double_array coordinatesCell(numBasis * spaceDim);
double_array refCoordsVertex(cohesiveDim);
// Allocate orientation field.
@@ -1410,16 +1416,16 @@
const ALE::Obj<RealSection>& orientationSection = orientation.section();
assert(!orientationSection.isNull());
// Create subspaces for along-strike, up-dip, and normal directions
- for (int iDim=0; iDim <= cohesiveDim; ++iDim)
+ for (int iDim = 0; iDim <= cohesiveDim; ++iDim)
orientationSection->addSpace();
- for (int iDim=0; iDim <= cohesiveDim; ++iDim)
+ for (int iDim = 0; iDim <= cohesiveDim; ++iDim)
orientationSection->setFiberDimension(vertices, spaceDim, iDim);
orientation.allocate();
orientation.zero();
-
+
// Get fault cells.
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ faultSieveMesh->heightStratum(0);
assert(!cells.isNull());
const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
@@ -1427,12 +1433,11 @@
// Compute orientation of fault at constraint vertices
// Get section containing coordinates of vertices
- const ALE::Obj<RealSection>& coordinatesSection =
- faultSieveMesh->getRealSection("coordinates");
+ const ALE::Obj<RealSection>& coordinatesSection =
+ faultSieveMesh->getRealSection("coordinates");
assert(!coordinatesSection.isNull());
topology::Mesh::RestrictVisitor coordinatesVisitor(*coordinatesSection,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ coordinatesCell.size(), &coordinatesCell[0]);
// Set orientation function
assert(cohesiveDim == _quadrature->cellDim());
@@ -1440,36 +1445,38 @@
// Loop over cohesive cells, computing orientation weighted by
// jacobian at constraint vertices
-
+
const ALE::Obj<SieveSubMesh::sieve_type>& sieve = faultSieveMesh->getSieve();
assert(!sieve.isNull());
typedef ALE::SieveAlg<SieveSubMesh> SieveAlg;
- ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type> ncV(*sieve, (size_t) pow(sieve->getMaxConeSize(), std::max(0, faultSieveMesh->depth())));
+ ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type>
+ ncV(*sieve, (size_t) pow(sieve->getMaxConeSize(), std::max(0,
+ faultSieveMesh->depth())));
- for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin; c_iter
+ != cellsEnd; ++c_iter) {
// Get orientations at fault cell's vertices.
coordinatesVisitor.clear();
faultSieveMesh->restrictClosure(*c_iter, coordinatesVisitor);
ncV.clear();
- ALE::ISieveTraversal<SieveSubMesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
- const int coneSize = ncV.getSize();
- const Mesh::point_type *cone = ncV.getPoints();
-
- for (int v=0; v < coneSize; ++v) {
+ ALE::ISieveTraversal<SieveSubMesh::sieve_type>::orientedClosure(*sieve,
+ *c_iter, ncV);
+ const int coneSize = ncV.getSize();
+ const Mesh::point_type *cone = ncV.getPoints();
+
+ for (int v = 0; v < coneSize; ++v) {
// Compute Jacobian and determinant of Jacobian at vertex
- memcpy(&refCoordsVertex[0], &verticesRef[v*cohesiveDim],
- cohesiveDim*sizeof(double));
+ memcpy(&refCoordsVertex[0], &verticesRef[v * cohesiveDim], cohesiveDim
+ * sizeof(double));
cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell,
- refCoordsVertex);
+ refCoordsVertex);
// Compute orientation
- cellGeometry.orientation(&orientationVertex, jacobian, jacobianDet,
- upDirArray);
-
+ cellGeometry.orientation(&orientationVertex, jacobian, jacobianDet,
+ upDirArray);
+
// Update orientation
orientationSection->updateAddPoint(cone[v], &orientationVertex[0]);
} // for
@@ -1483,20 +1490,19 @@
// Loop over vertices, make orientation information unit magnitude
double_array vertexDir(orientationSize);
int count = 0;
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++count) {
+ for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
+ != verticesEnd; ++v_iter, ++count) {
orientationVertex = 0.0;
orientationSection->restrictPoint(*v_iter, &orientationVertex[0],
- orientationVertex.size());
- for (int iDim=0; iDim < spaceDim; ++iDim) {
+ orientationVertex.size());
+ for (int iDim = 0; iDim < spaceDim; ++iDim) {
double mag = 0;
- for (int jDim=0, index=iDim*spaceDim; jDim < spaceDim; ++jDim)
- mag += pow(orientationVertex[index+jDim],2);
+ for (int jDim = 0, index = iDim * spaceDim; jDim < spaceDim; ++jDim)
+ mag += pow(orientationVertex[index + jDim], 2);
mag = sqrt(mag);
assert(mag > 0.0);
- for (int jDim=0, index=iDim*spaceDim; jDim < spaceDim; ++jDim)
- orientationVertex[index+jDim] /= mag;
+ for (int jDim = 0, index = iDim * spaceDim; jDim < spaceDim; ++jDim)
+ orientationVertex[index + jDim] /= mag;
} // for
orientationSection->updatePoint(*v_iter, &orientationVertex[0]);
@@ -1519,34 +1525,31 @@
// are used are the opposite of what we would want, but we cannot
// flip the fault normal direction because it is tied to how the
// cohesive cells are created.
-
+
assert(vertices->size() > 0);
- orientationSection->restrictPoint(*vertices->begin(), &orientationVertex[0],
- orientationVertex.size());
-
+ orientationSection->restrictPoint(*vertices->begin(),
+ &orientationVertex[0], orientationVertex.size());
+
assert(3 == spaceDim);
double_array normalDirVertex(&orientationVertex[6], 3);
- const double normalDot =
- normalDir[0]*normalDirVertex[0] +
- normalDir[1]*normalDirVertex[1] +
- normalDir[2]*normalDirVertex[2];
-
+ const double normalDot = normalDir[0] * normalDirVertex[0] + normalDir[1]
+ * normalDirVertex[1] + normalDir[2] * normalDirVertex[2];
+
const int istrike = 0;
const int idip = 3;
const int inormal = 6;
if (normalDot < 0.0) {
// Flip dip direction
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
- orientationSection->restrictPoint(*v_iter, &orientationVertex[0],
- orientationVertex.size());
- assert(9 == orientationSection->getFiberDimension(*v_iter));
- for (int iDim=0; iDim < 3; ++iDim) // flip dip
- orientationVertex[idip+iDim] *= -1.0;
-
- // Update direction
- orientationSection->updatePoint(*v_iter, &orientationVertex[0]);
+ for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
+ != verticesEnd; ++v_iter) {
+ orientationSection->restrictPoint(*v_iter, &orientationVertex[0],
+ orientationVertex.size());
+ assert(9 == orientationSection->getFiberDimension(*v_iter));
+ for (int iDim = 0; iDim < 3; ++iDim) // flip dip
+ orientationVertex[idip + iDim] *= -1.0;
+
+ // Update direction
+ orientationSection->updatePoint(*v_iter, &orientationVertex[0]);
} // for
PetscLogFlops(5 + count * 3);
} // if
@@ -1556,9 +1559,7 @@
} // _calcOrientation
// ----------------------------------------------------------------------
-void
-pylith::faults::FaultCohesiveKin::_calcArea(void)
-{ // _calcArea
+void pylith::faults::FaultCohesiveKin::_calcArea(void) { // _calcArea
assert(0 != _faultMesh);
assert(0 != _fields);
@@ -1576,11 +1577,12 @@
// Get vertices in fault mesh.
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
+ faultSieveMesh->depthStratum(0);
+ const SieveSubMesh::label_sequence::iterator verticesBegin =
+ vertices->begin();
const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
-
+
// Allocate area field.
_fields->add("area", "area");
@@ -1591,28 +1593,26 @@
area.zero();
const ALE::Obj<RealSection>& areaSection = area.section();
assert(!areaSection.isNull());
- topology::Mesh::UpdateAddVisitor areaVisitor(*areaSection, &areaCell[0]);
-
- double_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- faultSieveMesh->getRealSection("coordinates");
+ topology::Mesh::UpdateAddVisitor areaVisitor(*areaSection, &areaCell[0]);
+
+ double_array coordinatesCell(numBasis * spaceDim);
+ const ALE::Obj<RealSection>& coordinates = faultSieveMesh->getRealSection(
+ "coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ faultSieveMesh->heightStratum(0);
assert(!cells.isNull());
const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
// Loop over cells in fault mesh, compute area
- for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin; c_iter
+ != cellsEnd; ++c_iter) {
areaCell = 0.0;
-
+
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
@@ -1627,11 +1627,11 @@
const double_array& jacobianDet = _quadrature->jacobianDet();
// Compute area
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ for (int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- const double dArea = wt*basis[iQuad*numBasis+iBasis];
- areaCell[iBasis] += dArea;
+ for (int iBasis = 0; iBasis < numBasis; ++iBasis) {
+ const double dArea = wt * basis[iQuad * numBasis + iBasis];
+ areaCell[iBasis] += dArea;
} // for
} // for
areaVisitor.clear();
@@ -1653,11 +1653,8 @@
// ----------------------------------------------------------------------
// Compute change in tractions on fault surface using solution.
// NOTE: We must convert vertex labels to fault vertex labels
-void
-pylith::faults::FaultCohesiveKin::_calcTractionsChange(
- topology::Field<topology::SubMesh>* tractions,
- const topology::Field<topology::Mesh>& dispT)
-{ // _calcTractionsChange
+void pylith::faults::FaultCohesiveKin::_calcTractionsChange(topology::Field<
+ topology::SubMesh>* tractions, const topology::Field<topology::Mesh>& dispT) { // _calcTractionsChange
assert(0 != tractions);
assert(0 != _faultMesh);
assert(0 != _fields);
@@ -1669,8 +1666,8 @@
// Get vertices from mesh of domain.
const ALE::Obj<SieveMesh>& sieveMesh = dispT.mesh().sieveMesh();
assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
assert(!vertices.isNull());
const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
@@ -1678,38 +1675,38 @@
// Get fault vertices
const ALE::Obj<SieveMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& fvertices =
- faultSieveMesh->depthStratum(0);
- const SieveSubMesh::label_sequence::iterator fverticesBegin = fvertices->begin();
+ const ALE::Obj<SieveSubMesh::label_sequence>& fvertices =
+ faultSieveMesh->depthStratum(0);
+ const SieveSubMesh::label_sequence::iterator fverticesBegin =
+ fvertices->begin();
const SieveSubMesh::label_sequence::iterator fverticesEnd = fvertices->end();
// Get sections.
- const ALE::Obj<RealSection>& areaSection =
- _fields->get("area").section();
+ const ALE::Obj<RealSection>& areaSection = _fields->get("area").section();
assert(!areaSection.isNull());
const ALE::Obj<RealSection>& dispTSection = dispT.section();
- assert(!dispTSection.isNull());
+ assert(!dispTSection.isNull());
const int numFaultVertices = fvertices->size();
Mesh::renumbering_type& renumbering = faultSieveMesh->getRenumbering();
const SieveSubMesh::renumbering_type::const_iterator renumberingEnd =
- renumbering.end();
+ renumbering.end();
#if 0 // DEBUGGING, MOVE TO SEPARATE CHECK METHOD
// Check fault mesh and volume mesh coordinates
- const ALE::Obj<RealSection>& coordinates = mesh->getRealSection("coordinates");
+ const ALE::Obj<RealSection>& coordinates = mesh->getRealSection("coordinates");
const ALE::Obj<RealSection>& fCoordinates = _faultMesh->getRealSection("coordinates");
for (Mesh::label_sequence::iterator v_iter = verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
+ v_iter != verticesEnd;
+ ++v_iter) {
if (renumbering.find(*v_iter) != renumberingEnd) {
- const int v = *v_iter;
- const int dim = coordinates->getFiberDimension(*v_iter);
- const double *a = coordinates->restrictPoint(*v_iter);
- const int fv = renumbering[*v_iter];
- const int fDim = fCoordinates->getFiberDimension(fv);
- const double *fa = fCoordinates->restrictPoint(fv);
+ const int v = *v_iter;
+ const int dim = coordinates->getFiberDimension(*v_iter);
+ const double *a = coordinates->restrictPoint(*v_iter);
+ const int fv = renumbering[*v_iter];
+ const int fDim = fCoordinates->getFiberDimension(fv);
+ const double *fa = fCoordinates->restrictPoint(fv);
if (dim != fDim) throw ALE::Exception("Coordinate fiber dimensions do not match");
for(int d = 0; d < dim; ++d) {
@@ -1729,7 +1726,7 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
//logger.stagePush("Fault");
- const topology::Field<topology::SubMesh>& slip =_fields->get("slip");
+ const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
tractions->newSection(slip, fiberDim);
tractions->allocate();
@@ -1737,11 +1734,10 @@
} // if
assert(!tractionsSection.isNull());
tractions->zero();
-
+
const double pressureScale = tractions->scale();
- for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
- v_iter != verticesEnd;
- ++v_iter)
+ for (SieveMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
+ != verticesEnd; ++v_iter)
if (renumbering.find(*v_iter) != renumberingEnd) {
const int vertexMesh = *v_iter;
const int vertexFault = renumbering[*v_iter];
@@ -1754,8 +1750,8 @@
const double* areaVertex = areaSection->restrictPoint(vertexFault);
assert(0 != areaVertex);
- for (int i=0; i < fiberDim; ++i)
- tractionsVertex[i] = dispTVertex[i] / areaVertex[0];
+ for (int i = 0; i < fiberDim; ++i)
+ tractionsVertex[i] = dispTVertex[i] / areaVertex[0];
tractionsSection->updatePoint(vertexFault, &tractionsVertex[0]);
} // if
@@ -1769,9 +1765,7 @@
// ----------------------------------------------------------------------
// Allocate buffer for vector field.
-void
-pylith::faults::FaultCohesiveKin::_allocateBufferVectorField(void)
-{ // _allocateBufferVectorField
+void pylith::faults::FaultCohesiveKin::_allocateBufferVectorField(void) { // _allocateBufferVectorField
assert(0 != _fields);
if (_fields->hasField("buffer (vector)"))
return;
@@ -1782,10 +1776,8 @@
// Create vector field; use same shape/chart as slip field.
assert(0 != _faultMesh);
_fields->add("buffer (vector)", "buffer");
- topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
- const topology::Field<topology::SubMesh>& slip =
- _fields->get("slip");
+ topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
+ const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
buffer.cloneSection(slip);
buffer.zero();
@@ -1794,9 +1786,7 @@
// ----------------------------------------------------------------------
// Allocate buffer for scalar field.
-void
-pylith::faults::FaultCohesiveKin::_allocateBufferScalarField(void)
-{ // _allocateBufferScalarField
+void pylith::faults::FaultCohesiveKin::_allocateBufferScalarField(void) { // _allocateBufferScalarField
assert(0 != _fields);
if (_fields->hasField("buffer (scalar)"))
return;
@@ -1807,8 +1797,7 @@
// Create vector field; use same shape/chart as area field.
assert(0 != _faultMesh);
_fields->add("buffer (scalar)", "buffer");
- topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (scalar)");
+ topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (scalar)");
const topology::Field<topology::SubMesh>& area = _fields->get("area");
buffer.cloneSection(area);
buffer.zero();
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2010-02-01 23:06:25 UTC (rev 16207)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2010-02-02 02:40:38 UTC (rev 16208)
@@ -243,6 +243,34 @@
const topology::Fields<topology::Field<topology::SubMesh> >*
fields(void) const;
+ // PROTECTED STRUCTS //////////////////////////////////////////////////
+protected :
+
+ /** Data structure to hold relations between vertices in cohesive cell
+ * in mesh of domain and cell of fault mesh.
+ */
+ struct CohesiveInfo {
+ int lagrange; ///< Vertex associated with Lagrange multiplier.
+ int positive; ///< Vertex on positive side of the fault.
+ int negative; ///< Vertex on negative side of the fault.
+ int fault; ///< Vertex in fault mesh.
+ };
+
+ // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+ /** Initialize auxiliary cohesive cell information.
+ *
+ * @param mesh Finite-element mesh of the domain.
+ */
+ void _initializeCohesiveInfo(const topology::Mesh& mesh);
+
+ // PROTECTED MEMBERS //////////////////////////////////////////////////
+protected :
+
+ /// Array of cohesive vertex information.
+ std::vector<CohesiveInfo> _cohesiveVertices;
+
// PRIVATE METHODS ////////////////////////////////////////////////////
private :
@@ -284,6 +312,9 @@
typedef std::map<std::string, EqKinSrc*> srcs_type;
+ // PRIVATE STRUCTS ////////////////////////////////////////////////////
+private :
+
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
More information about the CIG-COMMITS
mailing list