[cig-commits] r13375 - in short/3D/PyLith/trunk: libsrc/faults pylith/faults
brad at geodynamics.org
brad at geodynamics.org
Fri Nov 21 17:06:06 PST 2008
Author: brad
Date: 2008-11-21 17:06:06 -0800 (Fri, 21 Nov 2008)
New Revision: 13375
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
Log:
Change output of change in tractions from increments to total values for quasi-static solutions.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-11-22 01:05:23 UTC (rev 13374)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-11-22 01:06:06 UTC (rev 13375)
@@ -122,10 +122,21 @@
_cumSlip = new real_section_type(_faultMesh->comm(), _faultMesh->debug());
_cumSlip->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(),
vertices->end()),
- *std::max_element(vertices->begin(), vertices->end())+1));
+ *std::max_element(vertices->begin(),
+ vertices->end())+1));
_cumSlip->setFiberDimension(vertices, cs->spaceDim());
_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
// ----------------------------------------------------------------------
@@ -425,6 +436,47 @@
} // integrateJacobianAssembled
// ----------------------------------------------------------------------
+// Update state variables as needed.
+void
+pylith::faults::FaultCohesiveKin::updateState(const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh)
+{ // updateState
+ assert(0 != fields);
+ assert(!mesh.isNull());
+ assert(!_faultMesh.isNull());
+ assert(!_cumSoln.isNull());
+
+ // Update cumulateive solution for calculating total change in tractions.
+
+ 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];
+ std::cout << "Fault vertex " << fv << " getting value " << vertexSoln[0] << " from domain vertex " << *v_iter << std::endl;
+ _cumSoln->updateAddPoint(fv, vertexSoln);
+ } // if
+
+#if 1
+ _faultMesh->view("FAULT MESH");
+ _cumSoln->view("CUMULATIVE SOLUTION");
+#endif
+} // updateState
+
+// ----------------------------------------------------------------------
// Verify configuration is acceptable.
void
pylith::faults::FaultCohesiveKin::verifyConfiguration(const ALE::Obj<Mesh>& mesh) const
@@ -528,8 +580,7 @@
return _bufferTmp;
} else if (0 == strcasecmp("traction_change", name)) {
*fieldType = VECTOR_FIELD;
- const ALE::Obj<real_section_type>& solution = fields->getSolution();
- _calcTractionsChange(&_bufferVertexVector, mesh, solution);
+ _calcTractionsChange(&_bufferVertexVector);
return _bufferVertexVector;
} else {
@@ -755,7 +806,10 @@
_pseudoStiffness = new real_section_type(_faultMesh->comm(),
_faultMesh->debug());
assert(!_pseudoStiffness.isNull());
- _pseudoStiffness->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(), vertices->end()), *std::max_element(vertices->begin(), vertices->end())+1));
+ _pseudoStiffness->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(),
+ vertices->end()),
+ *std::max_element(vertices->begin(),
+ vertices->end())+1));
_pseudoStiffness->setFiberDimension(vertices, 1);
_faultMesh->allocate(_pseudoStiffness);
@@ -814,7 +868,10 @@
_area = new real_section_type(_faultMesh->comm(),
_faultMesh->debug());
assert(!_area.isNull());
- _area->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(), vertices->end()), *std::max_element(vertices->begin(), vertices->end())+1));
+ _area->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(),
+ vertices->end()),
+ *std::max_element(vertices->begin(),
+ vertices->end())+1));
_area->setFiberDimension(vertices, 1);
_faultMesh->allocate(_area);
@@ -890,27 +947,20 @@
// NOTE: We must convert vertex labels to fault vertex labels
void
pylith::faults::FaultCohesiveKin::_calcTractionsChange(
- ALE::Obj<real_section_type>* tractions,
- const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<real_section_type>& solution)
+ ALE::Obj<real_section_type>* tractions)
{ // _calcTractionsChange
assert(0 != tractions);
- assert(!solution.isNull());
assert(!_faultMesh.isNull());
+ assert(!_cumSoln.isNull());
assert(!_pseudoStiffness.isNull());
assert(!_area.isNull());
- const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
- assert(!vertices.isNull());
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- const int numVertices = vertices->size();
- Mesh::renumbering_type& renumbering = _faultMesh->getRenumbering();
- Mesh::point_type firstFaultVertex = -1;
+ const ALE::Obj<Mesh::label_sequence>& vertices =
+ _faultMesh->depthStratum(0);
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+ const int numVertices = vertices->size();
- const int fiberDim = solution->getFiberDimension(*vertices->begin());
- double_array tractionValues(fiberDim);
-
-#if 0
+#if 0 // MOVE TO SEPARATE CHECK METHOD
// Check fault mesh and volume mesh coordinates
const ALE::Obj<real_section_type>& coordinates = mesh->getRealSection("coordinates");
const ALE::Obj<real_section_type>& fCoordinates = _faultMesh->getRealSection("coordinates");
@@ -932,63 +982,53 @@
}
#endif
+ // Fiber dimension of tractions matches fiber dimension of solution
+ const int fiberDim = _cumSoln->getFiberDimension(*vertices->begin());
+ double_array tractionValues(fiberDim);
+
// Allocate buffer for tractions field (if nec.).
- for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != verticesEnd; ++v_iter) {
- if (renumbering.find(*v_iter) != renumbering.end()) {
- firstFaultVertex = renumbering[*v_iter];
- }
- }
if (tractions->isNull() ||
- fiberDim != (*tractions)->getFiberDimension(firstFaultVertex)) {
+ fiberDim != (*tractions)->getFiberDimension(*vertices->begin())) {
*tractions = new real_section_type(_faultMesh->comm(), _faultMesh->debug());
- int minE = _faultMesh->getSieve()->getChart().min();
- int maxE = _faultMesh->getSieve()->getChart().max();
-
- for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != verticesEnd; ++v_iter) {
- if (renumbering.find(*v_iter) != renumbering.end()) {
- minE = std::min(minE, renumbering[*v_iter]);
- maxE = std::max(maxE, renumbering[*v_iter]);
- }
- }
- (*tractions)->setChart(real_section_type::chart_type(minE, maxE+1));
- for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != verticesEnd; ++v_iter) {
- if (renumbering.find(*v_iter) != renumbering.end()) {
- (*tractions)->setFiberDimension(renumbering[*v_iter], fiberDim);
- }
- }
- (*tractions)->allocatePoint();
+ (*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);
+ _faultMesh->allocate(*tractions);
+ assert(!tractions->isNull());
} // if
- for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != verticesEnd; ++v_iter) {
- if (renumbering.find(*v_iter) == renumbering.end()) continue;
- const int fv = renumbering[*v_iter];
- assert(fiberDim == solution->getFiberDimension(*v_iter));
- assert(fiberDim == (*tractions)->getFiberDimension(fv));
- assert(1 == _pseudoStiffness->getFiberDimension(fv));
- assert(1 == _area->getFiberDimension(fv));
+ 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));
const real_section_type::value_type* solutionValues =
- solution->restrictPoint(*v_iter);
+ _cumSoln->restrictPoint(*v_iter);
assert(0 != solutionValues);
const real_section_type::value_type* pseudoStiffValue =
- _pseudoStiffness->restrictPoint(fv);
+ _pseudoStiffness->restrictPoint(*v_iter);
assert(0 != _pseudoStiffness);
const real_section_type::value_type* areaValue =
- _area->restrictPoint(fv);
+ _area->restrictPoint(*v_iter);
assert(0 != _area);
const double scale = pseudoStiffValue[0] / areaValue[0];
for (int i=0; i < fiberDim; ++i)
tractionValues[i] = solutionValues[i] * scale;
- (*tractions)->updatePoint(fv, &tractionValues[0]);
+ (*tractions)->updatePoint(*v_iter, &tractionValues[0]);
} // for
PetscLogFlops(numVertices * (1 + fiberDim) );
-#if 0
+#if 1
_faultMesh->view("FAULT MESH");
- solution->view("SOLUTION");
+ _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-22 01:05:23 UTC (rev 13374)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-11-22 01:06:06 UTC (rev 13375)
@@ -137,6 +137,16 @@
topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh);
+ /** Update state variables as needed.
+ *
+ * @param t Current time
+ * @param fields Solution fields
+ * @param mesh Finite-element mesh
+ */
+ void updateState(const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh);
+
/** Verify configuration is acceptable.
*
* @param mesh Finite-element mesh
@@ -211,11 +221,8 @@
/** Compute change in tractions on fault surface using solution.
*
* @param tractions Field for tractions.
- * @param solution Solution field.
*/
- void _calcTractionsChange(ALE::Obj<real_section_type>* tractions,
- const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<real_section_type>& solution);
+ void _calcTractionsChange(ALE::Obj<real_section_type>* tractions);
/// Allocate scalar field for output of vertex information.
void _allocateBufferVertexScalar(void);
@@ -261,6 +268,10 @@
/// 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.
Modified: short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py 2008-11-22 01:05:23 UTC (rev 13374)
+++ short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py 2008-11-22 01:06:06 UTC (rev 13375)
@@ -161,6 +161,7 @@
logEvent = "%spoststep" % self._loggingPrefix
self._logger.eventBegin(logEvent)
+ Integrator.poststep(self, t, dt, totalTime, fields)
FaultCohesive.poststep(self, t, dt, totalTime, fields)
self._logger.eventEnd(logEvent)
More information about the CIG-COMMITS
mailing list