[cig-commits] r12976 - short/3D/PyLith/trunk/libsrc/faults
knepley at geodynamics.org
knepley at geodynamics.org
Tue Sep 30 20:38:35 PDT 2008
Author: knepley
Date: 2008-09-30 20:38:35 -0700 (Tue, 30 Sep 2008)
New Revision: 12976
Modified:
short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
Log:
Fixed 1 processor fault handling with renumbering
Modified: short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2008-09-30 23:47:09 UTC (rev 12975)
+++ short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2008-10-01 03:38:35 UTC (rev 12976)
@@ -856,7 +856,6 @@
const ALE::Obj<sieve_type> ifaultSieve = new sieve_type(sieve->comm(), sieve->debug());
ALE::Obj<ALE::Mesh> fault = new ALE::Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
ALE::Obj<ALE::Mesh::sieve_type> faultSieve = new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
- Mesh::renumbering_type& fRenumbering = (*ifault)->getRenumbering();
cohesiveToFault->clear();
const ALE::Obj<Mesh::label_sequence>& cohesiveCells = mesh->getLabelStratum("material-id", materialId);
@@ -886,7 +885,6 @@
// Use first vertices (negative side of the fault) for fault mesh
for(int i = 0; i < faceSize; ++i) {
faultSieve->addArrow(cone[i], face, color++);
- fRenumbering[cone[i]] = cone[i];
}
} else {
const int faceSize = coneSize / 3;
@@ -895,7 +893,6 @@
// Use last vertices (contraints) for fault mesh
for(int i = 2*faceSize; i < 3*faceSize; ++i) {
faultSieve->addArrow(cone[i], face, color++);
- fRenumbering[cone[i]] = cone[i];
}
} // if/else
(*cohesiveToFault)[*c_iter] = face;
@@ -906,10 +903,14 @@
fault->stratify();
// Convert fault to an IMesh
- std::map<Mesh::point_type,Mesh::point_type> renumbering;
+ Mesh::renumbering_type& fRenumbering = (*ifault)->getRenumbering();
(*ifault)->setSieve(ifaultSieve);
- ALE::ISieveConverter::convertMesh(*fault, *(*ifault), renumbering, false);
- renumbering.clear();
+ //ALE::ISieveConverter::convertMesh(*fault, *(*ifault), fRenumbering, true);
+ {
+ ALE::ISieveConverter::convertSieve(*fault->getSieve(), *(*ifault)->getSieve(), fRenumbering, true);
+ (*ifault)->stratify();
+ ALE::ISieveConverter::convertOrientation(*fault->getSieve(), *(*ifault)->getSieve(), fRenumbering, fault->getArrowSection("orientation").ptr());
+ }
fault = NULL;
faultSieve = NULL;
@@ -919,25 +920,29 @@
(*cohesiveToFault)[*c_iter] = *f_iter;
}
-#if 1
+#if 0
(*ifault)->setRealSection("coordinates", mesh->getRealSection("coordinates"));
#else
const ALE::Obj<Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
- const ALE::Obj<Mesh::real_section_type>& fCoordinates = (*fault)->getRealSection("coordinates");
- const ALE::Obj<Mesh::label_sequence>& vertices = (*fault)->depthStratum(0);
+ const ALE::Obj<Mesh::real_section_type>& fCoordinates = (*ifault)->getRealSection("coordinates");
+ const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
const Mesh::label_sequence::iterator vBegin = vertices->begin();
const Mesh::label_sequence::iterator vEnd = vertices->end();
+ fCoordinates->setChart(Mesh::real_section_type::chart_type((*ifault)->heightStratum(0)->size(),
+ (*ifault)->getSieve()->getChart().max()));
for(Mesh::label_sequence::iterator v_iter = vBegin;
v_iter != vEnd;
++v_iter) {
- fCoordinates->setFiberDimension(*v_iter, coordinates->getFiberDimension(*v_iter));
+ if (fRenumbering.find(*v_iter) == fRenumbering.end()) continue;
+ fCoordinates->setFiberDimension(fRenumbering[*v_iter], coordinates->getFiberDimension(*v_iter));
}
- (*fault)->allocate(fCoordinates);
+ fCoordinates->allocatePoint();
for(Mesh::label_sequence::iterator v_iter = vBegin;
v_iter != vEnd;
++v_iter) {
- fCoordinates->updatePoint(*v_iter, coordinates->restrictPoint(*v_iter));
+ if (fRenumbering.find(*v_iter) == fRenumbering.end()) continue;
+ fCoordinates->updatePoint(fRenumbering[*v_iter], coordinates->restrictPoint(*v_iter));
}
#endif
(*ifault)->view("Parallel fault mesh");
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-09-30 23:47:09 UTC (rev 12975)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2008-10-01 03:38:35 UTC (rev 12976)
@@ -85,6 +85,7 @@
CohesiveTopology::createParallel(&_faultMesh, &_cohesiveToFault, mesh, id(),
_useLagrangeConstraints());
+ _faultMesh->getLabel("height")->view("Fault mesh height");
//_faultMesh->view("FAULT MESH");
@@ -524,7 +525,7 @@
} else if (0 == strcasecmp("traction_change", name)) {
*fieldType = VECTOR_FIELD;
const ALE::Obj<real_section_type>& solution = fields->getSolution();
- _calcTractionsChange(&_bufferVertexVector, solution);
+ _calcTractionsChange(&_bufferVertexVector, mesh, solution);
return _bufferVertexVector;
} else {
@@ -916,9 +917,11 @@
// ----------------------------------------------------------------------
// Compute change in tractions on fault surface using solution.
+// 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)
{ // _calcTractionsChange
assert(0 != tractions);
@@ -927,53 +930,69 @@
assert(!_pseudoStiffness.isNull());
assert(!_area.isNull());
- const ALE::Obj<Mesh::label_sequence>& vertices =
- _faultMesh->depthStratum(0);
+ _pseudoStiffness->view("PSEUDOSTIFFNESS");
+ _area->view("AREA");
+
+ 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();
+ const Mesh::label_sequence::iterator verticesEnd = vertices->end();
+ const int numVertices = vertices->size();
+ Mesh::renumbering_type& renumbering = _faultMesh->getRenumbering();
const int fiberDim = solution->getFiberDimension(*vertices->begin());
double_array tractionValues(fiberDim);
// Allocate buffer for tractions field (if nec.).
if (tractions->isNull() ||
- fiberDim != (*tractions)->getFiberDimension(*vertices->begin())) {
+ fiberDim != (*tractions)->getFiberDimension(renumbering[*vertices->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);
- _faultMesh->allocate(*tractions);
+ 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();
} // if
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != verticesEnd;
- ++v_iter) {
+ 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(*v_iter));
- assert(1 == _pseudoStiffness->getFiberDimension(*v_iter));
- assert(1 == _area->getFiberDimension(*v_iter));
+ assert(fiberDim == (*tractions)->getFiberDimension(fv));
+ assert(1 == _pseudoStiffness->getFiberDimension(fv));
+ assert(1 == _area->getFiberDimension(fv));
const real_section_type::value_type* solutionValues =
solution->restrictPoint(*v_iter);
assert(0 != solutionValues);
const real_section_type::value_type* pseudoStiffValue =
- _pseudoStiffness->restrictPoint(*v_iter);
+ _pseudoStiffness->restrictPoint(fv);
assert(0 != _pseudoStiffness);
const real_section_type::value_type* areaValue =
- _area->restrictPoint(*v_iter);
+ _area->restrictPoint(fv);
assert(0 != _area);
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]);
+ (*tractions)->updatePoint(fv, &tractionValues[0]);
} // for
PetscLogFlops(numVertices * (1 + fiberDim) );
//solution->view("SOLUTION");
- //(*tractions)->view("TRACTIONS");
+ (*tractions)->view("TRACTIONS");
} // _calcTractionsChange
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-09-30 23:47:09 UTC (rev 12975)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-10-01 03:38:35 UTC (rev 12976)
@@ -217,6 +217,7 @@
* @param solution Solution field.
*/
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.
More information about the cig-commits
mailing list