[cig-commits] r14601 - in short/3D/PyLith/branches/pylith-swig: . libsrc libsrc/faults libsrc/meshio libsrc/topology unittests/libtests unittests/libtests/faults
brad at geodynamics.org
brad at geodynamics.org
Sun Apr 5 20:55:51 PDT 2009
Author: brad
Date: 2009-04-05 20:55:49 -0700 (Sun, 05 Apr 2009)
New Revision: 14601
Modified:
short/3D/PyLith/branches/pylith-swig/TODO
short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.hh
short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.cc
short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.hh
short/3D/PyLith/branches/pylith-swig/libsrc/meshio/UCDFaultFile.cc
short/3D/PyLith/branches/pylith-swig/libsrc/meshio/UCDFaultFile.hh
short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.cc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/Makefile.am
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/Makefile.am
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestEqKinSrc.cc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestEqKinSrc.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFault.cc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFault.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesive.cc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesive.hh
Log:
More updating fault stuff and some unit tests.
Modified: short/3D/PyLith/branches/pylith-swig/TODO
===================================================================
--- short/3D/PyLith/branches/pylith-swig/TODO 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/TODO 2009-04-06 03:55:49 UTC (rev 14601)
@@ -25,6 +25,7 @@
(3) SolverNonlinear (get it to run then turn it over to Matt)
(4) Faults
+ Use visitors in FaultCohesiveKin
(5) Add missing unit tests
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-04-06 03:55:49 UTC (rev 14601)
@@ -36,8 +36,9 @@
faults/LiuCosSlipFn.cc \
faults/EqKinSrc.cc \
faults/FaultCohesive.cc \
- faults/FaultCohesiveKin.cc \
faults/FaultCohesiveDyn.cc \
+ faults/TopologyOps.cc \
+ faults/CohesiveTopology.cc \
feassemble/CellGeometry.cc \
feassemble/Constraint.cc \
feassemble/GeometryPoint1D.cc \
@@ -86,6 +87,7 @@
meshio/PsetFileAscii.cc \
meshio/PsetFileBinary.cc \
meshio/OutputSolnSubset.cc \
+ meshio/UCDFaultFile.cc \
problems/Formulation.cc \
problems/Solver.cc \
problems/SolverLinear.cc \
@@ -99,10 +101,8 @@
utils/EventLogger.cc \
utils/TestArray.cc
-# faults/TopologyOps.cc \
-# faults/CohesiveTopology.cc \
+# faults/FaultCohesiveKin.cc \
# materials/GenMaxwellIsotropic3D.cc \
-# meshio/UCDFaultFile.cc \
# topology/Distributor.cc \
# topology/MeshRefiner.cc \
# topology/RefineUniform.cc
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/BruneSlipFn.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -36,6 +36,7 @@
pylith::faults::BruneSlipFn::BruneSlipFn(void) :
_slipTimeVertex(0),
_riseTimeVertex(0),
+ _parameters(0),
_dbFinalSlip(0),
_dbSlipTime(0),
_dbRiseTime(0)
@@ -138,7 +139,7 @@
_dbRiseTime->open();
const char* riseTimeValues[] = {"rise-time"};
- _dbSlipTime->queryVals(slipTimeValues, 1);
+ _dbRiseTime->queryVals(riseTimeValues, 1);
// Get coordinates of vertices
const ALE::Obj<RealSection>& coordinates =
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -13,626 +13,324 @@
#include <portinfo>
#include "CohesiveTopology.hh" // implementation of object methods
-#include <Selection.hh> // Algorithms for submeshes
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+#include "TopologyOps.hh" // USES TopologyOps
+#include "TopologyVisitors.hh" // USES TopologyVisitors
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include <Selection.hh> // Algorithms for submeshes
+
#include <cassert> // USES assert()
// ----------------------------------------------------------------------
-void
-pylith::faults::CohesiveTopology::createFaultSieveFromVertices(const int dim,
- const int firstCell,
- const PointSet& faultVertices,
- const Obj<Mesh>& mesh,
- const Obj<ALE::Mesh::arrow_section_type>& orientation,
- const Obj<ALE::Mesh::sieve_type>& faultSieve,
- const bool flipFault)
-{
- typedef ALE::Selection<ALE::Mesh> selection;
- const Obj<sieve_type>& sieve = mesh->getSieve();
- const PointSet::const_iterator fvBegin = faultVertices.begin();
- const PointSet::const_iterator fvEnd = faultVertices.end();
- int curCell = firstCell;
- int curVertex = 0;
- int newElement = curCell + dim*faultVertices.size();
- int o = 1;
- ALE::Mesh::point_type f = firstCell;
- const int debug = mesh->debug();
- Obj<PointSet> face = new PointSet();
- int numCorners = 0; // The number of vertices in a mesh cell
- int faceSize = 0; // The number of vertices in a mesh face
- int *indices = NULL; // The indices of a face vertex set in a cell
- std::map<int,int*> curElement;
- std::map<int,PointArray> bdVertices;
- std::map<int,PointArray> faultFaces;
- std::map<int,oPointArray> oFaultFaces;
- PointSet faultCells;
- PointArray origVertices;
- PointArray faceVertices;
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::SieveSubMesh SieveSubMesh;
+typedef pylith::topology::Mesh::IntSection IntSection;
- if (!faultSieve->commRank()) {
- numCorners = mesh->getNumCellCorners();
- faceSize = selection::numFaceVertices(mesh);
- indices = new int[faceSize];
- }
-
- curElement[0] = &curVertex;
- curElement[dim] = &curCell;
- for(int d = 1; d < dim; d++) {
- curElement[d] = &newElement;
- }
-
- // This only works for uninterpolated meshes
- assert((mesh->depth() == 1) || (mesh->depth() == -1));
- ALE::ISieveVisitor::PointRetriever<sieve_type> sV(std::max(1, sieve->getMaxSupportSize()));
- ALE::ISieveVisitor::PointRetriever<sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
- for(PointSet::const_iterator fv_iter = fvBegin; fv_iter != fvEnd; ++fv_iter) {
- sieve->support(*fv_iter, sV);
- const Mesh::point_type *support = sV.getPoints();
-
- if (debug) std::cout << "Checking fault vertex " << *fv_iter << std::endl;
- const int sVsize = sV.getSize();
- for (int i=0; i < sVsize; ++i) {
- const int s = (!flipFault) ? i : sVsize - i - 1;
- sieve->cone(support[s], cV);
- const Mesh::point_type *cone = cV.getPoints();
-
- if (debug) std::cout << " Checking cell " << support[s] << std::endl;
- if (faultCells.find(support[s]) != faultCells.end()) {
- cV.clear();
- continue;
- }
- face->clear();
- for(int c = 0; c < cV.getSize(); ++c) {
- if (faultVertices.find(cone[c]) != fvEnd) {
- if (debug) std::cout << " contains fault vertex " << cone[c] << std::endl;
- face->insert(face->end(), cone[c]);
- } // if
- } // for
- if (face->size() > faceSize)
- throw ALE::Exception("Invalid fault mesh: Too many vertices of an "
- "element on the fault");
- if (face->size() == faceSize) {
- if (debug) std::cout << " Contains a face on the fault" << std::endl;
- ALE::Obj<sieve_type::supportSet> preFace;
- if (dim < 2) {
- preFace = faultSieve->nJoin1(face);
- } else {
- preFace = faultSieve->nJoin(face, dim);
- }
-
- if (preFace->size() > 1) {
- throw ALE::Exception("Invalid fault sieve: Multiple faces from vertex set");
- } else if (preFace->size() == 1) {
- // Add the other cell neighbor for this face
- if (dim == 0) {
- faultSieve->addArrow(*faceVertices.begin(), support[s]);
- } else {
- faultSieve->addArrow(*preFace->begin(), support[s]);
- }
- } else if (preFace->size() == 0) {
- if (debug) std::cout << " Orienting face " << f << std::endl;
- selection::getOrientedFace(mesh, support[s], face, numCorners, indices, &origVertices, &faceVertices);
- bdVertices[dim].clear();
- for(PointArray::const_iterator v_iter = faceVertices.begin(); v_iter != faceVertices.end(); ++v_iter) {
- bdVertices[dim].push_back(*v_iter);
- if (debug) std::cout << " Boundary vertex " << *v_iter << std::endl;
- }
- if (dim == 0) {
- f = *faceVertices.begin();
- }
- if (faceSize != dim+1) {
- if (debug) std::cout << " Adding hex face " << f << std::endl;
- ALE::SieveBuilder<ALE::Mesh>::buildHexFaces(faultSieve, orientation, dim, curElement, bdVertices, oFaultFaces, f, o);
- } else {
- if (debug) std::cout << " Adding simplicial face " << f << std::endl;
- ALE::SieveBuilder<ALE::Mesh>::buildFaces(faultSieve, orientation, dim, curElement, bdVertices, oFaultFaces, f, o);
- }
- faultSieve->addArrow(f, support[s]);
- //faultSieve->view("");
- f++;
- } // if/else
- faultCells.insert(support[s]);
- } // if
- cV.clear();
- } // for
- sV.clear();
- } // for
- if (!faultSieve->commRank()) delete [] indices;
-}
-
// ----------------------------------------------------------------------
void
-pylith::faults::CohesiveTopology::createFaultSieveFromFaces(const int dim,
- const int firstCell,
- const int numFaces,
- const int faultVertices[],
- const int faultCells[],
- const Obj<Mesh>& mesh,
- const Obj<ALE::Mesh::arrow_section_type>& orientation,
- const Obj<ALE::Mesh::sieve_type>& faultSieve)
-{
- typedef ALE::Selection<ALE::Mesh> selection;
- int faceSize = 0; // The number of vertices in a mesh face
- int curCell = firstCell;
- int curVertex = 0;
- int newElement = curCell + dim*numFaces;
- int o = 1;
- int f = firstCell;
- const int debug = mesh->debug();
- std::map<int,int*> curElement;
- std::map<int,PointArray> bdVertices;
- std::map<int,oPointArray> oFaultFaces;
+pylith::faults::CohesiveTopology::createFault(topology::SubMesh* faultMesh,
+ ALE::Obj<ALE::Mesh>& faultBoundary,
+ const topology::Mesh& mesh,
+ const ALE::Obj<topology::Mesh::IntSection>& groupField,
+ const bool flipFault)
+{ // createFault
+ assert(0 != faultMesh);
+ assert(!groupField.isNull());
- if (!faultSieve->commRank()) {
- faceSize = selection::numFaceVertices(mesh);
- }
+ faultMesh->coordsys(mesh);
- curElement[0] = &curVertex;
- curElement[dim] = &curCell;
- for(int d = 1; d < dim; d++) {
- curElement[d] = &newElement;
- }
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
+ faultSieveMesh =
+ new SieveSubMesh(mesh.comm(), mesh.dimension()-1, mesh.debug());
- // Loop over fault faces
- for(int face = 0; face < numFaces; ++face) {
- // Push oriented vertices of face
- bdVertices[dim].clear();
- for(int i = 0; i < faceSize; ++i) {
- bdVertices[dim].push_back(faultVertices[face*faceSize+i]);
- if (debug) std::cout << " Boundary vertex " << faultVertices[face*faceSize+i] << std::endl;
- }
- // Create face
- if (faceSize != dim+1) {
- if (debug) std::cout << " Adding hex face " << f << std::endl;
- ALE::SieveBuilder<ALE::Mesh>::buildHexFaces(faultSieve, orientation, dim, curElement, bdVertices, oFaultFaces, f, o);
- } else {
- if (debug) std::cout << " Adding simplicial face " << f << std::endl;
- ALE::SieveBuilder<ALE::Mesh>::buildFaces(faultSieve, orientation, dim, curElement, bdVertices, oFaultFaces, f, o);
- }
- // Add arrow to cells
- faultSieve->addArrow(face, faultCells[face*2+0]);
- faultSieve->addArrow(face, faultCells[face*2+1]);
- }
-}
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ assert(!sieve.isNull());
+ const ALE::Obj<SieveSubMesh::sieve_type> ifaultSieve =
+ new SieveMesh::sieve_type(sieve->comm(), sieve->debug());
+ assert(!ifaultSieve.isNull());
+ ALE::Obj<ALE::Mesh> fault =
+ new ALE::Mesh(mesh.comm(), mesh.dimension()-1, mesh.debug());
+ assert(!fault.isNull());
+ ALE::Obj<ALE::Mesh::sieve_type> faultSieve =
+ new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
+ assert(!faultSieve.isNull());
+ const int debug = mesh.debug();
-// ----------------------------------------------------------------------
-void
-pylith::faults::CohesiveTopology::orientFaultSieve(const int dim,
- const Obj<Mesh>& mesh,
- const Obj<ALE::Mesh::arrow_section_type>& orientation,
- const Obj<ALE::Mesh>& fault)
-{
- // Must check the orientation here
- typedef ALE::Selection<ALE::Mesh> selection;
- const Obj<ALE::Mesh::sieve_type>& faultSieve = fault->getSieve();
- const Mesh::point_type firstFaultCell = *fault->heightStratum(1)->begin();
- const Obj<ALE::Mesh::label_sequence>& fFaces = fault->heightStratum(2);
- const int numFaultFaces = fFaces->size();
- const int faultDepth = fault->depth()-1; // Depth of fault cells
- int numFaultCorners = 0; // The number of vertices in a fault cell
- int faultFaceSize = 0; // The number of vertices in a face between fault cells
- int faceSize = 0; // The number of vertices in a mesh face
- const int debug = fault->debug();
- Obj<PointSet> newCells = new PointSet();
- Obj<PointSet> loopCells = new PointSet();
- PointSet flippedCells; // Incorrectly oriented fault cells
- PointSet facesSeen; // Fault faces already considered
- PointSet cellsSeen; // Fault cells already matched
- PointArray faceVertices;
-
- if (!fault->commRank()) {
- faceSize = selection::numFaceVertices(mesh);
- numFaultCorners = faultSieve->nCone(firstFaultCell, faultDepth)->size();
- if (debug) std::cout << " Fault corners " << numFaultCorners << std::endl;
- if (dim == 0) {
- assert(numFaultCorners == faceSize-1);
- } else {
- assert(numFaultCorners == faceSize);
- }
- if (faultDepth == 1) {
- faultFaceSize = 1;
- } else {
- faultFaceSize = faultSieve->nCone(*fFaces->begin(), faultDepth-1)->size();
- }
- }
- if (debug) std::cout << " Fault face size " << faultFaceSize << std::endl;
-
- newCells->insert(firstFaultCell);
- while(facesSeen.size() != numFaultFaces) {
- Obj<PointSet> tmp = newCells; newCells = loopCells; loopCells = tmp;
-
- newCells->clear();
- if (!loopCells->size()) {throw ALE::Exception("Fault surface not a single connected component.");}
- // Loop over new cells
- for(PointSet::iterator c_iter = loopCells->begin(); c_iter != loopCells->end(); ++c_iter) {
- // Loop over edges of this cell
- const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& cone = faultSieve->cone(*c_iter);
- const ALE::Mesh::sieve_type::traits::coneSequence::iterator eBegin = cone->begin();
- const ALE::Mesh::sieve_type::traits::coneSequence::iterator eEnd = cone->end();
-
- for(ALE::Mesh::sieve_type::traits::coneSequence::iterator e_iter = eBegin; e_iter != eEnd; ++e_iter) {
- if (facesSeen.find(*e_iter) != facesSeen.end()) continue;
- facesSeen.insert(*e_iter);
- if (debug) std::cout << " Checking orientation of fault face " << *e_iter << std::endl;
- const Obj<ALE::Mesh::sieve_type::traits::supportSequence>& support = faultSieve->support(*e_iter);
- ALE::Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin();
-
- // Throw out boundary fault faces
- if (support->size() < 2) continue;
- ALE::Mesh::point_type cellA = *s_iter; ++s_iter;
- ALE::Mesh::point_type cellB = *s_iter;
- bool flippedA = (flippedCells.find(cellA) != flippedCells.end());
- bool flippedB = (flippedCells.find(cellB) != flippedCells.end());
- bool seenA = (cellsSeen.find(cellA) != cellsSeen.end());
- bool seenB = (cellsSeen.find(cellB) != cellsSeen.end());
-
- if (!seenA) newCells->insert(cellA);
- if (!seenB) newCells->insert(cellB);
- if (debug) std::cout << " neighboring cells " << cellA << " and " << cellB << std::endl;
- // In 1D, just check that vertices match
- if (dim == 1) {
- const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneA = faultSieve->cone(cellA);
- ALE::Mesh::sieve_type::traits::coneSequence::iterator iterA = coneA->begin();
- const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneB = faultSieve->cone(cellB);
- ALE::Mesh::sieve_type::traits::coneSequence::iterator iterB = coneB->begin();
- int posA, posB;
-
- for(posA = 0; posA < 2; ++posA, ++iterA) if (*iterA == *e_iter) break;
- for(posB = 0; posB < 2; ++posB, ++iterB) if (*iterB == *e_iter) break;
- if (debug) std::cout << " with face positions " << posA << " and " << posB << std::endl;
- if ((posA == 2) || (posB == 2)) {throw ALE::Exception("Could not find fault face in cone");}
- if ((posA == posB) ^ (flippedA || flippedB)) {
- if (debug) {
- std::cout << "Invalid orientation in fault mesh" << std::endl;
- std::cout << " fault face: " << *e_iter << " cellA: " << cellA << " cellB: " << cellB << std::endl;
- }
- if (flippedA && flippedB) {throw ALE::Exception("Attempt to flip already flipped cell: Fault mesh is non-orientable");}
- if (seenA && seenB) {throw ALE::Exception("Previous cells do not match: Fault mesh is non-orientable");}
- if (!seenA && !flippedA) {
- flippedCells.insert(cellA);
- } else if (!seenB && !flippedB) {
- flippedCells.insert(cellB);
- } else {
- throw ALE::Exception("Inconsistent mesh orientation: Fault mesh is non-orientable");
- }
- }
- } else if (dim == 2) {
- // Check orientation
- ALE::MinimalArrow<ALE::Mesh::sieve_type::point_type,ALE::Mesh::sieve_type::point_type> arrowA(*e_iter, cellA);
- const int oA = orientation->restrictPoint(arrowA)[0];
- ALE::MinimalArrow<ALE::Mesh::sieve_type::point_type,ALE::Mesh::sieve_type::point_type> arrowB(*e_iter, cellB);
- const int oB = orientation->restrictPoint(arrowB)[0];
- const bool mismatch = (oA == oB);
-
- // Truth Table
- // mismatch flips action mismatch flipA ^ flipB action
- // F 0 flips no F F F
- // F 1 flip yes F T T
- // F 2 flips no T F T
- // T 0 flips yes T T F
- // T 1 flip no
- // T 2 flips yes
- if (mismatch ^ (flippedA ^ flippedB)) {
- if (debug) {
- std::cout << "Invalid orientation in fault mesh" << std::endl;
- std::cout << " fault face: " << *e_iter << " cellA: " << cellA << " cellB: " << cellB << std::endl;
- }
- if (flippedA && flippedB) {throw ALE::Exception("Attempt to flip already flipped cell: Fault mesh is non-orientable");}
- if (seenA && seenB) {throw ALE::Exception("Previous cells do not match: Fault mesh is non-orientable");}
- if (!seenA && !flippedA) {
- flippedCells.insert(cellA);
- if (debug) {std::cout << " Scheduling cell " << cellA << " for flipping" << std::endl;}
- } else if (!seenB && !flippedB) {
- flippedCells.insert(cellB);
- if (debug) {std::cout << " Scheduling cell " << cellB << " for flipping" << std::endl;}
- } else {
- throw ALE::Exception("Inconsistent mesh orientation: Fault mesh is non-orientable");
- }
- }
- }
- cellsSeen.insert(cellA);
- cellsSeen.insert(cellB);
- }
- }
- }
- for(PointSet::const_iterator f_iter = flippedCells.begin(); f_iter != flippedCells.end(); ++f_iter) {
- if (debug) std::cout << " Reversing fault face " << *f_iter << std::endl;
- faceVertices.clear();
- const ALE::Obj<ALE::Mesh::sieve_type::traits::coneSequence>& cone = faultSieve->cone(*f_iter);
- for(ALE::Mesh::sieve_type::traits::coneSequence::iterator v_iter = cone->begin(); v_iter != cone->end(); ++v_iter) {
- faceVertices.insert(faceVertices.begin(), *v_iter);
- }
- faultSieve->clearCone(*f_iter);
- int color = 0;
- for(PointArray::const_iterator v_iter = faceVertices.begin(); v_iter != faceVertices.end(); ++v_iter) {
- faultSieve->addArrow(*v_iter, *f_iter, color++);
- }
-
- if (dim > 1) {
- // Here, they are edges, not vertices
- for(PointArray::const_iterator e_iter = faceVertices.begin(); e_iter != faceVertices.end(); ++e_iter) {
- ALE::MinimalArrow<ALE::Mesh::sieve_type::point_type,ALE::Mesh::sieve_type::point_type> arrow(*e_iter, *f_iter);
- int o = orientation->restrictPoint(arrow)[0];
-
- if (debug) std::cout << " Reversing orientation of " << *e_iter <<"-->"<<*f_iter << " from " << o << " to " << -(o+1) << std::endl;
- o = -(o+1);
- orientation->updatePoint(arrow, &o);
- }
- }
- }
- flippedCells.clear();
- for(ALE::Mesh::label_sequence::iterator e_iter = fFaces->begin(); e_iter != fFaces->end(); ++e_iter) {
- if (debug) std::cout << " Checking orientation of fault face " << *e_iter << std::endl;
- // for each face get the support (2 fault cells)
- const Obj<ALE::Mesh::sieve_type::traits::supportSequence>& support = faultSieve->support(*e_iter);
- ALE::Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin();
-
- // Throw out boundary fault faces
- if (support->size() > 1) {
- ALE::Mesh::point_type cellA = *s_iter; ++s_iter;
- ALE::Mesh::point_type cellB = *s_iter;
-
- if (debug) std::cout << " neighboring cells " << cellA << " and " << cellB << std::endl;
- // In 1D, just check that vertices match
- if (dim == 1) {
- const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneA = faultSieve->cone(cellA);
- ALE::Mesh::sieve_type::traits::coneSequence::iterator iterA = coneA->begin();
- const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneB = faultSieve->cone(cellB);
- ALE::Mesh::sieve_type::traits::coneSequence::iterator iterB = coneB->begin();
- int posA, posB;
-
- for(posA = 0; posA < 2; ++posA, ++iterA) if (*iterA == *e_iter) break;
- for(posB = 0; posB < 2; ++posB, ++iterB) if (*iterB == *e_iter) break;
- if (debug) std::cout << " with face positions " << posA << " and " << posB << std::endl;
- if ((posA == 2) || (posB == 2)) {throw ALE::Exception("Could not find fault face in cone");}
- if (posA == posB) {
- std::cout << "Invalid orientation in fault mesh" << std::endl;
- std::cout << " fault face: " << *e_iter << " cellA: " << cellA << " cellB: " << cellB << std::endl;
- throw ALE::Exception("Invalid orientation in fault mesh");
- }
- } else {
- // Check orientation
- ALE::MinimalArrow<ALE::Mesh::sieve_type::point_type,ALE::Mesh::sieve_type::point_type> arrowA(*e_iter, cellA);
- const int oA = orientation->restrictPoint(arrowA)[0];
- ALE::MinimalArrow<ALE::Mesh::sieve_type::point_type,ALE::Mesh::sieve_type::point_type> arrowB(*e_iter, cellB);
- const int oB = orientation->restrictPoint(arrowB)[0];
-
- if (oA == oB) {
- std::cout << "Invalid orientation in fault mesh" << std::endl;
- std::cout << " fault face: " << *e_iter << " cellA: " << cellA << " cellB: " << cellB << std::endl;
- throw ALE::Exception("Invalid orientation in fault mesh");
- }
- }
- }
- }
- if (debug) fault->view("Oriented Fault mesh");
-}
-
-// ----------------------------------------------------------------------
-void
-pylith::faults::CohesiveTopology::createFault(Obj<SubMesh>& ifault,
- Obj<ALE::Mesh>& faultBd,
- const Obj<Mesh>& mesh,
- const Obj<Mesh::int_section_type>& groupField,
- const bool flipFault)
-{
- const Obj<sieve_type>& sieve = mesh->getSieve();
- const Obj<SubMesh::sieve_type> ifaultSieve = new Mesh::sieve_type(sieve->comm(), sieve->debug());
- Obj<ALE::Mesh> fault = new ALE::Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
- Obj<ALE::Mesh::sieve_type> faultSieve = new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
- const int debug = mesh->debug();
-
// Create set with vertices on fault
- const int_section_type::chart_type& chart = groupField->getChart();
- PointSet faultVertices; // Vertices on fault
+ const IntSection::chart_type& chart = groupField->getChart();
+ TopologyOps::PointSet faultVertices; // Vertices on fault
- for(int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
- assert(!mesh->depth(*c_iter));
- if (groupField->getFiberDimension(*c_iter)) {faultVertices.insert(*c_iter);}
+ const IntSection::chart_type::const_iterator chartEnd = chart.end();
+ for(IntSection::chart_type::const_iterator c_iter = chart.begin();
+ c_iter != chartEnd;
+ ++c_iter) {
+ assert(!sieveMesh->depth(*c_iter));
+ if (groupField->getFiberDimension(*c_iter))
+ faultVertices.insert(*c_iter);
} // for
// Create a sieve which captures the fault
- const bool vertexFault = true;
- const int firstFaultCell = sieve->getBaseSize() + sieve->getCapSize();
+ const bool vertexFault = true;
+ const int firstFaultCell = sieve->getBaseSize() + sieve->getCapSize();
- createFaultSieveFromVertices(fault->getDimension(), firstFaultCell,
- faultVertices, mesh,
- fault->getArrowSection("orientation"),
- faultSieve, flipFault);
+ TopologyOps::createFaultSieveFromVertices(fault->getDimension(), firstFaultCell,
+ faultVertices, sieveMesh,
+ fault->getArrowSection("orientation"),
+ faultSieve, flipFault);
fault->setSieve(faultSieve);
fault->stratify();
- if (debug) fault->view("Fault mesh");
+ if (debug)
+ fault->view("Fault mesh");
- faultBd = ALE::Selection<ALE::Mesh>::boundary(fault);
- if (debug) faultBd->view("Fault boundary mesh");
+ faultBoundary = ALE::Selection<ALE::Mesh>::boundary(fault);
+ if (debug)
+ faultBoundary->view("Fault boundary mesh");
// Orient the fault sieve
- orientFaultSieve(fault->getDimension(), mesh, fault->getArrowSection("orientation"), fault);
+ TopologyOps::orientFaultSieve(fault->getDimension(), sieveMesh,
+ fault->getArrowSection("orientation"), fault);
// Convert fault to an IMesh
- SubMesh::renumbering_type& renumbering = ifault->getRenumbering();
- ifault->setSieve(ifaultSieve);
- ALE::ISieveConverter::convertMesh(*fault, *ifault, renumbering, false);
+ SieveSubMesh::renumbering_type& renumbering = faultSieveMesh->getRenumbering();
+ faultSieveMesh->setSieve(ifaultSieve);
+ ALE::ISieveConverter::convertMesh(*fault, *faultSieveMesh, renumbering, false);
renumbering.clear();
-};
+} // createFault
// ----------------------------------------------------------------------
void
-pylith::faults::CohesiveTopology::create(Obj<SubMesh>& ifault,
- const Obj<ALE::Mesh>& faultBd,
- const Obj<Mesh>& mesh,
- const Obj<Mesh::int_section_type>& groupField,
+pylith::faults::CohesiveTopology::create(topology::Mesh* mesh,
+ const topology::SubMesh& faultMesh,
+ const ALE::Obj<ALE::Mesh>& faultBoundary,
+ const ALE::Obj<topology::Mesh::IntSection>& groupField,
const int materialId,
const bool constraintCell)
{ // create
- typedef ALE::SieveAlg<ALE::Mesh> sieveAlg;
+ assert(0 != mesh);
+ assert(!faultBoundary.isNull());
+ assert(!groupField.isNull());
+
+ typedef ALE::SieveAlg<ALE::Mesh> sieveAlg;
typedef ALE::Selection<ALE::Mesh> selection;
- const Obj<sieve_type>& sieve = mesh->getSieve();
- const Obj<SubMesh::sieve_type> ifaultSieve = ifault->getSieve();
- const int depth = mesh->depth();
- const int numCells = mesh->heightStratum(0)->size();
- int numCorners = 0; // The number of vertices in a mesh cell
- int faceSize = 0; // The number of vertices in a mesh face
- int numFaultCorners = 0; // The number of vertices in a fault cell
- int *indices = NULL; // The indices of a face vertex set in a cell
- const int debug = mesh->debug();
- int oppositeVertex; // For simplices, the vertex opposite a given face
- PointArray origVertices;
- PointArray faceVertices;
- PointArray neighborVertices;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
+ assert(!faultSieveMesh.isNull());
- if (!ifault->commRank()) {
- const SubMesh::point_type p = *ifault->heightStratum(1)->begin();
+ const ALE::Obj<sieve_type>& sieve = sieveMesh->getSieve();
+ assert(!sieve.isNull());
+ const ALE::Obj<SieveSubMesh::sieve_type> ifaultSieve =
+ faultSieveMesh->getSieve();
+ assert(!ifaultSieve.isNull());
- numCorners = mesh->getNumCellCorners();
- faceSize = selection::numFaceVertices(mesh);
- indices = new int[faceSize];
- numFaultCorners = ifault->getNumCellCorners(p, ifault->depth(p));
+ const int depth = sieveMesh->depth();
+ assert(!sieveMesh->heightStratum(0).isNull());
+ const int numCells = sieveMesh->heightStratum(0)->size();
+ int numCorners = 0; // The number of vertices in a mesh cell
+ int faceSize = 0; // The number of vertices in a mesh face
+ int numFaultCorners = 0; // The number of vertices in a fault cell
+ int* indices = 0; // The indices of a face vertex set in a cell
+ const int debug = mesh->debug();
+ int oppositeVertex = 0; // For simplices, the vertex opposite a given face
+ TopologyOps::PointArray origVertices;
+ TopologyOps::PointArray faceVertices;
+ TopologyOps::PointArray neighborVertices;
+
+ if (!faultSieveMesh->commRank()) {
+ assert(!faultSieveMesh->heightStratum(1).isNull());
+ const SieveSubMesh::point_type p = *faultSieveMesh->heightStratum(1)->begin();
+
+ numCorners = sieveMesh->getNumCellCorners();
+ faceSize = selection::numFaceVertices(sieveMesh);
+ indices = new int[faceSize];
+ numFaultCorners = faultSieveMesh->getNumCellCorners(p, faultSieveMesh->depth(p));
}
- //ifault->view("Serial fault mesh");
+ //faultSieveMesh->view("Serial fault mesh");
// Add new shadow vertices and possibly Lagrange multipler vertices
- const Obj<SubMesh::label_sequence>& fVertices = ifault->depthStratum(0);
- const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
- const Obj<std::set<std::string> >& groupNames = mesh->getIntSections();
- Mesh::point_type newPoint = sieve->getBaseSize() + sieve->getCapSize();
- const int numFaultVertices = fVertices->size();
- std::map<Mesh::point_type,Mesh::point_type> vertexRenumber;
- std::map<Mesh::point_type,Mesh::point_type> cellRenumber;
+ const ALE::Obj<SieveSubMesh::label_sequence>& fVertices = faultSieveMesh->depthStratum(0);
+ assert(!fVertices.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const ALE::Obj<std::set<std::string> >& groupNames =
+ sieveMesh->getIntSections();
+ assert(!groupNames.isNull());
+ point_type newPoint = sieve->getBaseSize() + sieve->getCapSize();
+ const int numFaultVertices = fVertices->size();
+ std::map<point_type,point_type> vertexRenumber;
+ std::map<point_type,point_type> cellRenumber;
- for(SubMesh::label_sequence::iterator v_iter = fVertices->begin(); v_iter != fVertices->end(); ++v_iter, ++newPoint) {
+ const SieveSubMesh::label_sequence::const_iterator fVerticesEnd =
+ fVertices->end();
+ for(SieveSubMesh::label_sequence::iterator v_iter = fVertices->begin();
+ v_iter != fVerticesEnd;
+ ++v_iter, ++newPoint) {
vertexRenumber[*v_iter] = newPoint;
- if (debug) std::cout << "Duplicating " << *v_iter << " to " << vertexRenumber[*v_iter] << std::endl;
+ if (debug)
+ std::cout << "Duplicating " << *v_iter << " to "
+ << vertexRenumber[*v_iter] << std::endl;
// Add shadow and constraint vertices (if they exist) to group
// associated with fault
groupField->addPoint(newPoint, 1);
- if (constraintCell) {
+ if (constraintCell)
groupField->addPoint(newPoint+numFaultVertices, 1);
- }
// Add shadow vertices to other groups, don't add constraint
// vertices (if they exist) because we don't want BC, etc to act
// on constraint vertices
+ const std::set<std::string>::const_iterator namesEnd = groupNames->end();
for(std::set<std::string>::const_iterator name = groupNames->begin();
- name != groupNames->end(); ++name) {
- const ALE::Obj<int_section_type>& group = mesh->getIntSection(*name);
+ name != namesEnd;
+ ++name) {
+ const ALE::Obj<IntSection>& group = sieveMesh->getIntSection(*name);
+ assert(!group.isNull());
if (group->getFiberDimension(*v_iter))
group->addPoint(newPoint, 1);
} // for
} // for
+ const std::set<std::string>::const_iterator namesEnd = groupNames->end();
for(std::set<std::string>::const_iterator name = groupNames->begin();
- name != groupNames->end(); ++name) {
- mesh->reallocate(mesh->getIntSection(*name));
+ name != namesEnd;
+ ++name) {
+ sieveMesh->reallocate(sieveMesh->getIntSection(*name));
} // for
#if 0
- for(SubMesh::label_sequence::iterator v_iter = fVertices->begin(); v_iter != fVertices->end(); ++v_iter, ++newPoint) {
+ for(SieveSubMesh::label_sequence::iterator v_iter = fVertices->begin();
+ v_iter != fVerticesEnd;
+ ++v_iter, ++newPoint) {
vertexRenumber[*v_iter] = newPoint;
// OPTIMIZATION
- mesh->setHeight(newPoint, 1);
- mesh->setDepth(newPoint, 0);
+ sieveMesh->setHeight(newPoint, 1);
+ sieveMesh->setDepth(newPoint, 0);
if (constraintCell) {
// OPTIMIZATION
- mesh->setHeight(newPoint+numFaultVertices, 1);
- mesh->setDepth(newPoint+numFaultVertices, 0);
+ sieveMesh->setHeight(newPoint+numFaultVertices, 1);
+ sieveMesh->setDepth(newPoint+numFaultVertices, 0);
}
}
#endif
- if (constraintCell) newPoint += numFaultVertices;
+ if (constraintCell)
+ newPoint += numFaultVertices;
// Split the mesh along the fault sieve and create cohesive elements
- const ALE::Obj<SubMesh::label_sequence>& faces = ifault->heightStratum(1);
- const ALE::Obj<Mesh::label_type>& material = mesh->getLabel("material-id");
+ const ALE::Obj<SieveSubMesh::label_sequence>& faces =
+ faultSieveMesh->heightStratum(1);
+ assert(!faces.isNull());
+ const ALE::Obj<Mesh::label_type>& material =
+ sieveMesh->getLabel("material-id");
+ assert(!material.isNull());
const int firstCohesiveCell = newPoint;
- PointSet replaceCells;
- PointSet noReplaceCells;
- PointSet replaceVertices;
+ TopologyOps::PointSet replaceCells;
+ TopologyOps::PointSet noReplaceCells;
+ TopologyOps::PointSet replaceVertices;
ALE::ISieveVisitor::PointRetriever<sieve_type> sV2(std::max(1, ifaultSieve->getMaxSupportSize()));
- ALE::ISieveVisitor::NConeRetriever<sieve_type> cV2(*ifaultSieve, (size_t) pow(std::max(1, ifaultSieve->getMaxConeSize()), ifault->depth()));
+ ALE::ISieveVisitor::NConeRetriever<sieve_type> cV2(*ifaultSieve, (size_t) pow(std::max(1, ifaultSieve->getMaxConeSize()), faultSieveMesh->depth()));
std::set<Mesh::point_type> faceSet;
- for(SubMesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter, ++newPoint) {
- const Mesh::point_type face = *f_iter;
- if (debug) std::cout << "Considering fault face " << face << std::endl;
+ const SieveSubMesh::label_sequence::const_iterator facesEnd = faces->end();
+ for(SieveSubMesh::label_sequence::iterator f_iter = faces->begin();
+ f_iter != facesEnd;
+ ++f_iter, ++newPoint) {
+ const point_type face = *f_iter;
+ if (debug)
+ std::cout << "Considering fault face " << face << std::endl;
ifaultSieve->support(face, sV2);
- const Mesh::point_type *cells = sV2.getPoints();
- Mesh::point_type cell = cells[0];
- Mesh::point_type otherCell = cells[1];
+ const point_type *cells = sV2.getPoints();
+ point_type cell = cells[0];
+ point_type otherCell = cells[1];
- if (debug) std::cout << " Checking orientation against cell " << cell << std::endl;
- ALE::ISieveTraversal<sieve_type>::orientedClosure(*ifaultSieve, face, cV2);
- const int coneSize = cV2.getSize();
- const Mesh::point_type *faceCone = cV2.getPoints();
+ if (debug)
+ std::cout << " Checking orientation against cell " << cell << std::endl;
+ ALE::ISieveTraversal<SieveMesh::sieve_type>::orientedClosure(*ifaultSieve,
+ face, cV2);
+ const int coneSize = cV2.getSize();
+ const point_type *faceCone = cV2.getPoints();
//ifaultSieve->cone(face, cV2);
- //const int coneSize = cV2.getSize() ? cV2.getSize() : 1;
- //const Mesh::point_type *faceCone = cV2.getSize() ? cV2.getPoints() : &face;
- bool found = true;
+ //const int coneSize = cV2.getSize() ? cV2.getSize() : 1;
+ //const point_type *faceCone = cV2.getSize() ? cV2.getPoints() : &face;
+ bool found = true;
- for(int i = 0; i < coneSize; ++i) faceSet.insert(faceCone[i]);
- selection::getOrientedFace(mesh, cell, &faceSet, numCorners, indices, &origVertices, &faceVertices);
+ for(int i = 0; i < coneSize; ++i)
+ faceSet.insert(faceCone[i]);
+ selection::getOrientedFace(sieveMesh, cell, &faceSet, numCorners, indices,
+ &origVertices, &faceVertices);
if (faceVertices.size() != coneSize) {
- std::cout << "Invalid size for faceVertices " << faceVertices.size() << " of face " << face << "should be " << coneSize << std::endl;
- std::cout << " firstCohesiveCell " << firstCohesiveCell << " newPoint " << newPoint << " numFaces " << faces->size() << std::endl;
+ std::cout << "Invalid size for faceVertices " << faceVertices.size()
+ << " of face " << face << "should be " << coneSize << std::endl;
+ std::cout << " firstCohesiveCell " << firstCohesiveCell << " newPoint "
+ << newPoint << " numFaces " << faces->size() << std::endl;
std::cout << " faceSet:" << std::endl;
- for(std::set<Mesh::point_type>::const_iterator p_iter = faceSet.begin(); p_iter != faceSet.end(); ++p_iter) {
+ for(std::set<Mesh::point_type>::const_iterator p_iter = faceSet.begin();
+ p_iter != faceSet.end();
+ ++p_iter) {
std::cout << " " << *p_iter << std::endl;
- }
+ } // if
std::cout << " cell cone:" << std::endl;
- ALE::ISieveVisitor::PointRetriever<sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
+ ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
sieve->cone(cell, cV);
- const int coneSize2 = cV.getSize();
- const Mesh::point_type *cellCone = cV.getPoints();
+ const int coneSize2 = cV.getSize();
+ const point_type *cellCone = cV.getPoints();
- for(int c = 0; c < coneSize2; ++c) {
+ for(int c = 0; c < coneSize2; ++c)
std::cout << " " << cellCone[c] << std::endl;
- }
std::cout << " fault cell support:" << std::endl;
ALE::ISieveVisitor::PointRetriever<sieve_type> sV(std::max(1, ifaultSieve->getMaxSupportSize()));
ifaultSieve->support(face, sV);
- const int supportSize2 = sV.getSize();
- const Mesh::point_type *cellSupport = sV.getPoints();
- for(int s = 0; s < supportSize2; ++s) {
+ const int supportSize2 = sV.getSize();
+ const point_type *cellSupport = sV.getPoints();
+ for(int s = 0; s < supportSize2; ++s)
std::cout << " " << cellSupport[s] << std::endl;
- }
- }
+ } // if
assert(faceVertices.size() == coneSize);
faceSet.clear();
- ///selection::getOrientedFace(mesh, cell, &vertexRenumber, numCorners, indices, &origVertices, &faceVertices);
+ //selection::getOrientedFace(sieveMesh, cell, &vertexRenumber, numCorners,
+ // indices, &origVertices, &faceVertices);
if (numFaultCorners == 0) {
found = false;
} else if (numFaultCorners == 2) {
- if (faceVertices[0] != faceCone[0]) found = false;
+ if (faceVertices[0] != faceCone[0])
+ found = false;
} else {
int v = 0;
// Locate first vertex
- while((v < numFaultCorners) && (faceVertices[v] != faceCone[0])) ++v;
+ while((v < numFaultCorners) && (faceVertices[v] != faceCone[0]))
+ ++v;
for(int c = 0; c < coneSize; ++c, ++v) {
- if (debug) std::cout << " Checking " << faceCone[c] << " against " << faceVertices[v%numFaultCorners] << std::endl;
+ if (debug) std::cout << " Checking " << faceCone[c] << " against "
+ << faceVertices[v%numFaultCorners] << std::endl;
if (faceVertices[v%numFaultCorners] != faceCone[c]) {
found = false;
break;
- }
- }
- }
+ } // if
+ } // for
+ } // if/else
if (found) {
- if (debug) std::cout << " Choosing other cell" << std::endl;
- Mesh::point_type tmpCell = otherCell;
+ if (debug)
+ std::cout << " Choosing other cell" << std::endl;
+ point_type tmpCell = otherCell;
otherCell = cell;
- cell = tmpCell;
+ cell = tmpCell;
} else {
- if (debug) std::cout << " Verifing reverse orientation" << std::endl;
+ if (debug)
+ std::cout << " Verifing reverse orientation" << std::endl;
found = true;
int v = 0;
if (numFaultCorners > 0) {
// Locate first vertex
- while((v < numFaultCorners) && (faceVertices[v] != faceCone[coneSize-1])) ++v;
+ while((v < numFaultCorners) && (faceVertices[v] != faceCone[coneSize-1]))
+ ++v;
for(int c = coneSize-1; c >= 0; --c, ++v) {
- if (debug) std::cout << " Checking " << faceCone[c] << " against " << faceVertices[v%numFaultCorners] << std::endl;
+ if (debug)
+ std::cout << " Checking " << faceCone[c] << " against "
+ << faceVertices[v%numFaultCorners] << std::endl;
if (faceVertices[v%numFaultCorners] != faceCone[c]) {
found = false;
break;
@@ -653,270 +351,354 @@
replaceVertices.insert(faceCone, &faceCone[coneSize]);
cellRenumber[cell] = newPoint;
// Adding cohesive cell (not interpolated)
- if (debug) std::cout << " Creating cohesive cell " << newPoint << std::endl;
- for(int c = 0; c < coneSize; ++c) {
- if (debug) std::cout << " vertex " << faceCone[c] << std::endl;
+ if (debug)
+ std::cout << " Creating cohesive cell " << newPoint << std::endl;
+ for (int c = 0; c < coneSize; ++c) {
+ if (debug)
+ std::cout << " vertex " << faceCone[c] << std::endl;
sieve->addArrow(faceCone[c], newPoint);
- }
- for(int c = 0; c < coneSize; ++c) {
- if (debug) std::cout << " shadow vertex " << vertexRenumber[faceCone[c]] << std::endl;
+ } // for
+ for (int c = 0; c < coneSize; ++c) {
+ if (debug)
+ std::cout << " shadow vertex " << vertexRenumber[faceCone[c]] << std::endl;
sieve->addArrow(vertexRenumber[faceCone[c]], newPoint);
- }
+ } // for
if (constraintCell) {
- for(int c = 0; c < coneSize; ++c) {
- if (debug) std::cout << " Lagrange vertex " << vertexRenumber[faceCone[c]]+numFaultVertices << std::endl;
+ for (int c = 0; c < coneSize; ++c) {
+ if (debug)
+ std::cout << " Lagrange vertex " << vertexRenumber[faceCone[c]]+numFaultVertices << std::endl;
sieve->addArrow(vertexRenumber[faceCone[c]]+numFaultVertices, newPoint);
- }
- }
+ } // for
+ } // if
// TODO: Need to reform the material label when sieve is reallocated
- mesh->setValue(material, newPoint, materialId);
+ sieveMesh->setValue(material, newPoint, materialId);
#if 0
// OPTIMIZATION
- mesh->setHeight(newPoint, 0);
- mesh->setDepth(newPoint, 1);
+ sieveMesh->setHeight(newPoint, 0);
+ sieveMesh->setDepth(newPoint, 1);
#endif
sV2.clear();
cV2.clear();
} // for
// This completes the set of cells scheduled to be replaced
- PointSet replaceCellsBase(replaceCells);
+ TopologyOps::PointSet replaceCellsBase(replaceCells);
- const ALE::Obj<ALE::Mesh::label_sequence>& faultBdVerts = faultBd->depthStratum(0);
- PointSet faultBdVertices;
+ const ALE::Obj<ALE::Mesh::label_sequence>& faultBdVerts =
+ faultBoundary->depthStratum(0);
+ assert(!faultBdVerts.isNull());
+ TopologyOps::PointSet faultBdVertices;
faultBdVertices.insert(faultBdVerts->begin(), faultBdVerts->end());
- for(PointSet::const_iterator v_iter = replaceVertices.begin(); v_iter != replaceVertices.end(); ++v_iter) {
- if (faultBdVertices.find(*v_iter) != faultBdVertices.end()) continue;
- classifyCells(sieve, *v_iter, depth, faceSize, firstCohesiveCell, replaceCells, noReplaceCells, debug);
- }
- for(PointSet::const_iterator v_iter = faultBdVertices.begin(); v_iter != faultBdVertices.end(); ++v_iter) {
- classifyCells(sieve, *v_iter, depth, faceSize, firstCohesiveCell, replaceCells, noReplaceCells, debug);
- }
+ TopologyOps::PointSet::const_iterator rVerticesEnd = replaceVertices.end();
+ for (TopologyOps::PointSet::const_iterator v_iter = replaceVertices.begin();
+ v_iter != rVerticesEnd; ++v_iter) {
+ if (faultBdVertices.find(*v_iter) != faultBdVertices.end())
+ continue;
+ TopologyOps::classifyCells(sieve, *v_iter, depth, faceSize,
+ firstCohesiveCell, replaceCells, noReplaceCells,
+ debug);
+ } // for
+ const TopologyOps::PointSet::const_iterator fbdVerticesEnd =
+ faultBdVertices.end();
+ for (TopologyOps::PointSet::const_iterator v_iter=faultBdVertices.begin();
+ v_iter != fbdVerticesEnd;
+ ++v_iter) {
+ TopologyOps::classifyCells(sieve, *v_iter, depth, faceSize,
+ firstCohesiveCell, replaceCells, noReplaceCells,
+ debug);
+ } // for
// Add new arrows for support of replaced vertices
- ALE::ISieveVisitor::PointRetriever<sieve_type> sV(std::max(1, sieve->getMaxSupportSize()));
+ ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> sV(std::max(1, sieve->getMaxSupportSize()));
- for(PointSet::const_iterator v_iter = replaceVertices.begin(); v_iter != replaceVertices.end(); ++v_iter) {
+ rVerticesEnd = replaceVertices.end();
+ for (TopologyOps::PointSet::const_iterator v_iter = replaceVertices.begin();
+ v_iter != rVerticesEnd;
+ ++v_iter) {
sieve->support(*v_iter, sV);
- const Mesh::point_type *support = sV.getPoints();
+ const point_type *support = sV.getPoints();
- if (debug) std::cout << " Checking support of " << *v_iter << std::endl;
- for(int s = 0; s < sV.getSize(); ++s) {
+ if (debug)
+ std::cout << " Checking support of " << *v_iter << std::endl;
+ const int sVSize = sV.getSize();
+ for (int s = 0; s < sVSize; ++s) {
if (replaceCells.find(support[s]) != replaceCells.end()) {
- if (debug) std::cout << " Adding new support " << vertexRenumber[*v_iter] << " --> " << support[s] << std::endl;
+ if (debug)
+ std::cout << " Adding new support " << vertexRenumber[*v_iter]
+ << " --> " << support[s] << std::endl;
sieve->addArrow(vertexRenumber[*v_iter], support[s]);
} else {
- if (debug) std::cout << " Keeping same support " << *v_iter<<","<<vertexRenumber[*v_iter] << " --> " << support[s] << std::endl;
- }
- }
+ if (debug)
+ std::cout << " Keeping same support " << *v_iter<<","
+ << vertexRenumber[*v_iter] << " --> " << support[s]
+ << std::endl;
+ } // if/else
+ } // for
sV.clear();
}
sieve->reallocate();
// More checking
- const bool firstFault = !mesh->hasRealSection("replaced_cells");
- const ALE::Obj<real_section_type>& replacedCells = mesh->getRealSection("replaced_cells");
- PointSet cellNeighbors;
-
+ const bool firstFault = !sieveMesh->hasRealSection("replaced_cells");
+ const ALE::Obj<topology::Mesh::RealSection>& replacedCells =
+ sieveMesh->getRealSection("replaced_cells");
+ assert(!replacedCells.isNull());
+ TopologyOps::PointSet cellNeighbors;
+
if (firstFault) {
- const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->heightStratum(0);
+ assert(!cells.isNull());
- replacedCells->setChart(real_section_type::chart_type(*std::min_element(cells->begin(), cells->end()), *std::max_element(cells->begin(), cells->end())+1));
+ replacedCells->setChart(topology::Mesh::RealSection::chart_type(*std::min_element(cells->begin(), cells->end()), *std::max_element(cells->begin(), cells->end())+1));
replacedCells->setFiberDimension(cells, 1);
replacedCells->allocatePoint();
- }
- for(PointSet::const_iterator c_iter = noReplaceCells.begin(); c_iter != noReplaceCells.end(); ++c_iter) {
+ } // if
+
+ const TopologyOps::PointSet::const_iterator noRCellsEnd = noReplaceCells.end();
+ for (TopologyOps::PointSet::const_iterator c_iter = noReplaceCells.begin();
+ c_iter != noRCellsEnd;
+ ++c_iter) {
const double minusOne = -1.0;
-
if (replacedCells->restrictPoint(*c_iter)[0] == 0.0) {
replacedCells->updatePoint(*c_iter, &minusOne);
} else {
const double minusTwo = -2.0;
-
replacedCells->updatePoint(*c_iter, &minusTwo);
- }
- }
- for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
+ } // if/else
+ } // for
+
+ TopologyOps::PointSet::const_iterator rCellsEnd = replaceCells.end();
+ for (TopologyOps::PointSet::const_iterator c_iter = replaceCells.begin();
+ c_iter != rCellsEnd;
+ ++c_iter) {
if (replaceCellsBase.find(*c_iter) != replaceCellsBase.end()) {
const double one = 1.0;
-
if (replacedCells->restrictPoint(*c_iter)[0] == 0.0) {
replacedCells->updatePoint(*c_iter, &one);
} else {
const double two = 2.0;
-
replacedCells->updatePoint(*c_iter, &two);
- }
+ } // if/else
continue;
- }
+ } // if
const double ten = 10.0;
-
if (replacedCells->restrictPoint(*c_iter)[0] == 0.0) {
replacedCells->updatePoint(*c_iter, &ten);
} else {
- const double twenty = 20.0;
-
- replacedCells->updatePoint(*c_iter, &twenty);
- }
+ const double twenty = 20.0;
+ replacedCells->updatePoint(*c_iter, &twenty);
+ } // if/else
// There should be a way to check for boundary elements
- if (mesh->getDimension() == 1) {
+ if (mesh->dimension() == 1) {
if (cellNeighbors.size() > 2) {
- std::cout << "Cell " << *c_iter << " has an invalid number of neighbors " << cellNeighbors.size() << std::endl;
+ std::cout << "Cell " << *c_iter
+ << " has an invalid number of neighbors "
+ << cellNeighbors.size() << std::endl;
throw ALE::Exception("Invalid number of neighbors");
- }
- } else if (mesh->getDimension() == 2) {
+ } // if
+ } else if (mesh->dimension() == 2) {
if (numCorners == 3) {
if (cellNeighbors.size() > 3) {
- std::cout << "Cell " << *c_iter << " has an invalid number of neighbors " << cellNeighbors.size() << std::endl;
+ std::cout << "Cell " << *c_iter
+ << " has an invalid number of neighbors "
+ << cellNeighbors.size() << std::endl;
throw ALE::Exception("Invalid number of neighbors");
- }
+ } // if
} else if (numCorners == 4) {
if (cellNeighbors.size() > 4) {
- std::cout << "Cell " << *c_iter << " has an invalid number of neighbors " << cellNeighbors.size() << std::endl;
+ std::cout << "Cell " << *c_iter
+ << " has an invalid number of neighbors "
+ << cellNeighbors.size() << std::endl;
throw ALE::Exception("Invalid number of neighbors");
- }
- }
- } else if (mesh->getDimension() == 3) {
+ } // if
+ } // if/else
+ } else if (mesh->dimension() == 3) {
if (numCorners == 4) {
if (cellNeighbors.size() > 4) {
- std::cout << "Cell " << *c_iter << " has an invalid number of neighbors " << cellNeighbors.size() << std::endl;
+ std::cout << "Cell " << *c_iter
+ << " has an invalid number of neighbors "
+ << cellNeighbors.size() << std::endl;
throw ALE::Exception("Invalid number of neighbors");
- }
+ } // if
} else if (numCorners == 8) {
if (cellNeighbors.size() > 6) {
- std::cout << "Cell " << *c_iter << " has an invalid number of neighbors " << cellNeighbors.size() << std::endl;
+ std::cout << "Cell " << *c_iter
+ << " has an invalid number of neighbors "
+ << cellNeighbors.size() << std::endl;
throw ALE::Exception("Invalid number of neighbors");
- }
- }
- }
- }
+ } // if
+ } // if/else
+ } // if/else
+ } // for
ReplaceVisitor<sieve_type,std::map<Mesh::point_type,Mesh::point_type> > rVc(vertexRenumber, std::max(1, sieve->getMaxConeSize()), debug);
-
- for(PointSet::const_iterator c_iter = replaceCells.begin(); c_iter != replaceCells.end(); ++c_iter) {
+
+ rCellsEnd = replaceCells.end();
+ for (TopologyOps::PointSet::const_iterator c_iter = replaceCells.begin();
+ c_iter != rCellsEnd;
+ ++c_iter) {
sieve->cone(*c_iter, rVc);
if (rVc.mappedPoint()) {
- if (debug) std::cout << " Replacing cell " << *c_iter << std::endl;
+ if (debug)
+ std::cout << " Replacing cell " << *c_iter << std::endl;
sieve->setCone(rVc.getPoints(), *c_iter);
- }
+ } // if
rVc.clear();
- }
+ } // for
ReplaceVisitor<sieve_type,std::map<Mesh::point_type,Mesh::point_type> > rVs(cellRenumber, std::max(1, sieve->getMaxSupportSize()), debug);
- for(PointSet::const_iterator v_iter = replaceVertices.begin(); v_iter != replaceVertices.end(); ++v_iter) {
+ rVerticesEnd = replaceVertices.end();
+ for (TopologyOps::PointSet::const_iterator v_iter = replaceVertices.begin();
+ v_iter != rVerticesEnd;
+ ++v_iter) {
sieve->support(*v_iter, rVs);
if (rVs.mappedPoint()) {
- if (debug) std::cout << " Replacing support for " << *v_iter << std::endl;
+ if (debug)
+ std::cout << " Replacing support for " << *v_iter << std::endl;
sieve->setSupport(*v_iter, rVs.getPoints());
} else {
- if (debug) std::cout << " Not replacing support for " << *v_iter << std::endl;
- }
+ if (debug)
+ std::cout << " Not replacing support for " << *v_iter << std::endl;
+ } // if/else
rVs.clear();
- }
- if (!ifault->commRank()) delete [] indices;
+ } // for
+ if (!faultSieveMesh->commRank())
+ delete [] indices;
#if 1
- mesh->stratify();
+ sieveMesh->stratify();
#endif
const std::string labelName("censored depth");
- if (!mesh->hasLabel(labelName)) {
- const ALE::Obj<Mesh::label_type>& label = mesh->createLabel(labelName);
+ if (!sieveMesh->hasLabel(labelName)) {
+ const ALE::Obj<Mesh::label_type>& label = sieveMesh->createLabel(labelName);
+ assert(!label.isNull());
- _computeCensoredDepth(label, mesh->getSieve(), firstCohesiveCell-(constraintCell?numFaultVertices:0));
+ TopologyOps::computeCensoredDepth(label, sieveMesh->getSieve(),
+ firstCohesiveCell-(constraintCell?numFaultVertices:0));
} else {
// Insert new shadow vertices into existing label
- const ALE::Obj<Mesh::label_type>& label = mesh->getLabel(labelName);
+ const ALE::Obj<Mesh::label_type>& label = sieveMesh->getLabel(labelName);
+ assert(!label.isNull());
- for(std::map<int,int>::const_iterator v_iter = vertexRenumber.begin(); v_iter != vertexRenumber.end(); ++v_iter) {
- mesh->setValue(label, v_iter->second, 0);
- }
- }
- if (debug) mesh->view("Mesh with Cohesive Elements");
+ const std::map<int,int>::const_iterator vRenumberEnd = vertexRenumber.end();
+ for (std::map<int,int>::const_iterator v_iter = vertexRenumber.begin();
+ v_iter != vRenumberEnd;
+ ++v_iter)
+ sieveMesh->setValue(label, v_iter->second, 0);
+ } // if/else
+ if (debug)
+ mesh->view("Mesh with Cohesive Elements");
// Fix coordinates
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- const ALE::Obj<SubMesh::label_sequence>& fVertices2 = ifault->depthStratum(0);
+ const ALE::Obj<topology::Mesh::RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& fVertices2 =
+ faultSieveMesh->depthStratum(0);
+ assert(!fVertices2.isNull());
- if (debug) coordinates->view("Coordinates without shadow vertices");
- for(SubMesh::label_sequence::iterator v_iter = fVertices2->begin();
- v_iter != fVertices2->end();
+ if (debug)
+ coordinates->view("Coordinates without shadow vertices");
+ SieveSubMesh::label_sequence::const_iterator fVertices2End =
+ fVertices2->end();
+ for (SieveSubMesh::label_sequence::iterator v_iter = fVertices2->begin();
+ v_iter != fVertices2End;
++v_iter) {
coordinates->addPoint(vertexRenumber[*v_iter],
coordinates->getFiberDimension(*v_iter));
- if (constraintCell) {
+ if (constraintCell)
coordinates->addPoint(vertexRenumber[*v_iter]+numFaultVertices,
- coordinates->getFiberDimension(*v_iter));
- }
+ coordinates->getFiberDimension(*v_iter));
} // for
- mesh->reallocate(coordinates);
- for(SubMesh::label_sequence::iterator v_iter = fVertices2->begin();
- v_iter != fVertices2->end();
+ sieveMesh->reallocate(coordinates);
+ fVertices2End = fVertices2->end();
+ for (SieveSubMesh::label_sequence::iterator v_iter = fVertices2->begin();
+ v_iter != fVertices2End;
++v_iter) {
coordinates->updatePoint(vertexRenumber[*v_iter],
coordinates->restrictPoint(*v_iter));
- if (constraintCell) {
+ if (constraintCell)
coordinates->updatePoint(vertexRenumber[*v_iter]+numFaultVertices,
coordinates->restrictPoint(*v_iter));
- }
- }
- if (debug) coordinates->view("Coordinates with shadow vertices");
+ } // for
+ if (debug)
+ coordinates->view("Coordinates with shadow vertices");
} // createCohesiveCells
// ----------------------------------------------------------------------
// Form a parallel fault mesh using the cohesive cell information
void
-pylith::faults::CohesiveTopology::createParallel(
- ALE::Obj<SubMesh>* ifault,
- std::map<Mesh::point_type, Mesh::point_type>* cohesiveToFault,
- const ALE::Obj<Mesh>& mesh,
- const int materialId,
- const bool constraintCell)
-{
- assert(0 != ifault);
+pylith::faults::CohesiveTopology::createFaultParallel(
+ topology::SubMesh* faultMesh,
+ std::map<point_type, point_type>* cohesiveToFault,
+ const topology::Mesh& mesh,
+ const int materialId,
+ const bool constraintCell)
+{ // createFaultParallel
+ assert(0 != faultMesh);
assert(0 != cohesiveToFault);
- const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
- *ifault = new SubMesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
- 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());
+ faultMesh->coordsys(mesh);
+
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+
+ const ALE::Obj<sieve_type>& sieve = sieveMesh->getSieve();
+ assert(!sieve.isNull());
+ faultSieveMesh =
+ new SieveSubMesh(mesh.comm(), mesh.dimension()-1, mesh.debug());
+ const ALE::Obj<SieveMesh::sieve_type> ifaultSieve =
+ new SieveMesh::sieve_type(sieve->comm(), sieve->debug());
+ assert(!ifaultSieve.isNull());
+ ALE::Obj<ALE::Mesh> fault =
+ new ALE::Mesh(mesh.comm(), mesh.dimension()-1, mesh.debug());
+ assert(!fault.isNull());
+ ALE::Obj<ALE::Mesh::sieve_type> faultSieve =
+ new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
+ assert(!faultSieve.isNull());
cohesiveToFault->clear();
- const ALE::Obj<Mesh::label_sequence>& cohesiveCells = mesh->getLabelStratum("material-id", materialId);
- const Mesh::label_sequence::iterator cBegin = cohesiveCells->begin();
- const Mesh::label_sequence::iterator cEnd = cohesiveCells->end();
+ const ALE::Obj<SieveMesh::label_sequence>& cohesiveCells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ assert(!cohesiveCells.isNull());
+ const SieveMesh::label_sequence::iterator cBegin = cohesiveCells->begin();
+ const SieveMesh::label_sequence::iterator cEnd = cohesiveCells->end();
const int sieveEnd = sieve->getBaseSize() + sieve->getCapSize();
const int numFaces = cohesiveCells->size();
int globalSieveEnd = 0;
int globalFaceOffset = 0;
- MPI_Allreduce((void *) &sieveEnd, (void *) &globalSieveEnd, 1, MPI_INT, MPI_SUM, sieve->comm());
- MPI_Scan((void *) &numFaces, (void *) &globalFaceOffset, 1, MPI_INT, MPI_SUM, sieve->comm());
+ MPI_Allreduce((void *) &sieveEnd, (void *) &globalSieveEnd, 1,
+ MPI_INT, MPI_SUM, sieve->comm());
+ MPI_Scan((void *) &numFaces, (void *) &globalFaceOffset, 1,
+ MPI_INT, MPI_SUM, sieve->comm());
int face = globalSieveEnd + globalFaceOffset - numFaces;
ALE::ISieveVisitor::PointRetriever<sieve_type> cV(sieve->getMaxConeSize());
- for(Mesh::label_sequence::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter) {
+ for(SieveMesh::label_sequence::iterator c_iter = cBegin;
+ c_iter != cEnd;
+ ++c_iter) {
sieve->cone(*c_iter, cV);
- const int coneSize = cV.getSize();
- const Mesh::point_type *cone = cV.getPoints();
- int color = 0;
+ const int coneSize = cV.getSize();
+ const Mesh::point_type *cone = cV.getPoints();
+ int color = 0;
if (!constraintCell) {
const int faceSize = coneSize / 2;
assert(0 == coneSize % faceSize);
// Use first vertices (negative side of the fault) for fault mesh
- for(int i = 0; i < faceSize; ++i) {
+ for (int i = 0; i < faceSize; ++i)
faultSieve->addArrow(cone[i], face, color++);
- }
} else {
const int faceSize = coneSize / 3;
assert(0 == coneSize % faceSize);
// Use last vertices (contraints) for fault mesh
- for(int i = 2*faceSize; i < 3*faceSize; ++i) {
+ for (int i = 2*faceSize; i < 3*faceSize; ++i)
faultSieve->addArrow(cone[i], face, color++);
- }
} // if/else
(*cohesiveToFault)[*c_iter] = face;
++face;
@@ -926,167 +708,98 @@
fault->stratify();
// Convert fault to an IMesh
- SubMesh::renumbering_type& fRenumbering = (*ifault)->getRenumbering();
- (*ifault)->setSieve(ifaultSieve);
- //ALE::ISieveConverter::convertMesh(*fault, *(*ifault), fRenumbering, true);
+ SieveSubMesh::renumbering_type& fRenumbering =
+ faultSieveMesh->getRenumbering();
+ faultSieveMesh->setSieve(ifaultSieve);
+ //ALE::ISieveConverter::convertMesh(*fault, *faultSieveMesh, 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());
+ ALE::ISieveConverter::convertSieve(*fault->getSieve(),
+ *faultSieveMesh->getSieve(),
+ fRenumbering, true);
+ faultSieveMesh->stratify();
+ ALE::ISieveConverter::convertOrientation(*fault->getSieve(),
+ *faultSieveMesh->getSieve(),
+ fRenumbering,
+ fault->getArrowSection("orientation").ptr());
}
fault = NULL;
faultSieve = NULL;
- const ALE::Obj<SubMesh::label_sequence>& faultCells = (*ifault)->heightStratum(0);
+ const ALE::Obj<SieveSubMesh::label_sequence>& faultCells =
+ faultSieveMesh->heightStratum(0);
assert(!faultCells.isNull());
- SubMesh::label_sequence::iterator f_iter = faultCells->begin();
+ SieveSubMesh::label_sequence::iterator f_iter = faultCells->begin();
- for(Mesh::label_sequence::iterator c_iter = cBegin; c_iter != cEnd; ++c_iter, ++f_iter) {
+ for(SieveMesh::label_sequence::iterator c_iter = cBegin;
+ c_iter != cEnd;
+ ++c_iter, ++f_iter)
(*cohesiveToFault)[*c_iter] = *f_iter;
- }
#if 0
- (*ifault)->setRealSection("coordinates", mesh->getRealSection("coordinates"));
+ faultSieveMesh->setRealSection("coordinates",
+ sieveMesh->getRealSection("coordinates"));
#else
- const ALE::Obj<Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
- 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();
+ const ALE::Obj<topology::Mesh::RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ const ALE::Obj<topology::Mesh::RealSection>& fCoordinates =
+ faultSieveMesh->getRealSection("coordinates");
+ assert(!fCoordinates.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ const SieveMesh::label_sequence::iterator vBegin = vertices->begin();
+ const SieveMesh::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;
+ fCoordinates->setChart(topology::Mesh::RealSection::chart_type(faultSieveMesh->heightStratum(0)->size(),
+ faultSieveMesh->getSieve()->getChart().max()));
+ for (SieveMesh::label_sequence::iterator v_iter = vBegin;
v_iter != vEnd;
++v_iter) {
- if (fRenumbering.find(*v_iter) == fRenumbering.end()) continue;
- fCoordinates->setFiberDimension(fRenumbering[*v_iter], coordinates->getFiberDimension(*v_iter));
- }
+ if (fRenumbering.find(*v_iter) == fRenumbering.end())
+ continue;
+ fCoordinates->setFiberDimension(fRenumbering[*v_iter],
+ coordinates->getFiberDimension(*v_iter));
+ } // for
fCoordinates->allocatePoint();
- for(Mesh::label_sequence::iterator v_iter = vBegin;
+ for(SieveMesh::label_sequence::iterator v_iter = vBegin;
v_iter != vEnd;
++v_iter) {
if (fRenumbering.find(*v_iter) == fRenumbering.end()) continue;
- fCoordinates->updatePoint(fRenumbering[*v_iter], coordinates->restrictPoint(*v_iter));
+ fCoordinates->updatePoint(fRenumbering[*v_iter],
+ coordinates->restrictPoint(*v_iter));
}
#endif
- //(*ifault)->view("Parallel fault mesh");
+ //faultSieveMesh->view("Parallel fault mesh");
// Create the parallel overlap
// Can I figure this out in a nicer way?
- Obj<SubMesh::send_overlap_type> sendParallelMeshOverlap = (*ifault)->getSendOverlap();
- Obj<SubMesh::recv_overlap_type> recvParallelMeshOverlap = (*ifault)->getRecvOverlap();
+ ALE::Obj<SieveSubMesh::send_overlap_type> sendParallelMeshOverlap =
+ faultSieveMesh->getSendOverlap();
+ assert(!sendParallelMeshOverlap.isNull());
+ ALE::Obj<SieveSubMesh::recv_overlap_type> recvParallelMeshOverlap =
+ faultSieveMesh->getRecvOverlap();
+ assert(!recvParallelMeshOverlap.isNull());
// Must process the renumbering local --> fault to global --> fault
- Mesh::renumbering_type& renumbering = mesh->getRenumbering();
- Mesh::renumbering_type gRenumbering;
+ SieveMesh::renumbering_type& renumbering = sieveMesh->getRenumbering();
+ SieveMesh::renumbering_type gRenumbering;
- for(Mesh::renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter) {
- if (fRenumbering.find(r_iter->second) != fRenumbering.end()) {
+ const SieveMesh::renumbering_type::const_iterator renumberingEnd =
+ renumbering.end();
+ for (SieveMesh::renumbering_type::const_iterator r_iter = renumbering.begin();
+ r_iter != renumberingEnd;
+ ++r_iter)
+ if (fRenumbering.find(r_iter->second) != fRenumbering.end())
gRenumbering[r_iter->first] = fRenumbering[r_iter->second];
- }
- }
- ALE::SetFromMap<Mesh::renumbering_type> globalPoints(gRenumbering);
- ALE::OverlapBuilder<>::constructOverlap(globalPoints, gRenumbering, sendParallelMeshOverlap, recvParallelMeshOverlap);
- (*ifault)->setCalculatedOverlap(true);
+ ALE::SetFromMap<SieveMesh::renumbering_type> globalPoints(gRenumbering);
+ ALE::OverlapBuilder<>::constructOverlap(globalPoints, gRenumbering,
+ sendParallelMeshOverlap,
+ recvParallelMeshOverlap);
+ faultSieveMesh->setCalculatedOverlap(true);
//sendParallelMeshOverlap->view("Send parallel fault overlap");
//recvParallelMeshOverlap->view("Recv parallel fault overlap");
}
-// ----------------------------------------------------------------------
-void
-pylith::faults::CohesiveTopology::classifyCells(const ALE::Obj<Mesh::sieve_type>& sieve,
- const Mesh::point_type& vertex,
- const int depth,
- const int faceSize,
- const Mesh::point_type& firstCohesiveCell,
- PointSet& replaceCells,
- PointSet& noReplaceCells,
- const int debug)
-{
- // Replace all cells on a given side of the fault with a vertex on the fault
- ClassifyVisitor<Mesh::sieve_type> cV(*sieve, replaceCells, noReplaceCells, firstCohesiveCell, faceSize, debug);
- const PointSet& vReplaceCells = cV.getReplaceCells();
- const PointSet& vNoReplaceCells = cV.getNoReplaceCells();
- if (debug) {std::cout << "Checking fault vertex " << vertex << std::endl;}
- sieve->support(vertex, cV);
- cV.setMode(false);
- const int classifyTotal = cV.getSize();
- int classifySize = vReplaceCells.size() + vNoReplaceCells.size();
-
- while(cV.getModified() && (classifySize < classifyTotal)) {
- cV.reset();
- sieve->support(vertex, cV);
- if (debug) {
- std::cout << "classifySize: " << classifySize << std::endl;
- std::cout << "classifyTotal: " << classifyTotal << std::endl;
- std::cout << "vReplaceCells.size: " << vReplaceCells.size() << std::endl;
- std::cout << "vNoReplaceCells.size: " << vNoReplaceCells.size() << std::endl;
- }
- assert(classifySize < vReplaceCells.size() + vNoReplaceCells.size());
- classifySize = vReplaceCells.size() + vNoReplaceCells.size();
- assert(classifySize <= classifyTotal);
- }
- replaceCells.insert(vReplaceCells.begin(), vReplaceCells.end());
- // More checking
- noReplaceCells.insert(vNoReplaceCells.begin(), vNoReplaceCells.end());
-}
-
-// ----------------------------------------------------------------------
-template<class InputPoints>
-bool
-pylith::faults::CohesiveTopology::_compatibleOrientation(const ALE::Obj<Mesh>& mesh,
- const Mesh::point_type& p,
- const Mesh::point_type& q,
- const int numFaultCorners,
- const int faultFaceSize,
- const int faultDepth,
- const Obj<InputPoints>& points,
- int indices[],
- PointArray *origVertices,
- PointArray *faceVertices,
- PointArray *neighborVertices)
-{
- typedef ALE::Selection<Mesh> selection;
- const int debug = mesh->debug();
- bool compatible;
-
- bool eOrient = selection::getOrientedFace(mesh, p, points, numFaultCorners, indices, origVertices, faceVertices);
- bool nOrient = selection::getOrientedFace(mesh, q, points, numFaultCorners, indices, origVertices, neighborVertices);
-
- if (faultFaceSize > 1) {
- if (debug) {
- for(PointArray::iterator v_iter = faceVertices->begin(); v_iter != faceVertices->end(); ++v_iter) {
- std::cout << " face vertex " << *v_iter << std::endl;
- }
- for(PointArray::iterator v_iter = neighborVertices->begin(); v_iter != neighborVertices->end(); ++v_iter) {
- std::cout << " neighbor vertex " << *v_iter << std::endl;
- }
- }
- compatible = !(*faceVertices->begin() == *neighborVertices->begin());
- } else {
- compatible = !(nOrient == eOrient);
- }
- return compatible;
-}
-
-// ----------------------------------------------------------------------
-void
-pylith::faults::CohesiveTopology::_computeCensoredDepth(const ALE::Obj<Mesh::label_type>& depth,
- const ALE::Obj<Mesh::sieve_type>& sieve,
- const Mesh::point_type& firstCohesiveCell)
-{
- Mesh::DepthVisitor d(*sieve, firstCohesiveCell, *depth);
-
- sieve->roots(d);
- while(d.isModified()) {
- // FIX: Avoid the copy here somehow by fixing the traversal
- std::vector<Mesh::point_type> modifiedPoints(d.getModifiedPoints().begin(), d.getModifiedPoints().end());
-
- d.clear();
- sieve->support(modifiedPoints, d);
- }
-};
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.hh 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/CohesiveTopology.hh 2009-04-06 03:55:49 UTC (rev 14601)
@@ -21,29 +21,32 @@
// Include directives ---------------------------------------------------
#include "faultsfwd.hh" // forward declarations
+#include "pylith/topology/Mesh.hh" // USES Mesh::IntSection
+#include "pylith/utils/sievetypes.hh" // USE ALE::Obj
+
// CohesiveTopology -----------------------------------------------------
class pylith::faults::CohesiveTopology
{ // class CohesiveTopology
-public :
- typedef std::set<SieveMesh::point_type> PointSet;
- typedef std::vector<sieve_type::point_type> PointArray;
- typedef std::pair<sieve_type::point_type, int> oPoint_type;
- typedef std::vector<oPoint_type> oPointArray;
+private :
+ typedef pylith::topology::Mesh::SieveMesh::point_type point_type;
+
// PUBLIC METHODS /////////////////////////////////////////////////////
public :
+
/** Create the fault mesh.
*
- * @param fault Finite-element mesh of fault (output)
- * @param mesh Finite-element mesh
+ * @param faultMesh Finite-element mesh of fault (output).
+ * @param faultBoundary Finite-element mesh of fault boundary (output).
+ * @param mesh Finite-element mesh of domain.
* @param faultVertices Vertices assocated with faces of cells defining
* fault surface
*/
static
- void createFault(topology::SubMesh* ifault,
- ALE::Obj<ALE::Mesh>& faultBd,
- const topology::Mesh& mesh,
- const ALE::Obj<topology::Mesh::IntSection>& groupField,
+ void createFault(topology::SubMesh* faultMesh,
+ ALE::Obj<ALE::Mesh>& faultBoundary,
+ const topology::Mesh& mesh,
+ const ALE::Obj<topology::Mesh::IntSection>& groupField,
const bool flipFault =false);
/** Create cohesive cells.
@@ -55,28 +58,29 @@
* Lagrange multipliers that require extra vertices, false otherwise
*/
static
- void create(topology::SubMesh* ifault,
- const ALE::Obj<ALE::Mesh>& faultBd,
- const topology::Mesh& mesh,
+ void create(topology::Mesh* mesh,
+ const topology::SubMesh& faultMesh,
+ const ALE::Obj<ALE::Mesh>& faultBoundary,
const ALE::Obj<topology::Mesh::IntSection>& groupField,
const int materialId,
const bool constraintCell =false);
/** Create (distributed) fault mesh from cohesive cells.
*
- * @param fault Finite-element mesh of fault (output).
- * @param cohesiveToFault Mapping of cohesive cell to fault mesh cell.
+ * @param faultMesh Finite-element mesh of fault (output).
+ * @param cohesiveToFault Mapping of cohesive cell to fault mesh
+ * cell (output).
* @param mesh Finite-element mesh.
* @param materialId Material id for cohesive elements.
* @param constraintCell True if creating cells constrained with
* Lagrange multipliers that require extra vertices, false otherwise.
*/
static
- void createParallel(topology::SubMesh* ifault,
- std::map<Mesh::point_type, Mesh::point_type>* cohesiveToFault,
- const topology::Mesh& mesh,
- const int materialId,
- const bool constraintCell =false);
+ void createFaultParallel(topology::SubMesh* faultMesh,
+ std::map<point_type, point_type>* cohesiveToFault,
+ const topology::Mesh& mesh,
+ const int materialId,
+ const bool constraintCell =false);
}; // class CohesiveTopology
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.hh 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/Fault.hh 2009-04-06 03:55:49 UTC (rev 14601)
@@ -78,7 +78,7 @@
* @param mesh PETSc mesh
*/
virtual
- void adjustTopology(const topology::Mesh& mesh,
+ void adjustTopology(topology::Mesh* mesh,
const bool flipFault =false) = 0;
/** Initialize fault. Determine orientation and setup boundary
@@ -127,7 +127,7 @@
virtual
const topology::Field<topology::SubMesh>&
cellField(const char* name,
- topology::SolutionFields& fields) = 0;
+ const topology::SolutionFields& fields) = 0;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -14,11 +14,10 @@
#include "FaultCohesive.hh" // implementation of object methods
-//#include "CohesiveTopology.hh" // USES CohesiveTopology::create()
+#include "CohesiveTopology.hh" // USES CohesiveTopology
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/meshio/UCDFaultFile.hh" // USES UCDFaultFile
-//#include "pylith/meshio/UCDFaultFile.hh" // USES UCDFaultFile::read()
-#include "pylith/utils/array.hh" // USES double_array
-
#include <cassert> // USES assert()
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
@@ -57,31 +56,37 @@
// ----------------------------------------------------------------------
// Adjust mesh topology for fault implementation.
void
-pylith::faults::FaultCohesive::adjustTopology(const topology::Mesh& mesh,
+pylith::faults::FaultCohesive::adjustTopology(topology::Mesh* const mesh,
const bool flipFault)
{ // adjustTopology
-#if 0
+ assert(0 != mesh);
assert(std::string("") != label());
- SubMesh faultMesh;
- Mesh fault
- Obj<SubMesh> faultMesh = NULL;
- Obj<ALE::Mesh> faultBd = NULL;
+ topology::SubMesh faultMesh;
+ ALE::Obj<ALE::Mesh> faultBoundary;
+
+ // Get group of vertices associated with fault
+ const ALE::Obj<topology::Mesh::SieveMesh>& sieveMesh = mesh->sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<topology::Mesh::IntSection>& groupField =
+ sieveMesh->getIntSection(label());
+ assert(!groupField.isNull());
+
if (_useFaultMesh) {
const int faultDim = 2;
- assert(3 == mesh.dimension());
+ assert(3 == mesh->dimension());
- //MPI_Bcast(&faultDim, 1, MPI_INT, 0, comm);
- faultMesh = new SubMesh(mesh->comm(), faultDim, mesh->debug());
- meshio::UCDFaultFile::readFault(_faultMeshFilename, mesh,
- faultMesh, faultBd);
+ meshio::UCDFaultFile::read(_faultMeshFilename.c_str(),
+ &faultMesh, faultBoundary, *mesh);
- // Get group of vertices associated with fault
- const ALE::Obj<IntSection>& groupField = sieveMesh->getIntSection(label());
- faultMesh->setRealSection("coordinates",
- mesh->getRealSection("coordinates"));
+ // Set coordinates in fault mesh
+ const ALE::Obj<topology::SubMesh::SieveMesh>& faultSieveMesh =
+ faultMesh.sieveMesh();
+ assert(!faultSieveMesh.isNull());
+ faultSieveMesh->setRealSection("coordinates",
+ sieveMesh->getRealSection("coordinates"));
- CohesiveTopology::create(faultMesh, faultBd, mesh, groupField, id(),
+ CohesiveTopology::create(mesh, faultMesh, faultBoundary, groupField, id(),
_useLagrangeConstraints());
} else {
if (!sieveMesh->hasIntSection(label())) {
@@ -91,20 +96,12 @@
throw std::runtime_error(msg.str());
} // if
- // Get group of vertices associated with fault
- const ALE::Obj<IntSection>& groupField = sieveMesh->getIntSection(label());
- assert(!groupField.isNull());
-
- faultMesh =
- new SubMesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
-
- CohesiveTopology::createFault(faultMesh, faultBd, mesh, groupField,
+ CohesiveTopology::createFault(&faultMesh, faultBoundary, *mesh, groupField,
flipFault);
- CohesiveTopology::create(faultMesh, faultBd, mesh, groupField, id(),
+ CohesiveTopology::create(mesh, faultMesh, faultBoundary, groupField, id(),
_useLagrangeConstraints());
} // if/else
-#endif
} // adjustTopology
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.hh 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesive.hh 2009-04-06 03:55:49 UTC (rev 14601)
@@ -56,7 +56,7 @@
* @param mesh PETSc mesh.
* @param flipFault Flip fault orientation.
*/
- void adjustTopology(const topology::Mesh& mesh,
+ void adjustTopology(topology::Mesh* const mesh,
const bool flipFault =false);
// PROTECTED METHODS //////////////////////////////////////////////////
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -41,7 +41,8 @@
// ----------------------------------------------------------------------
// Default constructor.
-pylith::faults::FaultCohesiveKin::FaultCohesiveKin(void)
+pylith::faults::FaultCohesiveKin::FaultCohesiveKin(void) :
+ _fields(0)
{ // constructor
} // constructor
@@ -49,7 +50,8 @@
// Destructor.
pylith::faults::FaultCohesiveKin::~FaultCohesiveKin(void)
{ // destructor
- // Don't manage memory for eq source pointers
+ delete _fields; _fields = 0;
+ // :TODO: Use shared pointers for earthquake sources
} // destructor
// ----------------------------------------------------------------------
@@ -59,24 +61,27 @@
EqKinSrc** sources,
const int numSources)
{ // eqsrcs
+ // :TODO: Use shared pointers for earthquake sources
_eqSrcs.clear();
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]; // Don't manage memory for eq source
+ _eqSrcs[std::string(names[i])] = sources[i];
} // for
} // eqsrcs
// ----------------------------------------------------------------------
// Initialize fault. Determine orientation and setup boundary
void
-pylith::faults::FaultCohesiveKin::initialize(const ALE::Obj<Mesh>& mesh,
- const spatialdata::geocoords::CoordSys* cs,
+pylith::faults::FaultCohesiveKin::initialize(const topology::Mesh& mesh,
const double_array& upDir,
const double_array& normalDir,
spatialdata::spatialdb::SpatialDB* matDB)
{ // initialize
assert(0 != _quadrature);
+
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+ assert(0 != cs);
if (3 != upDir.size())
throw std::runtime_error("Up direction for fault orientation must be "
@@ -90,6 +95,9 @@
//_faultMesh->getLabel("height")->view("Fault mesh height");
//_faultMesh->view("FAULT MESH");
+ delete _fields;
+ _fields = new topology::Fields<topology::Field<topology::SubMesh> >;
+
// Setup pseudo-stiffness of cohesive cells to improve conditioning
// of Jacobian matrix
_calcConditioning(cs, matDB);
@@ -106,28 +114,24 @@
++s_iter) {
EqKinSrc* src = s_iter->second;
assert(0 != src);
- src->initialize(_faultMesh, cs, *_normalizer);
+ src->initialize(_faultMesh, *_normalizer);
} // for
// Allocate slip field
- const ALE::Obj<SubMesh::label_sequence>& vertices = _faultMesh->depthStratum(0);
- _slip = new real_section_type(_faultMesh->comm(), _faultMesh->debug());
- _slip->setChart(real_section_type::chart_type(
- *std::min_element(vertices->begin(), vertices->end()),
- *std::max_element(vertices->begin(), vertices->end())+1));
- _slip->setFiberDimension(vertices, cs->spaceDim());
- _faultMesh->allocate(_slip);
- assert(!_slip.isNull());
+ const ALE::Obj<SubMesh::SieveMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ const ALE::Obj<SubMesh::label_sequence>& vertices =
+ _faultSieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+ _fields->add("slip");
+ topology::Field<topology::SubMesh>& slip = _fields->get("slip");
+ slip.newSection(vertices, cs->spaceDim());
+ slip.allocate();
// Allocate cumulative slip field
- _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));
- _cumSlip->setFiberDimension(vertices, cs->spaceDim());
- _faultMesh->allocate(_cumSlip);
- assert(!_cumSlip.isNull());
-
+ _fields->add("cumulative slip");
+ topology::Field<topology::SubMesh>& cumSlip = _fields->get("cumulative slip");
+ cumSlip.newSection(slip);
+ cumSlip.allocate();
} // initialize
// ----------------------------------------------------------------------
@@ -135,15 +139,11 @@
// require assembly across processors.
void
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)
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
{ // integrateResidual
- 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:
@@ -175,23 +175,26 @@
double_array cellAreaAssembled(numConstraintVert);
// Get cohesive cells
- const ALE::Obj<Mesh::label_sequence>& cellsCohesive =
- mesh->getLabelStratum("material-id", id());
+ const ALE::Obj<SieveMesh>& sieveMesh = residual.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cellsCohesive =
+ sieveMesh->getLabelStratum("material-id", id());
assert(!cellsCohesive.isNull());
- const Mesh::label_sequence::iterator cellsCohesiveBegin =
+ const SieveMesh::label_sequence::iterator cellsCohesiveBegin =
cellsCohesive->begin();
- const Mesh::label_sequence::iterator cellsCohesiveEnd =
+ const SieveMesh::label_sequence::iterator cellsCohesiveEnd =
cellsCohesive->end();
const int cellsCohesiveSize = cellsCohesive->size();
// Get section information
- const ALE::Obj<real_section_type>& solution = fields->getSolution();
- assert(!solution.isNull());
+ topology::Field<topology::Mesh>& solution = fields->solution();
+ const ALE::Obj<Mesh::RealSection>& solutionSection = solution.section();
+ assert(!solutionSection.isNull());
- for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
+ for (Mesh::SieveMesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
c_iter != cellsCohesiveEnd;
++c_iter) {
- const Mesh::point_type c_fault = _cohesiveToFault[*c_iter];
+ const Mesh::SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
cellArea = 0.0;
cellResidual = 0.0;
@@ -487,30 +490,32 @@
// ----------------------------------------------------------------------
// Update state variables as needed.
void
-pylith::faults::FaultCohesiveKin::updateState(const double t,
- topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh)
-{ // updateState
+pylith::faults::FaultCohesiveKin::updateStateVars(const double t,
+ topology::SolutionFields* const fields)
+{ // updateStateVars
assert(0 != fields);
- assert(!mesh.isNull());
- assert(!_faultMesh.isNull());
+ assert(0 != _fields);
// Update cumulative slip
- assert(!_cumSlip.isNull());
+ topology::Field<topology::SubMesh>& cumSlip = fields->get("cumulative slip");
+ topology::Field<topology::SubMesh>& slip = fields->get("slip");
if (!_useSolnIncr)
- _cumSlip->zero();
- _cumSlip->add(_cumSlip, _slip);
-} // updateState
+ cumSlip.zero();
+ cumSlip += slip;
+} // updateStateVars
// ----------------------------------------------------------------------
// Verify configuration is acceptable.
void
-pylith::faults::FaultCohesiveKin::verifyConfiguration(const ALE::Obj<Mesh>& mesh) const
+pylith::faults::FaultCohesiveKin::verifyConfiguration(
+ const topology::Mesh& mesh) const
{ // verifyConfiguration
- assert(!mesh.isNull());
assert(0 != _quadrature);
- if (!mesh->hasIntSection(label())) {
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ if (!sieveMesh->hasIntSection(label())) {
std::ostringstream msg;
msg << "Mesh missing group of vertices '" << label()
<< " for boundary condition.";
@@ -518,7 +523,7 @@
} // if
// check compatibility of mesh and quadrature scheme
- const int dimension = mesh->getDimension()-1;
+ const int dimension = mesh.dimension()-1;
if (_quadrature->cellDim() != dimension) {
std::ostringstream msg;
msg << "Dimension of reference cell in quadrature scheme ("
@@ -530,8 +535,8 @@
} // if
const int numCorners = _quadrature->refGeometry().numCorners();
- const ALE::Obj<Mesh::label_sequence>& cells =
- mesh->getLabelStratum("material-id", id());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", id());
assert(!cells.isNull());
const Mesh::label_sequence::iterator cellsBegin = cells->begin();
const Mesh::label_sequence::iterator cellsEnd = cells->end();
@@ -1094,44 +1099,5 @@
#endif
} // _calcTractionsChange
-// ----------------------------------------------------------------------
-// Allocate scalar field for output of vertex information.
-void
-pylith::faults::FaultCohesiveKin::_allocateBufferVertexScalar(void)
-{ // _allocateBufferVertexScalar
- const int fiberDim = 1;
- if (_bufferVertexScalar.isNull()) {
- _bufferVertexScalar = new real_section_type(_faultMesh->comm(),
- _faultMesh->debug());
- const ALE::Obj<SubMesh::label_sequence>& vertices =
- _faultMesh->depthStratum(0);
- _bufferVertexScalar->setChart(real_section_type::chart_type(
- *std::min_element(vertices->begin(), vertices->end()),
- *std::max_element(vertices->begin(), vertices->end())+1));
- _bufferVertexScalar->setFiberDimension(vertices, fiberDim);
- _faultMesh->allocate(_bufferVertexScalar);
- } // if
-} // _allocateBufferVertexScalar
-// ----------------------------------------------------------------------
-// Allocate vector field for output of vertex information.
-void
-pylith::faults::FaultCohesiveKin::_allocateBufferVertexVector(void)
-{ // _allocateBufferVertexVector
- assert(0 != _quadrature);
- const int fiberDim = _quadrature->spaceDim();
- if (_bufferVertexVector.isNull()) {
- _bufferVertexVector = new real_section_type(_faultMesh->comm(),
- _faultMesh->debug());
- const ALE::Obj<SubMesh::label_sequence>& vertices =
- _faultMesh->depthStratum(0);
- _bufferVertexVector->setChart(real_section_type::chart_type(
- *std::min_element(vertices->begin(), vertices->end()),
- *std::max_element(vertices->begin(), vertices->end())+1));
- _bufferVertexVector->setFiberDimension(vertices, fiberDim);
- _faultMesh->allocate(_bufferVertexVector);
- } // if
-} // _allocateBufferVertexVector
-
-
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh 2009-04-06 03:55:49 UTC (rev 14601)
@@ -207,12 +207,6 @@
void _calcTractionsChange(topology::Field<topology::SubMesh>* tractions,
const topology::Field<topology::Mesh>& solution);
- /// Allocate scalar field for output of vertex information.
- void _allocateBufferVertexScalar(void);
-
- /// Allocate vector field for output of vertex information.
- void _allocateBufferVertexVector(void);
-
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/LiuCosSlipFn.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -138,7 +138,7 @@
_dbRiseTime->open();
const char* riseTimeValues[] = {"rise-time"};
- _dbSlipTime->queryVals(slipTimeValues, 1);
+ _dbRiseTime->queryVals(riseTimeValues, 1);
// Get coordinates of vertices
const ALE::Obj<RealSection>& coordinates =
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/Makefile.am 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/Makefile.am 2009-04-06 03:55:49 UTC (rev 14601)
@@ -31,7 +31,8 @@
LiuCosSlipFn.icc \
SlipTimeFn.hh \
StepSlipFn.hh \
- StepSlipFn.icc
+ StepSlipFn.icc \
+ faultsfwd.hh
noinst_HEADERS = \
TopologyOps.hh \
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -13,28 +13,29 @@
#include <portinfo>
#include "TopologyOps.hh" // implementation of object methods
-//#include <Selection.hh> // Algorithms for submeshes
-//#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+#include "TopologyVisitors.hh" // USES ClassifyVisitor
+#include <Selection.hh> // Algorithms for submeshes
+
#include <cassert> // USES assert()
// ----------------------------------------------------------------------
template<class InputPoints>
bool
-pylith::faults::CohesiveTopology::compatibleOrientation(const ALE::Obj<Mesh>& mesh,
- const Mesh::point_type& p,
- const Mesh::point_type& q,
+pylith::faults::TopologyOps::compatibleOrientation(const ALE::Obj<SieveMesh>& mesh,
+ const SieveMesh::point_type& p,
+ const SieveMesh::point_type& q,
const int numFaultCorners,
const int faultFaceSize,
const int faultDepth,
- const Obj<InputPoints>& points,
+ const ALE::Obj<InputPoints>& points,
int indices[],
PointArray *origVertices,
PointArray *faceVertices,
PointArray *neighborVertices)
{
- typedef ALE::Selection<Mesh> selection;
+ typedef ALE::Selection<SieveMesh> selection;
const int debug = mesh->debug();
bool compatible;
@@ -59,16 +60,16 @@
// ----------------------------------------------------------------------
void
-pylith::faults::CohesiveTopology::computeCensoredDepth(const ALE::Obj<Mesh::label_type>& depth,
- const ALE::Obj<Mesh::sieve_type>& sieve,
- const Mesh::point_type& firstCohesiveCell)
+pylith::faults::TopologyOps::computeCensoredDepth(const ALE::Obj<SieveMesh::label_type>& depth,
+ const ALE::Obj<SieveMesh::sieve_type>& sieve,
+ const SieveMesh::point_type& firstCohesiveCell)
{
- Mesh::DepthVisitor d(*sieve, firstCohesiveCell, *depth);
+ SieveMesh::DepthVisitor d(*sieve, firstCohesiveCell, *depth);
sieve->roots(d);
while(d.isModified()) {
// FIX: Avoid the copy here somehow by fixing the traversal
- std::vector<Mesh::point_type> modifiedPoints(d.getModifiedPoints().begin(), d.getModifiedPoints().end());
+ std::vector<SieveMesh::point_type> modifiedPoints(d.getModifiedPoints().begin(), d.getModifiedPoints().end());
d.clear();
sieve->support(modifiedPoints, d);
@@ -77,17 +78,18 @@
// ----------------------------------------------------------------------
void
-pylith::faults::CohesiveTopology::classifyCells(const ALE::Obj<Mesh::sieve_type>& sieve,
- const Mesh::point_type& vertex,
+pylith::faults::TopologyOps::classifyCells(const ALE::Obj<SieveMesh::sieve_type>& sieve,
+ const SieveMesh::point_type& vertex,
const int depth,
const int faceSize,
- const Mesh::point_type& firstCohesiveCell,
+ const SieveMesh::point_type& firstCohesiveCell,
PointSet& replaceCells,
PointSet& noReplaceCells,
const int debug)
{
// Replace all cells on a given side of the fault with a vertex on the fault
- ClassifyVisitor<Mesh::sieve_type> cV(*sieve, replaceCells, noReplaceCells, firstCohesiveCell, faceSize, debug);
+ ClassifyVisitor<SieveMesh::sieve_type> cV(*sieve, replaceCells, noReplaceCells,
+ firstCohesiveCell, faceSize, debug);
const PointSet& vReplaceCells = cV.getReplaceCells();
const PointSet& vNoReplaceCells = cV.getNoReplaceCells();
@@ -117,16 +119,16 @@
// ----------------------------------------------------------------------
void
-pylith::faults::CohesiveTopology::createFaultSieveFromVertices(const int dim,
+pylith::faults::TopologyOps::createFaultSieveFromVertices(const int dim,
const int firstCell,
const PointSet& faultVertices,
- const Obj<Mesh>& mesh,
- const Obj<ALE::Mesh::arrow_section_type>& orientation,
- const Obj<ALE::Mesh::sieve_type>& faultSieve,
+ const ALE::Obj<SieveMesh>& mesh,
+ const ALE::Obj<ALE::Mesh::arrow_section_type>& orientation,
+ const ALE::Obj<ALE::Mesh::sieve_type>& faultSieve,
const bool flipFault)
{
typedef ALE::Selection<ALE::Mesh> selection;
- const Obj<sieve_type>& sieve = mesh->getSieve();
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = mesh->getSieve();
const PointSet::const_iterator fvBegin = faultVertices.begin();
const PointSet::const_iterator fvEnd = faultVertices.end();
int curCell = firstCell;
@@ -135,7 +137,7 @@
int o = 1;
ALE::Mesh::point_type f = firstCell;
const int debug = mesh->debug();
- Obj<PointSet> face = new PointSet();
+ ALE::Obj<PointSet> face = new PointSet();
int numCorners = 0; // The number of vertices in a mesh cell
int faceSize = 0; // The number of vertices in a mesh face
int *indices = NULL; // The indices of a face vertex set in a cell
@@ -161,18 +163,18 @@
// This only works for uninterpolated meshes
assert((mesh->depth() == 1) || (mesh->depth() == -1));
- ALE::ISieveVisitor::PointRetriever<sieve_type> sV(std::max(1, sieve->getMaxSupportSize()));
- ALE::ISieveVisitor::PointRetriever<sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
+ ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> sV(std::max(1, sieve->getMaxSupportSize()));
+ ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
for(PointSet::const_iterator fv_iter = fvBegin; fv_iter != fvEnd; ++fv_iter) {
sieve->support(*fv_iter, sV);
- const Mesh::point_type *support = sV.getPoints();
+ const SieveMesh::point_type *support = sV.getPoints();
if (debug) std::cout << "Checking fault vertex " << *fv_iter << std::endl;
const int sVsize = sV.getSize();
for (int i=0; i < sVsize; ++i) {
const int s = (!flipFault) ? i : sVsize - i - 1;
sieve->cone(support[s], cV);
- const Mesh::point_type *cone = cV.getPoints();
+ const SieveMesh::point_type *cone = cV.getPoints();
if (debug) std::cout << " Checking cell " << support[s] << std::endl;
if (faultCells.find(support[s]) != faultCells.end()) {
@@ -191,7 +193,7 @@
"element on the fault");
if (face->size() == faceSize) {
if (debug) std::cout << " Contains a face on the fault" << std::endl;
- ALE::Obj<sieve_type::supportSet> preFace;
+ ALE::Obj<SieveMesh::sieve_type::supportSet> preFace;
if (dim < 2) {
preFace = faultSieve->nJoin1(face);
} else {
@@ -240,14 +242,14 @@
// ----------------------------------------------------------------------
void
-pylith::faults::CohesiveTopology::createFaultSieveFromFaces(const int dim,
+pylith::faults::TopologyOps::createFaultSieveFromFaces(const int dim,
const int firstCell,
const int numFaces,
const int faultVertices[],
const int faultCells[],
- const Obj<Mesh>& mesh,
- const Obj<ALE::Mesh::arrow_section_type>& orientation,
- const Obj<ALE::Mesh::sieve_type>& faultSieve)
+ const ALE::Obj<SieveMesh>& mesh,
+ const ALE::Obj<ALE::Mesh::arrow_section_type>& orientation,
+ const ALE::Obj<ALE::Mesh::sieve_type>& faultSieve)
{
typedef ALE::Selection<ALE::Mesh> selection;
int faceSize = 0; // The number of vertices in a mesh face
@@ -295,24 +297,24 @@
// ----------------------------------------------------------------------
void
-pylith::faults::CohesiveTopology::orientFaultSieve(const int dim,
- const Obj<Mesh>& mesh,
- const Obj<ALE::Mesh::arrow_section_type>& orientation,
- const Obj<ALE::Mesh>& fault)
+pylith::faults::TopologyOps::orientFaultSieve(const int dim,
+ const ALE::Obj<SieveMesh>& mesh,
+ const ALE::Obj<ALE::Mesh::arrow_section_type>& orientation,
+ const ALE::Obj<ALE::Mesh>& fault)
{
// Must check the orientation here
typedef ALE::Selection<ALE::Mesh> selection;
- const Obj<ALE::Mesh::sieve_type>& faultSieve = fault->getSieve();
- const Mesh::point_type firstFaultCell = *fault->heightStratum(1)->begin();
- const Obj<ALE::Mesh::label_sequence>& fFaces = fault->heightStratum(2);
+ const ALE::Obj<ALE::Mesh::sieve_type>& faultSieve = fault->getSieve();
+ const SieveMesh::point_type firstFaultCell = *fault->heightStratum(1)->begin();
+ const ALE::Obj<ALE::Mesh::label_sequence>& fFaces = fault->heightStratum(2);
const int numFaultFaces = fFaces->size();
const int faultDepth = fault->depth()-1; // Depth of fault cells
int numFaultCorners = 0; // The number of vertices in a fault cell
int faultFaceSize = 0; // The number of vertices in a face between fault cells
int faceSize = 0; // The number of vertices in a mesh face
const int debug = fault->debug();
- Obj<PointSet> newCells = new PointSet();
- Obj<PointSet> loopCells = new PointSet();
+ ALE::Obj<PointSet> newCells = new PointSet();
+ ALE::Obj<PointSet> loopCells = new PointSet();
PointSet flippedCells; // Incorrectly oriented fault cells
PointSet facesSeen; // Fault faces already considered
PointSet cellsSeen; // Fault cells already matched
@@ -337,14 +339,14 @@
newCells->insert(firstFaultCell);
while(facesSeen.size() != numFaultFaces) {
- Obj<PointSet> tmp = newCells; newCells = loopCells; loopCells = tmp;
+ ALE::Obj<PointSet> tmp = newCells; newCells = loopCells; loopCells = tmp;
newCells->clear();
if (!loopCells->size()) {throw ALE::Exception("Fault surface not a single connected component.");}
// Loop over new cells
for(PointSet::iterator c_iter = loopCells->begin(); c_iter != loopCells->end(); ++c_iter) {
// Loop over edges of this cell
- const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& cone = faultSieve->cone(*c_iter);
+ const ALE::Obj<ALE::Mesh::sieve_type::traits::coneSequence>& cone = faultSieve->cone(*c_iter);
const ALE::Mesh::sieve_type::traits::coneSequence::iterator eBegin = cone->begin();
const ALE::Mesh::sieve_type::traits::coneSequence::iterator eEnd = cone->end();
@@ -352,7 +354,7 @@
if (facesSeen.find(*e_iter) != facesSeen.end()) continue;
facesSeen.insert(*e_iter);
if (debug) std::cout << " Checking orientation of fault face " << *e_iter << std::endl;
- const Obj<ALE::Mesh::sieve_type::traits::supportSequence>& support = faultSieve->support(*e_iter);
+ const ALE::Obj<ALE::Mesh::sieve_type::traits::supportSequence>& support = faultSieve->support(*e_iter);
ALE::Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin();
// Throw out boundary fault faces
@@ -369,9 +371,9 @@
if (debug) std::cout << " neighboring cells " << cellA << " and " << cellB << std::endl;
// In 1D, just check that vertices match
if (dim == 1) {
- const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneA = faultSieve->cone(cellA);
+ const ALE::Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneA = faultSieve->cone(cellA);
ALE::Mesh::sieve_type::traits::coneSequence::iterator iterA = coneA->begin();
- const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneB = faultSieve->cone(cellB);
+ const ALE::Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneB = faultSieve->cone(cellB);
ALE::Mesh::sieve_type::traits::coneSequence::iterator iterB = coneB->begin();
int posA, posB;
@@ -462,7 +464,7 @@
for(ALE::Mesh::label_sequence::iterator e_iter = fFaces->begin(); e_iter != fFaces->end(); ++e_iter) {
if (debug) std::cout << " Checking orientation of fault face " << *e_iter << std::endl;
// for each face get the support (2 fault cells)
- const Obj<ALE::Mesh::sieve_type::traits::supportSequence>& support = faultSieve->support(*e_iter);
+ const ALE::Obj<ALE::Mesh::sieve_type::traits::supportSequence>& support = faultSieve->support(*e_iter);
ALE::Mesh::sieve_type::traits::supportSequence::iterator s_iter = support->begin();
// Throw out boundary fault faces
@@ -473,9 +475,9 @@
if (debug) std::cout << " neighboring cells " << cellA << " and " << cellB << std::endl;
// In 1D, just check that vertices match
if (dim == 1) {
- const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneA = faultSieve->cone(cellA);
+ const ALE::Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneA = faultSieve->cone(cellA);
ALE::Mesh::sieve_type::traits::coneSequence::iterator iterA = coneA->begin();
- const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneB = faultSieve->cone(cellB);
+ const ALE::Obj<ALE::Mesh::sieve_type::traits::coneSequence>& coneB = faultSieve->cone(cellB);
ALE::Mesh::sieve_type::traits::coneSequence::iterator iterB = coneB->begin();
int posA, posB;
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.hh 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyOps.hh 2009-04-06 03:55:49 UTC (rev 14601)
@@ -21,13 +21,18 @@
// Include directives ---------------------------------------------------
#include "faultsfwd.hh" // forward declarations
+#include "pylith/topology/Mesh.hh" // USES Mesh
+
// TopologyOps ----------------------------------------------------------
class pylith::faults::TopologyOps
{ // class TopologyOps
+
+ // PUBLIC TYPEDEFS ////////////////////////////////////////////////////
public :
+ typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef std::set<SieveMesh::point_type> PointSet;
- typedef std::vector<sieve_type::point_type> PointArray;
- typedef std::pair<sieve_type::point_type, int> oPoint_type;
+ typedef std::vector<SieveMesh::sieve_type::point_type> PointArray;
+ typedef std::pair<SieveMesh::sieve_type::point_type, int> oPoint_type;
typedef std::vector<oPoint_type> oPointArray;
// PUBLIC METHODS /////////////////////////////////////////////////////
@@ -35,7 +40,7 @@
template<class InputPoints>
static
- bool compatibleOrientation(const topology::& mesh,
+ bool compatibleOrientation(const ALE::Obj<SieveMesh>& mesh,
const SieveMesh::point_type& p,
const SieveMesh::point_type& q,
const int numFaultCorners,
@@ -48,16 +53,16 @@
PointArray *neighborVertices);
static
- void computeCensoredDepth(const ALE::Obj<Mesh::label_type>& depth,
- const ALE::Obj<Mesh::sieve_type>& sieve,
- const Mesh::point_type& firstCohesiveCell);
+ void computeCensoredDepth(const ALE::Obj<SieveMesh::label_type>& depth,
+ const ALE::Obj<SieveMesh::sieve_type>& sieve,
+ const SieveMesh::point_type& firstCohesiveCell);
static
- void classifyCells(const ALE::Obj<Mesh::sieve_type>& sieve,
- const Mesh::point_type& vertex,
+ void classifyCells(const ALE::Obj<SieveMesh::sieve_type>& sieve,
+ const SieveMesh::point_type& vertex,
const int depth,
const int faceSize,
- const Mesh::point_type& firstCohesiveCell,
+ const SieveMesh::point_type& firstCohesiveCell,
PointSet& replaceCells,
PointSet& noReplaceCells,
const int debug);
@@ -66,7 +71,7 @@
void createFaultSieveFromVertices(const int dim,
const int firstCell,
const PointSet& faultVertices,
- const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<SieveMesh>& mesh,
const ALE::Obj<ALE::Mesh::arrow_section_type>& orientation,
const ALE::Obj<ALE::Mesh::sieve_type>& faultSieve,
const bool flipFault);
@@ -77,13 +82,13 @@
const int numFaces,
const int faultVertices[],
const int faultCells[],
- const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<SieveMesh>& mesh,
const ALE::Obj<ALE::Mesh::arrow_section_type>& orientation,
const ALE::Obj<ALE::Mesh::sieve_type>& faultSieve);
static
void orientFaultSieve(const int dim,
- const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<SieveMesh>& mesh,
const ALE::Obj<ALE::Mesh::arrow_section_type>& orientation,
const ALE::Obj<ALE::Mesh>& fault);
}; // class CohesiveTopology
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -17,7 +17,7 @@
pylith::faults::ReplaceVisitor<Sieve,Renumbering>::ReplaceVisitor(
Renumbering& r,
const int size,
- const int debug = 0) :
+ const int debug) :
renumbering(r),
size(size),
i(0),
@@ -65,7 +65,7 @@
// ----------------------------------------------------------------------
template<typename Sieve, typename Renumbering>
inline
-const point_type*
+const typename Sieve::point_type*
pylith::faults::ReplaceVisitor<Sieve,Renumbering>::getPoints(void)
{ // getPoints
return this->points;
@@ -97,7 +97,7 @@
const PointSet& nrC,
const point_type& fC,
const int fS,
- const int debug =0) :
+ const int debug) :
sieve(s),
replaceCells(rC),
noReplaceCells(nrC),
@@ -152,7 +152,7 @@
return;
} // if
// If neighbor shares a face with anyone in replaceCells, then add
- for(PointSet::const_iterator c_iter = vReplaceCells.begin();
+ for (typename PointSet::const_iterator c_iter = vReplaceCells.begin();
c_iter != vReplaceCells.end();
++c_iter) {
sieve.meet(*c_iter, point, pR);
@@ -172,7 +172,7 @@
if (classified)
return;
// It is unclear whether taking out the noReplace cells will speed this up
- for(PointSet::const_iterator c_iter = vNoReplaceCells.begin();
+ for (typename PointSet::const_iterator c_iter = vNoReplaceCells.begin();
c_iter != vNoReplaceCells.end();
++c_iter) {
sieve.meet(*c_iter, point, pR);
@@ -202,7 +202,7 @@
// ----------------------------------------------------------------------
template<typename Sieve>
inline
-const PointSet&
+const std::set<typename Sieve::point_type>&
pylith::faults::ClassifyVisitor<Sieve>::getReplaceCells(void) const
{ // getReplaceCells
return this->vReplaceCells;
@@ -211,7 +211,7 @@
// ----------------------------------------------------------------------
template<typename Sieve>
inline
-const PointSet&
+const std::set<typename Sieve::point_type>&
pylith::faults::ClassifyVisitor<Sieve>::getNoReplaceCells() const
{ // getNoReplaceCells
return this->vNoReplaceCells;
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.hh 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/TopologyVisitors.hh 2009-04-06 03:55:49 UTC (rev 14601)
@@ -24,7 +24,7 @@
// ReplaceVisitor -------------------------------------------------------
template<typename Sieve, typename Renumbering>
class pylith::faults::ReplaceVisitor {
-public:
+private:
typedef typename Sieve::point_type point_type;
protected:
Renumbering& renumbering;
@@ -50,6 +50,7 @@
class pylith::faults::ClassifyVisitor {
public:
typedef typename Sieve::point_type point_type;
+ typedef std::set<point_type> PointSet;
protected:
const Sieve& sieve;
const PointSet& replaceCells;
@@ -75,12 +76,13 @@
void visitArrow(const typename Sieve::arrow_type&);
const PointSet& getReplaceCells() const;
const PointSet& getNoReplaceCells() const;
- const bool getModified() const;
- const int getSize() const;
- void setMode(const bool isSetup);
- void reset();
+ bool getModified() const;
+ int getSize() const;
+ void setMode(const bool isSetup);
+ void reset();
};
+#include "TopologyVisitors.cc" // template definitions
#endif // pylith_faults_topologyvisitors_hh
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/meshio/UCDFaultFile.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/meshio/UCDFaultFile.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/meshio/UCDFaultFile.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -16,6 +16,7 @@
#include "MeshBuilder.hh" // USES MeshBuilder
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/utils/array.hh" // USES double_array, int_array
#include <petsc.h> // USES MPI_Comm
@@ -26,10 +27,12 @@
// ----------------------------------------------------------------------
void
pylith::meshio::UCDFaultFile::read(const char* filename,
- const ALE::Obj<SieveMesh>& mesh,
- const ALE::Obj<SieveMesh>& fault,
- ALE::Obj<ALE::Mesh>& faultBd)
+ topology::SubMesh* faultMesh,
+ ALE::Obj<ALE::Mesh>& faultBoundary,
+ const topology::Mesh& mesh)
{ // read
+ assert(0 != faultMesh);
+
int faultDim = 2;
int fSpaceDim = 0;
int numFVertices = 0;
@@ -40,16 +43,23 @@
int_array fMaterialIds;
int_array faceCells;
int_array vertexIDs;
-
- if (!fault->commRank()) {
+
+ const ALE::Obj<topology::Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ int rank = 0;
+ MPI_Comm_rank(mesh.comm(), &rank);
+ if (0 == rank) {
std::ifstream fin(filename, std::ios::in);
if (!(fin.is_open() && fin.good())) {
std::ostringstream msg;
msg << "Could not open ASCII INP file '" << filename << "' for reading.";
throw std::runtime_error(msg.str());
} // if
- int numNodeData, numCellData, numModelData;
-
+ int numNodeData = 0;
+ int numCellData = 0;
+ int numModelData = 0;
+
fSpaceDim = 3;
// Section 1: <num_nodes> <num_cells> <num_ndata> <num_cdata> <num_mdata>
@@ -91,7 +101,8 @@
} // for
// Section 4: <num_comp for node data> <size comp 1> <size comp 2>...<size comp n>
- int numComponents, totalSize;
+ int numComponents = 0;
+ int totalSize = 0;
fin >> numComponents;
totalSize = 0;
@@ -153,7 +164,8 @@
// Only do this once since we add cohesive cells after that
static int numCells = -1;
- if (numCells == -1) {numCells = mesh->heightStratum(0)->size();}
+ assert(!sieveMesh->heightStratum(0).isNull());
+ if (numCells == -1) {numCells = sieveMesh->heightStratum(0)->size();}
// Renumber vertices and use zero based indices
// UCD file has one-based indices for both vertexIDs and fCells
@@ -162,16 +174,17 @@
for(int corner = 0; corner < numFCorners; ++corner)
fCells[c*numFCorners+corner] =
vertexIDs[fCells[c*numFCorners+corner]-1] - 1 + numCells;
-
+
// Switch to zero based index for global cell numbering
for (int c=0; c < numFCells; ++c)
for (int i=0; i < 2; ++i)
faceCells[c*2+i] -= 1;
} // if
-
+
+ assert(!sieveMesh->getSieve().isNull());
const int firstFaultCell =
- mesh->getSieve()->getBaseSize() + mesh->getSieve()->getCapSize();
- MeshBuilder::buildFaultMesh(fault, faultBd,
+ sieveMesh->getSieve()->getBaseSize() + sieveMesh->getSieve()->getCapSize();
+ MeshBuilder::buildFaultMesh(faultMesh->sieveMesh(), faultBoundary,
fCoordinates, numFVertices, fSpaceDim, fCells,
numFCells, numFCorners, firstFaultCell,
faceCells, faultDim);
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/meshio/UCDFaultFile.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/meshio/UCDFaultFile.hh 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/meshio/UCDFaultFile.hh 2009-04-06 03:55:49 UTC (rev 14601)
@@ -26,19 +26,13 @@
#define pylith_meshio_ucdfaultfile_hh
// Include directives ---------------------------------------------------
-#define NEWPYLITHMESH 1
-#include "pylith/utils/sievetypes.hh" // USES SieveMesh, ALE::Mesh
+#include "meshiofwd.hh" // forward declarations
-// Forward declarations -------------------------------------------------
-namespace pylith {
- namespace meshio {
- class UCDFaultFile;
+#include "pylith/topology/topologyfwd.hh" // USES Mesh, SubMesh
- class TestFaultFile; // unit testing
- } // meshio
-} // pylith
+#include "pylith/utils/sievetypes.hh" // USES ALE::Obj, ALE::Mesh
-// MeshIOLagrit ---------------------------------------------------------
+// UCDFaultFile ---------------------------------------------------------
class pylith::meshio::UCDFaultFile
{ // UCDFaultFile
friend class TestUCDFaultFile; // unit testing
@@ -48,9 +42,9 @@
static
void read(const char* filename,
- const ALE::Obj<SieveMesh>& mesh,
- const ALE::Obj<SieveMesh>& fault,
- ALE::Obj<ALE::Mesh>& faultBd);
+ topology::SubMesh* faultMesh,
+ ALE::Obj<ALE::Mesh>& faultBoundary,
+ const topology::Mesh& mesh);
}; // UCDFaultFile
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -90,11 +90,19 @@
// Set data from mesh.
_mesh->setDebug(mesh.debug());
+ coordsys(mesh);
+} // createSubMesh
+
+// ----------------------------------------------------------------------
+// Set coordinate system using mesh.
+void
+pylith::topology::SubMesh::coordsys(const Mesh& mesh)
+{ // coordsys
delete _coordsys; _coordsys = 0;
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
if (0 != cs)
_coordsys = cs->clone();
-} // createSubMesh
+} // coordsys
// ----------------------------------------------------------------------
// Initialize the finite-element mesh.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh 2009-04-06 03:55:49 UTC (rev 14601)
@@ -70,7 +70,7 @@
* @param label Label for vertices marking boundary.
*/
void createSubMesh(const Mesh& mesh,
- const char* label);
+ const char* label);
/** Get Sieve mesh.
*
@@ -84,6 +84,12 @@
*/
ALE::Obj<SieveMesh>& sieveMesh(void);
+ /** Set coordinate system using mesh.
+ *
+ * @param mesh Finite-element mesh over domain.
+ */
+ void coordsys(const Mesh& mesh);
+
/** Get coordinate system.
*
* @returns Coordinate system.
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/Makefile.am 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/Makefile.am 2009-04-06 03:55:49 UTC (rev 14601)
@@ -12,12 +12,12 @@
SUBDIRS = \
bc \
+ faults \
feassemble \
materials \
meshio \
topology \
utils
-# faults
# End of file
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/Makefile.am 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/Makefile.am 2009-04-06 03:55:49 UTC (rev 14601)
@@ -21,33 +21,34 @@
# Primary source files
testfaults_SOURCES = \
- TestBruneSlipFn.cc \
- TestLiuCosSlipFn.cc \
- TestConstRateSlipFn.cc \
- TestStepSlipFn.cc \
+ TestFaultCohesive.cc \
TestEqKinSrc.cc \
- TestFault.cc \
- TestFaultCohesive.cc \
- TestFaultCohesiveKin.cc \
- TestFaultCohesiveKinLine2.cc \
- TestFaultCohesiveKinTri3.cc \
- TestFaultCohesiveKinTri3d.cc \
- TestFaultCohesiveKinQuad4.cc \
- TestFaultCohesiveKinQuad4e.cc \
- TestFaultCohesiveKinTet4.cc \
- TestFaultCohesiveKinTet4e.cc \
- TestFaultCohesiveKinTet4f.cc \
- TestFaultCohesiveKinHex8.cc \
- TestFaultCohesiveKinSrcs.cc \
- TestFaultCohesiveKinSrcsLine2.cc \
- TestFaultCohesiveKinSrcsTri3.cc \
- TestFaultCohesiveKinSrcsQuad4.cc \
- TestFaultCohesiveKinSrcsTet4.cc \
- TestFaultCohesiveKinSrcsHex8.cc \
test_faults.cc
+# TestBruneSlipFn.cc \
+# TestLiuCosSlipFn.cc \
+# TestConstRateSlipFn.cc \
+# TestStepSlipFn.cc \
+# TestFault.cc \
+# TestFaultCohesiveKin.cc \
+# TestFaultCohesiveKinLine2.cc \
+# TestFaultCohesiveKinTri3.cc \
+# TestFaultCohesiveKinTri3d.cc \
+# TestFaultCohesiveKinQuad4.cc \
+# TestFaultCohesiveKinQuad4e.cc \
+# TestFaultCohesiveKinTet4.cc \
+# TestFaultCohesiveKinTet4e.cc \
+# TestFaultCohesiveKinTet4f.cc \
+# TestFaultCohesiveKinHex8.cc \
+# TestFaultCohesiveKinSrcs.cc \
+# TestFaultCohesiveKinSrcsLine2.cc \
+# TestFaultCohesiveKinSrcsTri3.cc \
+# TestFaultCohesiveKinSrcsQuad4.cc \
+# TestFaultCohesiveKinSrcsTet4.cc \
+# TestFaultCohesiveKinSrcsHex8.cc
+
noinst_HEADERS = \
TestBruneSlipFn.hh \
TestLiuCosSlipFn.hh \
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestEqKinSrc.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestEqKinSrc.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestEqKinSrc.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -16,10 +16,13 @@
#include "pylith/faults/EqKinSrc.hh" // USES EqKinSrc
+
#include "pylith/faults/BruneSlipFn.hh" // USES BruneSlipFn
#include "pylith/faults/CohesiveTopology.hh" // USES CohesiveTopology
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
-#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
@@ -30,6 +33,11 @@
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestEqKinSrc );
// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::SieveSubMesh SieveSubMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
// Test constructor.
void
pylith::faults::TestEqKinSrc::testConstructor(void)
@@ -55,11 +63,12 @@
void
pylith::faults::TestEqKinSrc::testInitialize(void)
{ // testInitialize
- ALE::Obj<Mesh> faultMesh;
+ topology::Mesh mesh;
+ topology::SubMesh faultMesh;
EqKinSrc eqsrc;
BruneSlipFn slipfn;
const double originTime = 2.45;
- _initialize(&faultMesh, &eqsrc, &slipfn, originTime);
+ _initialize(&mesh, &faultMesh, &eqsrc, &slipfn, originTime);
// Don't have access to details of slip time function, so we can't
// check parameters. Have to rely on test of slip() for verification
@@ -74,46 +83,51 @@
const double finalSlipE[] = { 2.3, 0.1,
2.4, 0.2};
const double slipTimeE[] = { 1.2, 1.3 };
- const double peakRateE[] = { 1.4, 1.5 };
+ const double riseTimeE[] = { 1.4, 1.5 };
const double originTime = 2.42;
- ALE::Obj<Mesh> faultMesh;
+ topology::Mesh mesh;
+ topology::SubMesh faultMesh;
EqKinSrc eqsrc;
BruneSlipFn slipfn;
- _initialize(&faultMesh, &eqsrc, &slipfn, originTime);
+ _initialize(&mesh, &faultMesh, &eqsrc, &slipfn, originTime);
- const int spaceDim = faultMesh->getDimension() + 1;
- const ALE::Obj<Mesh::label_sequence>& vertices = faultMesh->depthStratum(0);
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- ALE::Obj<real_section_type> slip =
- new real_section_type(faultMesh->comm(), faultMesh->debug());
- slip->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(),
- vertices->end()),
- *std::max_element(vertices->begin(), vertices->end())+1));
- slip->setFiberDimension(vertices, spaceDim);
- faultMesh->allocate(slip);
- CPPUNIT_ASSERT(!slip.isNull());
+ const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
+ CPPUNIT_ASSERT(0 != cs);
+ const int spaceDim = cs->spaceDim();
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
+ CPPUNIT_ASSERT(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ faultSieveMesh->depthStratum(0);
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ topology::Field<topology::SubMesh> slip(faultMesh);
+ slip.newSection(vertices, spaceDim);
+ slip.allocate();
+
const double t = 2.134;
- eqsrc.slip(slip, originTime+t, faultMesh);
+ eqsrc.slip(&slip, originTime+t);
const double tolerance = 1.0e-06;
int iPoint = 0;
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+
+ const ALE::Obj<RealSection>& slipSection = slip.section();
+ CPPUNIT_ASSERT(!slipSection.isNull());
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter, ++iPoint) {
double slipMag = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim)
slipMag += pow(finalSlipE[iPoint*spaceDim+iDim], 2);
slipMag = sqrt(slipMag);
- const double tau = slipMag / (exp(1.0) * peakRateE[iPoint]);
+ const double peakRate = slipMag / riseTimeE[iPoint] * 1.745;
+ const double tau = slipMag / (exp(1.0) * peakRate);
const double t0 = slipTimeE[iPoint];
const double slipNorm = 1.0 - exp(-(t-t0)/tau) * (1.0 + (t-t0)/tau);
- const int fiberDim = slip->getFiberDimension(*v_iter);
+ const int fiberDim = slipSection->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const real_section_type::value_type* vals =
- slip->restrictPoint(*v_iter);
+ const double* vals = slipSection->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != vals);
for (int iDim=0; iDim < fiberDim; ++iDim) {
@@ -131,50 +145,54 @@
const double finalSlipE[] = { 2.3, 0.1,
2.4, 0.2};
const double slipTimeE[] = { 1.2, 1.3 };
- const double peakRateE[] = { 1.4, 1.5 };
+ const double riseTimeE[] = { 1.4, 1.5 };
const double originTime = -4.29;
- ALE::Obj<Mesh> faultMesh;
+ topology::Mesh mesh;
+ topology::SubMesh faultMesh;
EqKinSrc eqsrc;
BruneSlipFn slipfn;
- _initialize(&faultMesh, &eqsrc, &slipfn, originTime);
+ _initialize(&mesh, &faultMesh, &eqsrc, &slipfn, originTime);
- const int spaceDim = faultMesh->getDimension() + 1;
- const ALE::Obj<Mesh::label_sequence>& vertices = faultMesh->depthStratum(0);
- const Mesh::label_sequence::iterator verticesEnd = vertices->end();
- ALE::Obj<real_section_type> slip =
- new real_section_type(faultMesh->comm(), faultMesh->debug());
- slip->setChart(real_section_type::chart_type(*std::min_element(vertices->begin(),
- vertices->end()),
- *std::max_element(vertices->begin(), vertices->end())+1));
- slip->setFiberDimension(vertices, spaceDim);
- faultMesh->allocate(slip);
- CPPUNIT_ASSERT(!slip.isNull());
+ const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
+ CPPUNIT_ASSERT(0 != cs);
+ const int spaceDim = cs->spaceDim();
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
+ CPPUNIT_ASSERT(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ faultSieveMesh->depthStratum(0);
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ topology::Field<topology::SubMesh> slip(faultMesh);
+ slip.newSection(vertices, spaceDim);
+ slip.allocate();
+
const double t0 = 1.234;
const double t1 = 2.525;
- eqsrc.slipIncr(slip, originTime+t0, originTime+t1, faultMesh);
+ eqsrc.slipIncr(&slip, originTime+t0, originTime+t1);
const double tolerance = 1.0e-06;
int iPoint = 0;
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
+ const ALE::Obj<RealSection>& slipSection = slip.section();
+ CPPUNIT_ASSERT(!slipSection.isNull());
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
v_iter != verticesEnd;
++v_iter, ++iPoint) {
double slipMag = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim)
slipMag += pow(finalSlipE[iPoint*spaceDim+iDim], 2);
slipMag = sqrt(slipMag);
- const double tau = slipMag / (exp(1.0) * peakRateE[iPoint]);
+ const double peakRate = slipMag / riseTimeE[iPoint] * 1.745;
+ const double tau = slipMag / (exp(1.0) * peakRate);
const double tRef = slipTimeE[iPoint];
const double slipNorm0 =
(t0 > tRef) ? 1.0 - exp(-(t0-tRef)/tau) * (1.0 + (t0-tRef)/tau) : 0.0;
const double slipNorm1 =
(t1 > tRef) ? 1.0 - exp(-(t1-tRef)/tau) * (1.0 + (t1-tRef)/tau) : 0.0;
- const int fiberDim = slip->getFiberDimension(*v_iter);
+ const int fiberDim = slipSection->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const real_section_type::value_type* vals =
- slip->restrictPoint(*v_iter);
+ const double* vals = slipSection->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != vals);
for (int iDim=0; iDim < fiberDim; ++iDim) {
@@ -188,45 +206,52 @@
// ----------------------------------------------------------------------
// Initialize EqKinSrc.
void
-pylith::faults::TestEqKinSrc::_initialize(ALE::Obj<Mesh>* faultMesh,
+pylith::faults::TestEqKinSrc::_initialize(topology::Mesh* mesh,
+ topology::SubMesh* faultMesh,
EqKinSrc* eqsrc,
BruneSlipFn* slipfn,
const double originTime)
{ // _initialize
+ CPPUNIT_ASSERT(0 != faultMesh);
+ CPPUNIT_ASSERT(0 != eqsrc);
+ CPPUNIT_ASSERT(0 != slipfn);
+
const char* meshFilename = "data/tri3.mesh";
const char* faultLabel = "fault";
const int faultId = 2;
const char* finalSlipFilename = "data/tri3_finalslip.spatialdb";
const char* slipTimeFilename = "data/tri3_sliptime.spatialdb";
- const char* peakRateFilename = "data/tri3_peakrate.spatialdb";
+ const char* peakRateFilename = "data/tri3_risetime.spatialdb";
- ALE::Obj<Mesh> mesh;
meshio::MeshIOAscii meshIO;
meshIO.filename(meshFilename);
meshIO.debug(false);
meshIO.interpolate(false);
- meshIO.read(&mesh);
- CPPUNIT_ASSERT(!mesh.isNull());
- const int spaceDim = mesh->getDimension();
+ meshIO.read(mesh);
+
+ // Set up coordinates
spatialdata::geocoords::CSCart cs;
- cs.setSpaceDim(spaceDim);
+ cs.setSpaceDim(mesh->dimension());
+ cs.initialize();
+ mesh->coordsys(&cs);
// Create fault mesh
const bool useLagrangeConstraints = true;
- (*faultMesh) = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
- ALE::Obj<ALE::Mesh> faultBd = NULL;
- CohesiveTopology::createFault(*faultMesh, faultBd,
- mesh,
- mesh->getIntSection(faultLabel));
- CohesiveTopology::create(*faultMesh, faultBd, mesh,
- mesh->getIntSection(faultLabel),
+ ALE::Obj<ALE::Mesh> faultBoundary = 0;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ CohesiveTopology::createFault(faultMesh, faultBoundary,
+ *mesh, sieveMesh->getIntSection(faultLabel));
+ CohesiveTopology::create(mesh, *faultMesh, faultBoundary,
+ sieveMesh->getIntSection(faultLabel),
faultId,
useLagrangeConstraints);
- CPPUNIT_ASSERT(!faultMesh->isNull());
// Need to copy coordinates from mesh to fault mesh since we are not
// using create() instead of createParallel().
- (*faultMesh)->setRealSection("coordinates",
- mesh->getRealSection("coordinates"));
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
+ CPPUNIT_ASSERT(!faultSieveMesh.isNull());
+ faultSieveMesh->setRealSection("coordinates",
+ sieveMesh->getRealSection("coordinates"));
// Setup databases
spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
@@ -239,21 +264,21 @@
ioSlipTime.filename(slipTimeFilename);
dbSlipTime.ioHandler(&ioSlipTime);
- spatialdata::spatialdb::SimpleDB dbPeakRate("peak rate");
- spatialdata::spatialdb::SimpleIOAscii ioPeakRate;
- ioPeakRate.filename(peakRateFilename);
- dbPeakRate.ioHandler(&ioPeakRate);
+ spatialdata::spatialdb::SimpleDB dbRiseTime("rise time");
+ spatialdata::spatialdb::SimpleIOAscii ioRiseTime;
+ ioRiseTime.filename(peakRateFilename);
+ dbRiseTime.ioHandler(&ioRiseTime);
spatialdata::units::Nondimensional normalizer;
// setup EqKinSrc
slipfn->dbFinalSlip(&dbFinalSlip);
slipfn->dbSlipTime(&dbSlipTime);
- slipfn->dbPeakRate(&dbPeakRate);
+ slipfn->dbRiseTime(&dbRiseTime);
eqsrc->originTime(originTime);
eqsrc->slipfn(slipfn);
- eqsrc->initialize(*faultMesh, &cs, normalizer);
+ eqsrc->initialize(*faultMesh, normalizer);
} // _initialize
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestEqKinSrc.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestEqKinSrc.hh 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestEqKinSrc.hh 2009-04-06 03:55:49 UTC (rev 14601)
@@ -21,7 +21,8 @@
#if !defined(pylith_faults_testeqkinsrc_hh)
#define pylith_faults_testeqkinsrc_hh
-#include "pylith/utils/sievetypes.hh" // USES Mesh
+#include "pylith/faults/faultsfwd.hh" // USES EqKinSrc, BruneSlipFn
+#include "pylith/topology/topologyfwd.hh" // USES Mesh
#include <cppunit/extensions/HelperMacros.h>
@@ -29,12 +30,6 @@
namespace pylith {
namespace faults {
class TestEqKinSrc;
- class EqKinSrc;
- class BruneSlipFn;
-
- namespace _TestEqKinSrc {
- struct DataStruct;
- } // _TestEqKinSrc
} // faults
} // pylith
@@ -79,13 +74,15 @@
/** Initialize EqKinSrc.
*
- * @param faultMesh Fault mesh.
+ * @param mesh Finite-element mesh of domain.
+ * @param faultMesh Finite-element mesh of fault.
* @param eqsrc Earthquake source.
* @param slipfn Slip time function.
* @param originTime Origin time for earthquake rupture.
*/
static
- void _initialize(ALE::Obj<Mesh>* faultMesh,
+ void _initialize(topology::Mesh* mesh,
+ topology::SubMesh* faultMesh,
EqKinSrc* eqsrc,
BruneSlipFn* slipfn,
const double originTime);
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFault.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFault.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFault.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -30,7 +30,7 @@
FaultCohesiveKin fault;
fault.id(id);
- CPPUNIT_ASSERT(id == fault.id());
+ CPPUNIT_ASSERT_EQUAL(id, fault.id());
} // testID
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFault.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFault.hh 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFault.hh 2009-04-06 03:55:49 UTC (rev 14601)
@@ -26,7 +26,6 @@
/// Namespace for pylith package
namespace pylith {
namespace faults {
- class Fault;
class TestFault;
} // faults
} // pylith
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesive.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesive.cc 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesive.cc 2009-04-06 03:55:49 UTC (rev 14601)
@@ -14,13 +14,12 @@
#include "TestFaultCohesive.hh" // Implementation of class methods
-#include "pylith/faults/FaultCohesiveKin.hh" // USES FaultsCohesiveKin
+//#include "pylith/faults/FaultCohesiveKin.hh" // USES FaultsCohesiveKin
#include "pylith/faults/FaultCohesiveDyn.hh" // USES FaultsCohesiveDyn
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/utils/array.hh" // USES int_array, double_array
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
-#include "pylith/feassemble/Quadrature2Din3D.hh" // USES Quadrature2Din3D
#include "data/CohesiveDataLine2.hh" // USES CohesiveDataLine2
#include "data/CohesiveDataTri3.hh" // USES CohesiveDataTri3
@@ -66,6 +65,9 @@
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestFaultCohesive );
// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+
+// ----------------------------------------------------------------------
void
pylith::faults::TestFaultCohesive::testUseFaultMesh(void)
{ // testUseFaultMesh
@@ -420,6 +422,7 @@
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyHex8i
+#if 0
// ----------------------------------------------------------------------
// Test adjustTopology() with 1-D line element for Lagrange
// multipliers.
@@ -474,6 +477,7 @@
FaultCohesiveKin fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyHex8Lagrange
+#endif
// ----------------------------------------------------------------------
// Test adjustTopology().
@@ -482,7 +486,7 @@
const CohesiveData& data,
const bool flipFault)
{ // _testAdjustTopology
- ALE::Obj<Mesh> mesh;
+ topology::Mesh mesh;
meshio::MeshIOAscii iohandler;
iohandler.filename(data.filename);
iohandler.debug(false);
@@ -492,27 +496,31 @@
CPPUNIT_ASSERT(0 != fault);
fault->id(1);
fault->label("fault");
- fault->adjustTopology(mesh, flipFault);
+ fault->adjustTopology(&mesh, flipFault);
//mesh->view(data.filename);
- CPPUNIT_ASSERT_EQUAL(data.cellDim, mesh->getDimension());
+ CPPUNIT_ASSERT_EQUAL(data.cellDim, mesh.dimension());
// Check vertices
- const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
- const ALE::Obj<Mesh::real_section_type>& coordsField =
- mesh->getRealSection("coordinates");
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const ALE::Obj<topology::Mesh::RealSection>& coordsSection =
+ sieveMesh->getRealSection("coordinates");
+ CPPUNIT_ASSERT(!coordsSection.isNull());
const int numVertices = vertices->size();
CPPUNIT_ASSERT_EQUAL(data.numVertices, numVertices);
CPPUNIT_ASSERT_EQUAL(data.spaceDim,
- coordsField->getFiberDimension(*vertices->begin()));
+ coordsSection->getFiberDimension(*vertices->begin()));
int i = 0;
const int spaceDim = data.spaceDim;
- for(Mesh::label_sequence::iterator v_iter =
+ for (SieveMesh::label_sequence::iterator v_iter =
vertices->begin();
v_iter != vertices->end();
++v_iter) {
- const Mesh::real_section_type::value_type *vertexCoords =
- coordsField->restrictPoint(*v_iter);
+ const double* vertexCoords = coordsSection->restrictPoint(*v_iter);
const double tolerance = 1.0e-06;
for (int iDim=0; iDim < spaceDim; ++iDim)
if (data.vertices[i] < 1.0)
@@ -524,22 +532,24 @@
} // for
// check cells
- const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
- const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ CPPUNIT_ASSERT(!sieve.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells = sieveMesh->heightStratum(0);
+ CPPUNIT_ASSERT(!cells.isNull());
const int numCells = cells->size();
CPPUNIT_ASSERT_EQUAL(data.numCells, numCells);
- ALE::ISieveVisitor::PointRetriever<Mesh::sieve_type> pV(sieve->getMaxConeSize());
+ ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> pV(sieve->getMaxConeSize());
int iCell = 0;
i = 0;
- for(Mesh::label_sequence::iterator c_iter = cells->begin();
+ for(SieveMesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
++c_iter) {
- const int numCorners = mesh->getNumCellCorners(*c_iter);
+ const int numCorners = sieveMesh->getNumCellCorners(*c_iter);
CPPUNIT_ASSERT_EQUAL(data.numCorners[iCell++], numCorners);
sieve->cone(*c_iter, pV);
- const Mesh::point_type *cone = pV.getPoints();
+ const SieveMesh::point_type *cone = pV.getPoints();
for(int p = 0; p < pV.getSize(); ++p, ++i) {
CPPUNIT_ASSERT_EQUAL(data.cells[i], cone[p]);
}
@@ -547,41 +557,44 @@
} // for
// check materials
- const ALE::Obj<Mesh::label_type>& labelMaterials =
- mesh->getLabel("material-id");
+ const ALE::Obj<SieveMesh::label_type>& labelMaterials =
+ sieveMesh->getLabel("material-id");
+ CPPUNIT_ASSERT(!labelMaterials.isNull());
const int idDefault = -999;
const int size = numCells;
int_array materialIds(size);
i = 0;
- for(Mesh::label_sequence::iterator c_iter = cells->begin();
+ for (SieveMesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
++c_iter)
- materialIds[i++] = mesh->getValue(labelMaterials, *c_iter, idDefault);
+ materialIds[i++] = sieveMesh->getValue(labelMaterials, *c_iter, idDefault);
for (int iCell=0; iCell < numCells; ++iCell)
CPPUNIT_ASSERT_EQUAL(data.materialIds[iCell], materialIds[iCell]);
// Check groups
const ALE::Obj<std::set<std::string> >& groupNames =
- mesh->getIntSections();
+ sieveMesh->getIntSections();
+ CPPUNIT_ASSERT(!groupNames.isNull());
int iGroup = 0;
int index = 0;
for (std::set<std::string>::const_iterator name=groupNames->begin();
name != groupNames->end();
++name, ++iGroup) {
- const ALE::Obj<int_section_type>& groupField = mesh->getIntSection(*name);
+ const ALE::Obj<topology::Mesh::IntSection>& groupField =
+ sieveMesh->getIntSection(*name);
CPPUNIT_ASSERT(!groupField.isNull());
- const int_section_type::chart_type& chart = groupField->getChart();
- Mesh::point_type firstPoint;
- for(int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
+ const topology::Mesh::IntSection::chart_type& chart = groupField->getChart();
+ SieveMesh::point_type firstPoint;
+ for (topology::Mesh::IntSection::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
if (groupField->getFiberDimension(*c_iter)) {firstPoint = *c_iter; break;}
}
std::string groupType =
- (mesh->height(firstPoint) == 0) ? "cell" : "vertex";
+ (sieveMesh->height(firstPoint) == 0) ? "cell" : "vertex";
const int numPoints = groupField->size();
int_array points(numPoints);
int i = 0;
- for(int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
+ for (topology::Mesh::IntSection::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
if (groupField->getFiberDimension(*c_iter)) points[i++] = *c_iter;
}
@@ -602,7 +615,7 @@
const bool flipFaultA,
const bool flipFaultB)
{ // _testAdjustTopology
- ALE::Obj<Mesh> mesh;
+ topology::Mesh mesh;
meshio::MeshIOAscii iohandler;
iohandler.filename(data.filename);
iohandler.debug(false);
@@ -612,60 +625,67 @@
CPPUNIT_ASSERT(0 != faultA);
faultA->id(1);
faultA->label("faultA");
- faultA->adjustTopology(mesh, flipFaultA);
+ faultA->adjustTopology(&mesh, flipFaultA);
CPPUNIT_ASSERT(0 != faultB);
faultB->id(2);
faultB->label("faultB");
- faultB->adjustTopology(mesh, flipFaultB);
+ faultB->adjustTopology(&mesh, flipFaultB);
- //mesh->view(data.filename);
- CPPUNIT_ASSERT_EQUAL(data.cellDim, mesh->getDimension());
+ //sieveMesh->view(data.filename);
+ CPPUNIT_ASSERT_EQUAL(data.cellDim, mesh.dimension());
// Check vertices
- const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
- const ALE::Obj<Mesh::real_section_type>& coordsField =
- mesh->getRealSection("coordinates");
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const ALE::Obj<topology::Mesh::RealSection>& coordsSection =
+ sieveMesh->getRealSection("coordinates");
+ CPPUNIT_ASSERT(!coordsSection.isNull());
const int numVertices = vertices->size();
CPPUNIT_ASSERT_EQUAL(data.numVertices, numVertices);
CPPUNIT_ASSERT_EQUAL(data.spaceDim,
- coordsField->getFiberDimension(*vertices->begin()));
+ coordsSection->getFiberDimension(*vertices->begin()));
int i = 0;
const int spaceDim = data.spaceDim;
- for(Mesh::label_sequence::iterator v_iter =
+ for(SieveMesh::label_sequence::iterator v_iter =
vertices->begin();
v_iter != vertices->end();
++v_iter) {
- const Mesh::real_section_type::value_type *vertexCoords =
- coordsField->restrictPoint(*v_iter);
+ const double* coordsVertex = coordsSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != coordsVertex);
const double tolerance = 1.0e-06;
for (int iDim=0; iDim < spaceDim; ++iDim)
if (data.vertices[i] < 1.0)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(data.vertices[i++], vertexCoords[iDim],
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.vertices[i++], coordsVertex[iDim],
tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vertexCoords[iDim]/data.vertices[i++],
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, coordsVertex[iDim]/data.vertices[i++],
tolerance);
} // for
// check cells
- const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
- const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ CPPUNIT_ASSERT(!sieve.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cells = sieveMesh->heightStratum(0);
+ CPPUNIT_ASSERT(!cells.isNull());
const int numCells = cells->size();
CPPUNIT_ASSERT_EQUAL(data.numCells, numCells);
- ALE::ISieveVisitor::PointRetriever<Mesh::sieve_type> pV(sieve->getMaxConeSize());
+ ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> pV(sieve->getMaxConeSize());
int iCell = 0;
i = 0;
- //mesh->view(data.filename);
- for(Mesh::label_sequence::iterator c_iter = cells->begin();
+ //sieveMesh->view(data.filename);
+ for(SieveMesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
++c_iter) {
- const int numCorners = mesh->getNumCellCorners(*c_iter);
+ const int numCorners = sieveMesh->getNumCellCorners(*c_iter);
CPPUNIT_ASSERT_EQUAL(data.numCorners[iCell++], numCorners);
sieve->cone(*c_iter, pV);
- const Mesh::point_type *cone = pV.getPoints();
+ const SieveMesh::point_type *cone = pV.getPoints();
for(int p = 0; p < pV.getSize(); ++p, ++i) {
CPPUNIT_ASSERT_EQUAL(data.cells[i], cone[p]);
}
@@ -673,41 +693,44 @@
} // for
// check materials
- const ALE::Obj<Mesh::label_type>& labelMaterials =
- mesh->getLabel("material-id");
+ const ALE::Obj<SieveMesh::label_type>& labelMaterials =
+ sieveMesh->getLabel("material-id");
+ CPPUNIT_ASSERT(!labelMaterials.isNull());
const int idDefault = -999;
const int size = numCells;
int_array materialIds(size);
i = 0;
- for(Mesh::label_sequence::iterator c_iter = cells->begin();
+ for(SieveMesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
++c_iter)
- materialIds[i++] = mesh->getValue(labelMaterials, *c_iter, idDefault);
+ materialIds[i++] = sieveMesh->getValue(labelMaterials, *c_iter, idDefault);
for (int iCell=0; iCell < numCells; ++iCell)
CPPUNIT_ASSERT_EQUAL(data.materialIds[iCell], materialIds[iCell]);
// Check groups
const ALE::Obj<std::set<std::string> >& groupNames =
- mesh->getIntSections();
+ sieveMesh->getIntSections();
+ CPPUNIT_ASSERT(!groupNames.isNull());
int iGroup = 0;
int index = 0;
for (std::set<std::string>::const_iterator name=groupNames->begin();
name != groupNames->end();
++name, ++iGroup) {
- const ALE::Obj<int_section_type>& groupField = mesh->getIntSection(*name);
+ const ALE::Obj<topology::Mesh::IntSection>& groupField =
+ sieveMesh->getIntSection(*name);
CPPUNIT_ASSERT(!groupField.isNull());
- const int_section_type::chart_type& chart = groupField->getChart();
- Mesh::point_type firstPoint;
- for(int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
+ const topology::Mesh::IntSection::chart_type& chart = groupField->getChart();
+ SieveMesh::point_type firstPoint;
+ for(topology::Mesh::IntSection::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
if (groupField->getFiberDimension(*c_iter)) {firstPoint = *c_iter; break;}
}
std::string groupType =
- (mesh->height(firstPoint) == 0) ? "cell" : "vertex";
+ (sieveMesh->height(firstPoint) == 0) ? "cell" : "vertex";
const int numPoints = groupField->size();
int_array points(numPoints);
int i = 0;
- for(int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
+ for(topology::Mesh::IntSection::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
if (groupField->getFiberDimension(*c_iter)) points[i++] = *c_iter;
}
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesive.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesive.hh 2009-04-05 23:12:03 UTC (rev 14600)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/faults/TestFaultCohesive.hh 2009-04-06 03:55:49 UTC (rev 14601)
@@ -23,12 +23,11 @@
#include <cppunit/extensions/HelperMacros.h>
-#include "pylith/utils/sievefwd.hh" // USES PETSc Mesh
+#include "pylith/faults/faultsfwd.hh" // USES PETSc Mesh
/// Namespace for pylith package
namespace pylith {
namespace faults {
- class Fault;
class TestFaultCohesive;
class CohesiveData;
} // faults
@@ -78,11 +77,13 @@
CPPUNIT_TEST( testAdjustTopologyHex8h );
CPPUNIT_TEST( testAdjustTopologyHex8i );
+#if 0
CPPUNIT_TEST( testAdjustTopologyLine2Lagrange );
CPPUNIT_TEST( testAdjustTopologyTri3Lagrange );
CPPUNIT_TEST( testAdjustTopologyQuad4Lagrange );
CPPUNIT_TEST( testAdjustTopologyTet4Lagrange );
CPPUNIT_TEST( testAdjustTopologyHex8Lagrange );
+#endif
CPPUNIT_TEST_SUITE_END();
@@ -194,6 +195,7 @@
/// Test adjustTopology() with 3-D hexahedral element (edge/vertex on fault).
void testAdjustTopologyHex8i(void);
+#if 0
/// Test adjustTopology() with 1-D line element for Lagrange
/// multipliers.
void testAdjustTopologyLine2Lagrange(void);
@@ -213,6 +215,7 @@
/// Test adjustTopology() with 3-D hexahedral element for Lagrange
/// multipliers.
void testAdjustTopologyHex8Lagrange(void);
+#endif
// PROTECTED METHODS //////////////////////////////////////////////////
public :
More information about the CIG-COMMITS
mailing list