[cig-commits] r13395 - short/3D/PyLith/trunk/libsrc/faults
brad at geodynamics.org
brad at geodynamics.org
Tue Nov 25 10:26:05 PST 2008
Author: brad
Date: 2008-11-25 10:26:05 -0800 (Tue, 25 Nov 2008)
New Revision: 13395
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
Log:
Fixed integrating residual for fault in parallel. Removed obsolete members from fault.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-11-25 18:16:44 UTC (rev 13394)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-11-25 18:26:05 UTC (rev 13395)
@@ -128,28 +128,19 @@
_faultMesh->allocate(_cumSlip);
assert(!_cumSlip.isNull());
- // Allocate cumulative solution field
- _cumSoln = new real_section_type(_faultMesh->comm(), _faultMesh->debug());
- _cumSoln->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(),
- vertices->end()),
- *std::max_element(vertices->begin(),
- vertices->end())+1));
- _cumSoln->setFiberDimension(vertices, cs->spaceDim());
- _faultMesh->allocate(_cumSoln);
- assert(!_cumSoln.isNull());
} // initialize
// ----------------------------------------------------------------------
-// Integrate contribution of cohesive cells to residual term that do
-// not require assembly across cells, vertices, or processors.
+// Integrate contribution of cohesive cells to residual term that
+// require assembly across processors.
void
-pylith::faults::FaultCohesiveKin::integrateResidualAssembled(
+pylith::faults::FaultCohesiveKin::integrateResidual(
const ALE::Obj<real_section_type>& residual,
const double t,
topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh,
const spatialdata::geocoords::CoordSys* cs)
-{ // integrateResidualAssembled
+{ // integrateResidual
assert(!residual.isNull());
assert(0 != fields);
assert(!mesh.isNull());
@@ -161,16 +152,27 @@
// slip
// * DOF k: slip values
+ if (!_useSolnIncr)
+ return;
+
// Get cell information and setup storage for cell data
const int spaceDim = _quadrature->spaceDim();
const int orientationSize = spaceDim*spaceDim;
- const int numConstraintVert = _quadrature->numBasis();
+ const int numBasis = _quadrature->numBasis();
+ const int numConstraintVert = numBasis;
const int numCorners = 3*numConstraintVert; // cohesive cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
double_array cellOrientation(numConstraintVert*orientationSize);
- double_array cellResidual(numCorners*spaceDim);
double_array cellSoln(numCorners*spaceDim);
- double_array cellSlip(numConstraintVert*spaceDim);
double_array cellStiffness(numConstraintVert);
+ double_array cellResidual(numCorners*spaceDim);
+ double_array cellArea(numConstraintVert);
+ double_array cellAreaAssembled(numConstraintVert);
// Get cohesive cells
const ALE::Obj<Mesh::label_sequence>& cellsCohesive =
@@ -186,42 +188,22 @@
const ALE::Obj<real_section_type>& solution = fields->getSolution();
assert(!solution.isNull());
- assert(!_slip.isNull());
- _slip->zero();
- if (!_useSolnIncr) {
- // 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) {
- EqKinSrc* src = s_iter->second;
- assert(0 != src);
- if (t >= src->originTime())
- src->slip(_slip, t, _faultMesh);
- } // for
- } else {
- // Compute increment of 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) {
- EqKinSrc* src = s_iter->second;
- assert(0 != src);
- if (t >= src->originTime())
- src->slipIncr(_slip, t-_dt, t, _faultMesh);
- } // for
- } // else
-
for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
c_iter != cellsCohesiveEnd;
++c_iter) {
const Mesh::point_type c_fault = _cohesiveToFault[*c_iter];
+ cellArea = 0.0;
+ cellResidual = 0.0;
- // Get residual at fault cell's vertices. We want to only
- // overwrite values that need updating.
- mesh->restrictClosure(residual, *c_iter, &cellResidual[0],
- cellResidual.size());
-
+ // Compute contributory area for cell (to weight contributions)
+ 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];
+ cellArea[iBasis] += dArea;
+ } // for
+ } // for
+
// Get orientations at fault cell's vertices.
_faultMesh->restrictClosure(_orientation, c_fault, &cellOrientation[0],
cellOrientation.size());
@@ -230,9 +212,10 @@
_faultMesh->restrictClosure(_pseudoStiffness, c_fault, &cellStiffness[0],
cellStiffness.size());
- // Get slip at fault cell's vertices.
- _faultMesh->restrictClosure(_slip, c_fault, &cellSlip[0], cellSlip.size());
-
+ // Get pseudo stiffness at fault cell's vertices.
+ _faultMesh->restrictClosure(_area, c_fault, &cellAreaAssembled[0],
+ cellAreaAssembled.size());
+
// Get solution at cohesive cell's vertices.
mesh->restrictClosure(solution, *c_iter, &cellSoln[0], cellSoln.size());
@@ -244,40 +227,31 @@
const int indexK = iConstraint + 2*numConstraintVert;
const double pseudoStiffness = cellStiffness[iConstraint];
-
- // Set slip values in residual vector; contributions are at DOF of
- // constraint vertices (k) of the cohesive cells
- for (int iDim=0; iDim < spaceDim; ++iDim)
- cellResidual[indexK*spaceDim+iDim] =
- cellSlip[iConstraint*spaceDim+iDim];
+ const double wt = pseudoStiffness *
+ cellArea[iConstraint] / cellAreaAssembled[iConstraint];
- if (_useSolnIncr) {
- // Get orientation at constraint vertex
- const real_section_type::value_type* constraintOrient =
- &cellOrientation[iConstraint*orientationSize];
- assert(0 != constraintOrient);
-
- // Entries associated with constraint forces applied at node i
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- cellResidual[indexI*spaceDim+iDim] = 0.0; // ?
- for (int kDim=0; kDim < spaceDim; ++kDim)
- cellResidual[indexI*spaceDim+iDim] -=
- cellSoln[indexK*spaceDim+kDim] *
- -constraintOrient[kDim*spaceDim+iDim] * pseudoStiffness;
- } // for
-
+ // Get orientation at constraint vertex
+ const real_section_type::value_type* constraintOrient =
+ &cellOrientation[iConstraint*orientationSize];
+ assert(0 != constraintOrient);
+
+ // Entries associated with constraint forces applied at node i
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ for (int kDim=0; kDim < spaceDim; ++kDim)
+ cellResidual[indexI*spaceDim+iDim] -=
+ cellSoln[indexK*spaceDim+kDim] *
+ -constraintOrient[kDim*spaceDim+iDim] * wt;
+ } // for
+
// Entries associated with constraint forces applied at node j
- for (int jDim=0; jDim < spaceDim; ++jDim) {
- cellResidual[indexJ*spaceDim+jDim] = 0.0; // ?
- for (int kDim=0; kDim < spaceDim; ++kDim)
- cellResidual[indexJ*spaceDim+jDim] -=
- cellSoln[indexK*spaceDim+kDim] *
- constraintOrient[kDim*spaceDim+jDim] * pseudoStiffness;
- } // for
- } // if
+ for (int jDim=0; jDim < spaceDim; ++jDim) {
+ for (int kDim=0; kDim < spaceDim; ++kDim)
+ cellResidual[indexJ*spaceDim+jDim] -=
+ cellSoln[indexK*spaceDim+kDim] *
+ constraintOrient[kDim*spaceDim+jDim] * wt;
+ } // for
} // for
-
#if 0 // DEBUGGING
std::cout << "Updating fault residual for cell " << *c_iter << std::endl;
for(int i = 0; i < numConstraintVert; ++i) {
@@ -294,10 +268,79 @@
}
#endif
- // Assemble cell contribution into field
- mesh->update(residual, *c_iter, &cellResidual[0]);
+ mesh->updateAdd(residual, *c_iter, &cellResidual[0]);
} // for
+
+ // FIX THIS
PetscLogFlops(cellsCohesiveSize*numConstraintVert*spaceDim*spaceDim*7);
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Integrate contribution of cohesive cells to residual term that do
+// not require assembly across cells, vertices, or processors.
+void
+pylith::faults::FaultCohesiveKin::integrateResidualAssembled(
+ const ALE::Obj<real_section_type>& residual,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs)
+{ // integrateResidualAssembled
+ assert(!residual.isNull());
+ assert(0 != fields);
+ assert(!mesh.isNull());
+
+ // Cohesive cells with normal vertices i and j, and constraint
+ // vertex k make 2 contributions to the residual:
+ //
+ // * DOF i and j: internal forces in soln field associated with
+ // slip
+ // * DOF k: slip values
+
+ assert(!_slip.isNull());
+ _slip->zero();
+ if (!_useSolnIncr) {
+ // 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) {
+ EqKinSrc* src = s_iter->second;
+ assert(0 != src);
+ if (t >= src->originTime())
+ src->slip(_slip, t, _faultMesh);
+ } // for
+ } else {
+ // Compute increment of 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) {
+ EqKinSrc* src = s_iter->second;
+ assert(0 != src);
+ if (t >= src->originTime())
+ src->slipIncr(_slip, t-_dt, t, _faultMesh);
+ } // for
+ } // else
+
+ const int spaceDim = _quadrature->spaceDim();
+ const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+ Mesh::renumbering_type& renumbering = _faultMesh->getRenumbering();
+ for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != verticesEnd;
+ ++v_iter)
+ if (renumbering.find(*v_iter) != renumbering.end()) {
+ const int vertexFault = renumbering[*v_iter];
+ const int vertexMesh = *v_iter;
+ const real_section_type::value_type* slip =
+ _slip->restrictPoint(vertexFault);
+ assert(spaceDim == _slip->getFiberDimension(vertexFault));
+ assert(spaceDim == residual->getFiberDimension(vertexMesh));
+ assert(0 != slip);
+ residual->updatePoint(vertexMesh, slip);
+ } // if
} // integrateResidualAssembled
// ----------------------------------------------------------------------
@@ -457,33 +500,6 @@
if (!_useSolnIncr)
_cumSlip->zero();
_cumSlip->add(_cumSlip, _slip);
-
- // Update cumulative solution for calculating total change in tractions.
- assert(!_cumSoln.isNull());
- if (!_useSolnIncr)
- _cumSoln->zero();
-
- const ALE::Obj<real_section_type>& solution = fields->getSolution();
-
- const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
- assert(!vertices.isNull());
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- Mesh::renumbering_type& renumbering = _faultMesh->getRenumbering();
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != verticesEnd;
- ++v_iter)
- if (renumbering.find(*v_iter) != renumbering.end()) {
- const real_section_type::value_type* vertexSoln =
- solution->restrictPoint(*v_iter);
- assert(0 != vertexSoln);
- const int fv = renumbering[*v_iter];
- _cumSoln->updateAddPoint(fv, vertexSoln);
- } // if
-
-#if 0 // DEBUGGING
- _faultMesh->view("FAULT MESH");
- _cumSoln->view("CUMULATIVE SOLUTION");
-#endif
} // updateState
// ----------------------------------------------------------------------
@@ -590,7 +606,8 @@
return _bufferTmp;
} else if (0 == strcasecmp("traction_change", name)) {
*fieldType = VECTOR_FIELD;
- _calcTractionsChange(&_bufferVertexVector);
+ const ALE::Obj<real_section_type>& solution = fields->getSolution();
+ _calcTractionsChange(&_bufferVertexVector, mesh, solution);
return _bufferVertexVector;
} else {
@@ -711,20 +728,11 @@
ncV.clear();
} // for
-#if 0 // DEBUGGING
- _orientation->view("ORIENTATION Before complete");
- _orientation->setDebug(2);
-#endif
// Assemble orientation information
-
ALE::Completion::completeSectionAdd(_faultMesh->getSendOverlap(),
_faultMesh->getRecvOverlap(),
_orientation, _orientation);
-#if 0 // DEBUGGING
- _orientation->view("ORIENTATION After complete");
-#endif
-
// Loop over vertices, make orientation information unit magnitude
double_array vertexDir(orientationSize);
int count = 0;
@@ -903,8 +911,6 @@
const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
const double_array& quadWts = _quadrature->quadWts();
assert(quadWts.size() == numQuadPts);
- const int jacobianSize = (cellDim > 0) ? spaceDim * cellDim : 1;
- double_array jacobian(jacobianSize);
double jacobianDet = 0;
double_array cellArea(numBasis);
double_array cellVertices(numBasis*spaceDim);
@@ -920,7 +926,7 @@
const double_array& basis = _quadrature->basis();
const double_array& jacobianDet = _quadrature->jacobianDet();
- // Compute action for traction bc terms
+ // Compute area
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
for (int iBasis=0; iBasis < numBasis; ++iBasis) {
@@ -948,19 +954,27 @@
// NOTE: We must convert vertex labels to fault vertex labels
void
pylith::faults::FaultCohesiveKin::_calcTractionsChange(
- ALE::Obj<real_section_type>* tractions)
+ ALE::Obj<real_section_type>* tractions,
+ const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<real_section_type>& solution)
{ // _calcTractionsChange
assert(0 != tractions);
+ assert(!mesh.isNull());
+ assert(!solution.isNull());
assert(!_faultMesh.isNull());
- assert(!_cumSoln.isNull());
assert(!_pseudoStiffness.isNull());
assert(!_area.isNull());
const ALE::Obj<Mesh::label_sequence>& vertices =
- _faultMesh->depthStratum(0);
+ mesh->depthStratum(0);
const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- const int numVertices = vertices->size();
+ const ALE::Obj<Mesh::label_sequence>& fvertices =
+ _faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator fverticesEnd = fvertices->end();
+ const int numFaultVertices = fvertices->size();
+ Mesh::renumbering_type& renumbering = _faultMesh->getRenumbering();
+
#if 0 // MOVE TO SEPARATE CHECK METHOD
// Check fault mesh and volume mesh coordinates
const ALE::Obj<real_section_type>& coordinates = mesh->getRealSection("coordinates");
@@ -983,53 +997,55 @@
}
#endif
- // Fiber dimension of tractions matches fiber dimension of solution
- const int fiberDim = _cumSoln->getFiberDimension(*vertices->begin());
+ // Fiber dimension of tractions matches spatial dimension.
+ const int fiberDim = _quadrature->spaceDim();
double_array tractionValues(fiberDim);
// Allocate buffer for tractions field (if nec.).
if (tractions->isNull() ||
- fiberDim != (*tractions)->getFiberDimension(*vertices->begin())) {
+ fiberDim != (*tractions)->getFiberDimension(*fvertices->begin())) {
*tractions = new real_section_type(_faultMesh->comm(), _faultMesh->debug());
- (*tractions)->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(),
- vertices->end()),
- *std::max_element(vertices->begin(),
- vertices->end())+1));
- (*tractions)->setFiberDimension(vertices, fiberDim);
+ (*tractions)->setChart(real_section_type::chart_type(*std::min_element(fvertices->begin(),
+ fvertices->end()),
+ *std::max_element(fvertices->begin(),
+ fvertices->end())+1));
+ (*tractions)->setFiberDimension(fvertices, fiberDim);
_faultMesh->allocate(*tractions);
assert(!tractions->isNull());
} // if
for (Mesh::label_sequence::iterator v_iter = vertices->begin();
v_iter != verticesEnd;
- ++v_iter) {
- assert(fiberDim == _cumSoln->getFiberDimension(*v_iter));
- assert(fiberDim == (*tractions)->getFiberDimension(*v_iter));
- assert(1 == _pseudoStiffness->getFiberDimension(*v_iter));
- assert(1 == _area->getFiberDimension(*v_iter));
+ ++v_iter)
+ if (renumbering.find(*v_iter) != renumbering.end()) {
+ const int vertexMesh = *v_iter;
+ const int vertexFault = renumbering[*v_iter];
+ assert(fiberDim == solution->getFiberDimension(vertexMesh));
+ assert(fiberDim == (*tractions)->getFiberDimension(vertexFault));
+ assert(1 == _pseudoStiffness->getFiberDimension(vertexFault));
+ assert(1 == _area->getFiberDimension(vertexFault));
- const real_section_type::value_type* solutionValues =
- _cumSoln->restrictPoint(*v_iter);
- assert(0 != solutionValues);
- const real_section_type::value_type* pseudoStiffValue =
- _pseudoStiffness->restrictPoint(*v_iter);
- assert(0 != _pseudoStiffness);
- const real_section_type::value_type* areaValue =
- _area->restrictPoint(*v_iter);
- assert(0 != _area);
+ const real_section_type::value_type* solutionValues =
+ solution->restrictPoint(vertexMesh);
+ assert(0 != solutionValues);
+ const real_section_type::value_type* pseudoStiffValue =
+ _pseudoStiffness->restrictPoint(vertexFault);
+ assert(0 != _pseudoStiffness);
+ const real_section_type::value_type* areaValue =
+ _area->restrictPoint(vertexFault);
+ assert(0 != _area);
- const double scale = pseudoStiffValue[0] / areaValue[0];
- for (int i=0; i < fiberDim; ++i)
- tractionValues[i] = solutionValues[i] * scale;
+ const double scale = pseudoStiffValue[0] / areaValue[0];
+ for (int i=0; i < fiberDim; ++i)
+ tractionValues[i] = solutionValues[i] * scale;
- (*tractions)->updatePoint(*v_iter, &tractionValues[0]);
- } // for
+ (*tractions)->updatePoint(vertexFault, &tractionValues[0]);
+ } // if
- PetscLogFlops(numVertices * (1 + fiberDim) );
+ PetscLogFlops(numFaultVertices * (1 + fiberDim) );
#if 0 // DEBUGGING
_faultMesh->view("FAULT MESH");
- _cumSoln->view("CUMULATIVE SOLUTION");
_area->view("AREA");
_pseudoStiffness->view("CONDITIONING");
(*tractions)->view("TRACTIONS");
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-11-25 18:16:44 UTC (rev 13394)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-11-25 18:26:05 UTC (rev 13395)
@@ -110,6 +110,20 @@
spatialdata::spatialdb::SpatialDB* matDB);
/** Integrate contributions to residual term (r) for operator that
+ * require assembly across processors.
+ *
+ * @param residual Field containing values for residual
+ * @param t Current time
+ * @param fields Solution fields
+ * @param mesh Finite-element mesh
+ */
+ void integrateResidual(const ALE::Obj<real_section_type>& residual,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs);
+
+ /** Integrate contributions to residual term (r) for operator that
* do not require assembly across cells, vertices, or processors.
*
* @param residual Field containing values for residual
@@ -221,8 +235,12 @@
/** Compute change in tractions on fault surface using solution.
*
* @param tractions Field for tractions.
+ * @param mesh Finite-element mesh for domain
+ * @param solution Solution over domain
*/
- void _calcTractionsChange(ALE::Obj<real_section_type>* tractions);
+ void _calcTractionsChange(ALE::Obj<real_section_type>* tractions,
+ const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<real_section_type>& solution);
/// Allocate scalar field for output of vertex information.
void _allocateBufferVertexScalar(void);
@@ -268,10 +286,6 @@
/// Field over the fault mesh vertices of vector field of cumulative slip.
ALE::Obj<real_section_type> _cumSlip;
- /// Field over the fault mesh vertices of vector field of cumulative
- /// solution. This permits computing the total change in tractions.
- ALE::Obj<real_section_type> _cumSoln;
-
std::map<Mesh::point_type, Mesh::point_type> _cohesiveToFault;
/// Scalar field for vertex information over fault mesh.
More information about the CIG-COMMITS
mailing list