[cig-commits] r12776 - in short/3D/PyLith/trunk: libsrc/faults libsrc/meshio pylith/topology unittests/libtests/faults
knepley at geodynamics.org
knepley at geodynamics.org
Tue Sep 2 13:37:54 PDT 2008
Author: knepley
Date: 2008-09-02 13:37:54 -0700 (Tue, 02 Sep 2008)
New Revision: 12776
Modified:
short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh
short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.cc
short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.hh
short/3D/PyLith/trunk/pylith/topology/MeshImporter.py
short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc
Log:
Refactoring to allow faults to be input as entire meshes
Modified: short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.cc 2008-09-02 20:37:54 UTC (rev 12776)
@@ -19,73 +19,45 @@
#include <assert.h> // USES assert()
-// ----------------------------------------------------------------------
void
-pylith::faults::CohesiveTopology::create(ALE::Obj<Mesh>* ifault,
- const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<Mesh::int_section_type>& groupField,
- const int materialId,
- const bool constraintCell)
-{ // create
- assert(0 != ifault);
-
- typedef ALE::SieveAlg<ALE::Mesh> sieveAlg;
+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)
+{
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;
- const int_section_type::chart_type& chart = groupField->getChart();
- PointSet faultVertices; // Vertices on fault
- const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
- *ifault = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
- const ALE::Obj<Mesh::sieve_type> ifaultSieve = new Mesh::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());
- 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 *indices = NULL; // The indices of a face vertex set in a cell
- int oppositeVertex; // For simplices, the vertex opposite a given face
- PointArray origVertices;
- PointArray faceVertices;
- PointArray neighborVertices;
-
- if (!fault->commRank()) {
+ if (!faultSieve->commRank()) {
numCorners = mesh->getNumCellCorners();
faceSize = selection::numFaceVertices(mesh);
indices = new int[faceSize];
}
- // Create set with 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);}
- } // for
-
- const PointSet::const_iterator fvBegin = faultVertices.begin();
- const PointSet::const_iterator fvEnd = faultVertices.end();
-
- int f = sieve->getBaseSize() + sieve->getCapSize();
- int debug = mesh->debug();
- ALE::Obj<PointSet> face = new PointSet();
- PointSet faultCells;
-
- // Create a sieve which captures the fault
- const int fDim = fault->getDimension();
- const Obj<Mesh::arrow_section_type>& orientation = fault->getArrowSection("orientation");
- std::map<int,int*> curElement;
- std::map<int,PointArray> bdVertices;
- std::map<int,PointArray> faultFaces;
- std::map<int,oPointArray> oFaultFaces;
- int curCell = f;
- int curVertex = 0;
- int newElement = curCell + fDim*faultVertices.size();
- int o = 1;
-
curElement[0] = &curVertex;
- curElement[fDim] = &curCell;
- for(int d = 1; d < fDim; d++) {
+ curElement[dim] = &curCell;
+ for(int d = 1; d < dim; d++) {
curElement[d] = &newElement;
}
@@ -120,17 +92,17 @@
if (face->size() == faceSize) {
if (debug) std::cout << " Contains a face on the fault" << std::endl;
ALE::Obj<sieve_type::supportSet> preFace;
- if (fDim < 2) {
+ if (dim < 2) {
preFace = faultSieve->nJoin1(face);
} else {
- preFace = faultSieve->nJoin(face, fDim);
+ 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 (fDim == 0) {
+ if (dim == 0) {
faultSieve->addArrow(*faceVertices.begin(), support[s]);
} else {
faultSieve->addArrow(*preFace->begin(), support[s]);
@@ -138,20 +110,20 @@
} 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[fDim].clear();
+ bdVertices[dim].clear();
for(PointArray::const_iterator v_iter = faceVertices.begin(); v_iter != faceVertices.end(); ++v_iter) {
- bdVertices[fDim].push_back(*v_iter);
+ bdVertices[dim].push_back(*v_iter);
if (debug) std::cout << " Boundary vertex " << *v_iter << std::endl;
}
- if (fDim == 0) {
+ if (dim == 0) {
f = *faceVertices.begin();
}
- if (faceSize != fDim+1) {
+ if (faceSize != dim+1) {
if (debug) std::cout << " Adding hex face " << f << std::endl;
- ALE::SieveBuilder<ALE::Mesh>::buildHexFaces(faultSieve, orientation, fDim, curElement, bdVertices, oFaultFaces, f, o);
+ 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, fDim, curElement, bdVertices, oFaultFaces, f, o);
+ ALE::SieveBuilder<ALE::Mesh>::buildFaces(faultSieve, orientation, dim, curElement, bdVertices, oFaultFaces, f, o);
}
faultSieve->addArrow(f, support[s]);
//faultSieve->view("");
@@ -163,31 +135,92 @@
} // for
sV.clear();
} // for
- fault->setSieve(faultSieve);
- fault->stratify();
- faultCells.clear();
- if (debug) fault->view("Fault mesh");
- const ALE::Obj<ALE::Mesh> faultBd = ALE::Selection<ALE::Mesh>::boundary(fault);
- if (debug) faultBd->view("Fault boundary mesh");
+ if (!faultSieve->commRank()) delete [] indices;
+}
- // Orient the fault sieve
+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;
+
+ if (!faultSieve->commRank()) {
+ faceSize = selection::numFaceVertices(mesh);
+ }
+
+ curElement[0] = &curVertex;
+ curElement[dim] = &curCell;
+ for(int d = 1; d < dim; d++) {
+ curElement[d] = &newElement;
+ }
+
+ // 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]);
+ }
+}
+
+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
- const Mesh::point_type firstFaultCell = *fault->heightStratum(1)->begin();
- const ALE::Obj<Mesh::label_sequence>& fFaces = fault->heightStratum(2);
- const int numFaultFaces = fFaces->size();
- 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
- PointSet flippedCells; // Incorrectly oriented fault cells
- PointSet facesSeen; // Fault faces already considered
- PointSet cellsSeen; // Fault cells already matched
- Obj<PointSet> newCells = new PointSet();
- Obj<PointSet> loopCells = new PointSet();
+ 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<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 (fDim == 0) {
+ if (dim == 0) {
assert(numFaultCorners == faceSize-1);
} else {
assert(numFaultCorners == faceSize);
@@ -233,7 +266,7 @@
if (!seenB) newCells->insert(cellB);
if (debug) std::cout << " neighboring cells " << cellA << " and " << cellB << std::endl;
// In 1D, just check that vertices match
- if (fDim == 1) {
+ 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);
@@ -259,7 +292,7 @@
throw ALE::Exception("Inconsistent mesh orientation: Fault mesh is non-orientable");
}
}
- } else if (fDim == 2) {
+ } 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];
@@ -311,7 +344,7 @@
faultSieve->addArrow(*v_iter, *f_iter, color++);
}
- if (fDim > 1) {
+ 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);
@@ -337,7 +370,7 @@
if (debug) std::cout << " neighboring cells " << cellA << " and " << cellB << std::endl;
// In 1D, just check that vertices match
- if (fDim == 1) {
+ 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);
@@ -369,19 +402,109 @@
}
}
if (debug) fault->view("Oriented Fault mesh");
+}
+// ----------------------------------------------------------------------
+void
+pylith::faults::CohesiveTopology::createFault(Obj<Mesh>& ifault,
+ Obj<ALE::Mesh>& faultBd,
+ const Obj<Mesh>& mesh,
+ const Obj<Mesh::int_section_type>& groupField)
+{ // create
+ typedef ALE::SieveAlg<ALE::Mesh> sieveAlg;
+ typedef ALE::Selection<ALE::Mesh> selection;
+
+ const Obj<sieve_type>& sieve = mesh->getSieve();
+ const Obj<Mesh::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 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 *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;
+
+ if (!fault->commRank()) {
+ numCorners = mesh->getNumCellCorners();
+ faceSize = selection::numFaceVertices(mesh);
+ indices = new int[faceSize];
+ }
+
+ // Create set with vertices on fault
+ const int_section_type::chart_type& chart = groupField->getChart();
+ 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);}
+ } // for
+
+ // Create a sieve which captures the fault
+ const bool vertexFault = true;
+ const int firstFaultCell = sieve->getBaseSize() + sieve->getCapSize();
+
+ createFaultSieveFromVertices(fault->getDimension(), firstFaultCell, faultVertices, mesh, fault->getArrowSection("orientation"), faultSieve);
+ fault->setSieve(faultSieve);
+ fault->stratify();
+ if (debug) fault->view("Fault mesh");
+
+ faultBd = ALE::Selection<ALE::Mesh>::boundary(fault);
+ if (debug) faultBd->view("Fault boundary mesh");
+
+ // Orient the fault sieve
+ orientFaultSieve(fault->getDimension(), mesh, fault->getArrowSection("orientation"), fault);
+
// Convert fault to an IMesh
- Mesh::renumbering_type& renumbering = (*ifault)->getRenumbering();
- (*ifault)->setSieve(ifaultSieve);
- ALE::ISieveConverter::convertMesh(*fault, *(*ifault), renumbering, false);
+ Mesh::renumbering_type& renumbering = ifault->getRenumbering();
+ ifault->setSieve(ifaultSieve);
+ ALE::ISieveConverter::convertMesh(*fault, *ifault, renumbering, false);
renumbering.clear();
- fault = NULL;
- faultSieve = NULL;
+};
+// ----------------------------------------------------------------------
+void
+pylith::faults::CohesiveTopology::create(Obj<Mesh>& ifault,
+ const Obj<ALE::Mesh>& faultBd,
+ const Obj<Mesh>& mesh,
+ const Obj<Mesh::int_section_type>& groupField,
+ const int materialId,
+ const bool constraintCell)
+{ // create
+ typedef ALE::SieveAlg<ALE::Mesh> sieveAlg;
+ typedef ALE::Selection<ALE::Mesh> selection;
+
+ const Obj<sieve_type>& sieve = mesh->getSieve();
+ const Obj<Mesh::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;
+
+ if (!ifault->commRank()) {
+ const Mesh::point_type p = *ifault->heightStratum(1)->begin();
+
+ numCorners = mesh->getNumCellCorners();
+ faceSize = selection::numFaceVertices(mesh);
+ indices = new int[faceSize];
+ numFaultCorners = ifault->getNumCellCorners(p, ifault->depth(p));
+ }
+
// Add new shadow vertices and possibly Lagrange multipler vertices
- const ALE::Obj<Mesh::label_sequence>& fVertices = (*ifault)->depthStratum(0);
- const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
- const ALE::Obj<std::set<std::string> >& groupNames = mesh->getIntSections();
+ const Obj<Mesh::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;
@@ -421,14 +544,14 @@
if (constraintCell) newPoint += numFaultVertices;
// Split the mesh along the fault sieve and create cohesive elements
- const ALE::Obj<Mesh::label_sequence>& faces = (*ifault)->heightStratum(1);
+ const ALE::Obj<Mesh::label_sequence>& faces = ifault->heightStratum(1);
const ALE::Obj<Mesh::label_type>& material = mesh->getLabel("material-id");
const int firstCohesiveCell = newPoint;
PointSet replaceCells;
PointSet noReplaceCells;
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()), ifault->depth()));
for(Mesh::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter, ++newPoint) {
const Mesh::point_type face = *f_iter;
@@ -523,12 +646,14 @@
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, faultBdVertices, replaceCells, noReplaceCells, debug);
+ 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, faultBdVertices, replaceCells, noReplaceCells, debug);
+ classifyCells(sieve, *v_iter, depth, faceSize, firstCohesiveCell, replaceCells, noReplaceCells, debug);
}
// Add new arrows for support of replaced vertices
+ ALE::ISieveVisitor::PointRetriever<sieve_type> sV(std::max(1, sieve->getMaxSupportSize()));
+
for(PointSet::const_iterator v_iter = replaceVertices.begin(); v_iter != replaceVertices.end(); ++v_iter) {
sieve->support(*v_iter, sV);
const Mesh::point_type *support = sV.getPoints();
@@ -645,7 +770,7 @@
}
rVs.clear();
}
- if (!(*ifault)->commRank()) delete [] indices;
+ if (!ifault->commRank()) delete [] indices;
/// THIS IS TOO SLOW mesh->stratify();
const std::string labelName("censored depth");
@@ -666,7 +791,7 @@
// Fix coordinates
const ALE::Obj<real_section_type>& coordinates =
mesh->getRealSection("coordinates");
- const ALE::Obj<Mesh::label_sequence>& fVertices2 = (*ifault)->depthStratum(0);
+ const ALE::Obj<Mesh::label_sequence>& fVertices2 = ifault->depthStratum(0);
if (debug) coordinates->view("Coordinates without shadow vertices");
for(Mesh::label_sequence::iterator v_iter = fVertices2->begin();
@@ -810,7 +935,6 @@
const int depth,
const int faceSize,
const Mesh::point_type& firstCohesiveCell,
- const PointSet& faultBdVertices,
PointSet& replaceCells,
PointSet& noReplaceCells,
const int debug)
Modified: short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.hh 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/libsrc/faults/CohesiveTopology.hh 2008-09-02 20:37:54 UTC (rev 12776)
@@ -160,20 +160,32 @@
// PUBLIC METHODS /////////////////////////////////////////////////////
public :
- /** Create cohesive cells.
+ /** Create the fault mesh.
*
* @param fault Finite-element mesh of fault (output)
* @param mesh Finite-element mesh
* @param faultVertices Vertices assocated with faces of cells defining
* fault surface
+ */
+ static
+ void createFault(Obj<Mesh>& fault,
+ Obj<ALE::Mesh>& faultBd,
+ const Obj<Mesh>& mesh,
+ const Obj<Mesh::int_section_type>& groupField);
+
+ /** Create cohesive cells.
+ *
+ * @param fault Finite-element mesh of fault (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 create(ALE::Obj<Mesh>* fault,
- const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<Mesh::int_section_type>& groupField,
+ void create(Obj<Mesh>& fault,
+ const Obj<ALE::Mesh>& faultBd,
+ const Obj<Mesh>& mesh,
+ const Obj<Mesh::int_section_type>& groupField,
const int materialId,
const bool constraintCell =false);
@@ -259,11 +271,30 @@
const int depth,
const int faceSize,
const Mesh::point_type& firstCohesiveCell,
- const PointSet& faultBdVertices,
PointSet& replaceCells,
PointSet& noReplaceCells,
const int debug);
+ static void 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);
+
+ static void 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);
+
+ static void orientFaultSieve(const int dim,
+ const Obj<Mesh>& mesh,
+ const Obj<ALE::Mesh::arrow_section_type>& orientation,
+ const Obj<ALE::Mesh>& fault);
}; // class CohesiveTopology
#endif // pylith_faults_cohesivetopology_hh
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2008-09-02 20:37:54 UTC (rev 12776)
@@ -16,6 +16,8 @@
#include "CohesiveTopology.hh" // USES CohesiveTopology::create()
+#include "pylith/meshio/MeshIOLagrit.hh" // USES MeshIOLagrit::readFault()
+
#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
#include "pylith/utils/array.hh" // USES double_array
@@ -60,28 +62,41 @@
pylith::faults::FaultCohesive::adjustTopology(const ALE::Obj<Mesh>& mesh)
{ // adjustTopology
assert(std::string("") != label());
+ Obj<Mesh> faultMesh = NULL;
+ Obj<ALE::Mesh> faultBd = NULL;
- if (!_useFaultMesh) {
- // Use group of vertices to define fault.
+ if (_useFaultMesh) {
+ const int faultDim = 2;
+
+ //MPI_Bcast(&faultDim, 1, MPI_INT, 0, comm);
+ faultMesh = new Mesh(mesh->comm(), faultDim, mesh->debug());
+ pylith::meshio::MeshIOLagrit::readFault(_faultMeshFilename, faultMesh, faultBd);
+
+ // Get group of vertices associated with fault
+ const ALE::Obj<int_section_type>& groupField =
+ mesh->getIntSection(label());
+ Obj<ALE::Mesh> faultBd = NULL;
+
+ CohesiveTopology::create(faultMesh, faultBd, mesh, groupField, id(), _useLagrangeConstraints());
+ } else {
if (!mesh->hasIntSection(label())) {
std::ostringstream msg;
msg << "Mesh missing group of vertices '" << label()
- << " for fault interface condition.";
+ << " for fault interface condition.";
throw std::runtime_error(msg.str());
} // if
-
+
// Get group of vertices associated with fault
const ALE::Obj<int_section_type>& groupField =
mesh->getIntSection(label());
assert(!groupField.isNull());
-
- ALE::Obj<Mesh> faultMesh;
- CohesiveTopology::create(&faultMesh, mesh, groupField, id(),
- _useLagrangeConstraints());
- } else {
- // Use fault mesh to define fault.
- std::cout << "ADD FAULT MESH ADJUSTING TOPOLOGY STUFF HERE." << std::endl;
- } // else
+
+ faultMesh = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
+
+ CohesiveTopology::createFault(faultMesh, faultBd, mesh, groupField);
+
+ CohesiveTopology::create(faultMesh, faultBd, mesh, groupField, id(), _useLagrangeConstraints());
+ }
} // adjustTopology
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc 2008-09-02 20:37:54 UTC (rev 12776)
@@ -14,6 +14,8 @@
#include "MeshIO.hh" // implementation of class methods
+#include "Selection.hh" // USES boundary()
+
#include "pylith/utils/array.hh" // USES double_array, int_array
#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
@@ -101,26 +103,29 @@
MPI_Comm_rank(comm, &rank);
// Memory debugging
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.setDebug(_debug);
+ logger.setDebug(_debug/2);
logger.stagePush("MeshCreation");
if (!rank) {
assert(coordinates.size() == numVertices*spaceDim);
assert(cells.size() == numCells*numCorners);
- //if (!_interpolate) {
if (0) {
+ ///if (!_interpolate) {
// Create the ISieve
sieve->setChart(Mesh::sieve_type::chart_type(0, numCells+numVertices));
// Set cone and support sizes
for(int c = 0; c < numCells; ++c) {sieve->setConeSize(c, numCorners);}
- sieve->symmetrizeSizes(numCells, numCorners, const_cast<int*>(&cells[0]));
+ sieve->symmetrizeSizes(numCells, numCorners, const_cast<int*>(&cells[0]), numCells);
// Allocate point storage
sieve->allocate();
// Fill up cones
- int *cone = new int[numCorners];
+ int *cone = new int[numCorners];
+ int *coneO = new int[numCorners];
+ for(int v = 0; v < numCorners; ++v) {coneO[v] = 1;}
for(int c = 0; c < numCells; ++c) {
for(int v = 0; v < numCorners; ++v) cone[v] = cells[c*numCorners+v]+numCells;
sieve->setCone(cone, c);
+ sieve->setConeOrientation(coneO, c);
}
delete [] cone;
// Symmetrize to fill up supports
@@ -190,6 +195,98 @@
} // _buildMesh
// ----------------------------------------------------------------------
+// Set vertices in fault mesh.
+void
+pylith::meshio::MeshIO::_buildFaultMesh(const double_array& coordinates,
+ const int numVertices,
+ const int spaceDim,
+ const int_array& cells,
+ const int numCells,
+ const int numCorners,
+ const int_array& faceCells,
+ const int meshDim,
+ const Obj<Mesh>& fault,
+ Obj<ALE::Mesh>& faultBd)
+{ // _buildFaultMesh
+ MPI_Comm comm = PETSC_COMM_WORLD;
+ int dim = meshDim;
+ int rank;
+
+ assert(!fault.isNull());
+
+ ALE::Obj<sieve_type> sieve = new sieve_type(fault->comm());
+ fault->setDebug(fault->debug());
+ fault->setSieve(sieve);
+
+ MPI_Comm_rank(comm, &rank);
+ // Memory debugging
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.setDebug(fault->debug()/2);
+
+ logger.stagePush("FaultCreation");
+ if (!rank) {
+ assert(coordinates.size() == numVertices*spaceDim);
+ assert(cells.size() == numCells*numCorners);
+ {
+ ALE::Obj<ALE::Mesh::sieve_type> s = new ALE::Mesh::sieve_type(sieve->comm(), sieve->debug());
+
+ ALE::SieveBuilder<ALE::Mesh>::buildTopology(s, meshDim,
+ numCells,
+ const_cast<int*>(&cells[0]),
+ numVertices,
+ true,
+ numCorners,
+ 0,
+ fault->getArrowSection("orientation"));
+
+ // Add in cells
+ for(int c = 0; c < numCells; ++c) {
+ s->addArrow(c, faceCells[c*2+0]);
+ s->addArrow(c, faceCells[c*2+1]);
+ }
+
+ Mesh::renumbering_type& renumbering = fault->getRenumbering();
+ ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering, false);
+
+ Obj<ALE::Mesh> tmpMesh = new ALE::Mesh(fault->comm(), dim, fault->debug());
+ faultBd = ALE::Selection<ALE::Mesh>::boundary(tmpMesh);
+ }
+ logger.stagePop();
+ logger.stagePush("FaultStratification");
+ fault->stratify();
+ logger.stagePop();
+ } else {
+ logger.stagePush("FaultStratification");
+ fault->getSieve()->setChart(sieve_type::chart_type());
+ fault->getSieve()->allocate();
+ fault->stratify();
+ logger.stagePop();
+ }
+
+#if defined(ALE_MEM_LOGGING)
+ std::cout
+ << std::endl
+ << "FaultCreation " << logger.getNumAllocations("FaultCreation")
+ << " allocations " << logger.getAllocationTotal("FaultCreation")
+ << " bytes" << std::endl
+
+ << "FaultCreation " << logger.getNumDeallocations("FaultCreation")
+ << " deallocations " << logger.getDeallocationTotal("FaultCreation")
+ << " bytes" << std::endl
+
+ << "FaultStratification " << logger.getNumAllocations("FaultStratification")
+ << " allocations " << logger.getAllocationTotal("FaultStratification")
+ << " bytes" << std::endl
+
+ << "FaultStratification " << logger.getNumDeallocations("FaultStratification")
+ << " deallocations " << logger.getDeallocationTotal("FaultStratification")
+ << " bytes" << std::endl << std::endl;
+#endif
+
+ ALE::SieveBuilder<Mesh>::buildCoordinates(fault, spaceDim, &coordinates[0]);
+} // _buildFaultMesh
+
+// ----------------------------------------------------------------------
// Get coordinates of vertices in mesh.
void
pylith::meshio::MeshIO::_getVertices(double_array* coordinates,
@@ -278,6 +375,8 @@
assert(0 != _mesh);
assert(!_mesh->isNull());
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("MaterialsCreation");
const ALE::Obj<Mesh::label_type>& labelMaterials =
(*_mesh)->createLabel("material-id");
if (!(*_mesh)->commRank()) {
@@ -299,6 +398,12 @@
++e_iter)
(*_mesh)->setValue(labelMaterials, *e_iter, materialIds[i++]);
} // if
+ logger.stagePop();
+
+#if defined(ALE_MEM_LOGGING)
+ std::cout << "MaterialsCreation " << logger.getNumAllocations("MaterialsCreation") << " allocations " << logger.getAllocationTotal("MaterialsCreation") << " bytes" << std::endl;
+ std::cout << "MaterialsCreation " << logger.getNumDeallocations("MaterialsCreation") << " deallocations " << logger.getDeallocationTotal("MaterialsCreation") << " bytes" << std::endl;
+#endif
} // _setMaterials
// ----------------------------------------------------------------------
@@ -339,6 +444,8 @@
assert(0 != _mesh);
assert(!_mesh->isNull());
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("GroupCreation");
const ALE::Obj<int_section_type>& groupField = (*_mesh)->getIntSection(name);
assert(!groupField.isNull());
@@ -355,6 +462,12 @@
groupField->setFiberDimension(numCells+points[i], 1);
} // if/else
(*_mesh)->allocate(groupField);
+ logger.stagePop();
+
+#if defined(ALE_MEM_LOGGING)
+ std::cout << "GroupCreation " << logger.getNumAllocations("GroupCreation") << " allocations " << logger.getAllocationTotal("GroupCreation") << " bytes" << std::endl;
+ std::cout << "GroupCreation " << logger.getNumDeallocations("GroupCreation") << " deallocations " << logger.getDeallocationTotal("GroupCreation") << " bytes" << std::endl;
+#endif
} // _setGroup
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIO.hh 2008-09-02 20:37:54 UTC (rev 12776)
@@ -18,6 +18,8 @@
#include "pylith/utils/sievetypes.hh" // USES Obj, PETSc Mesh
+#include "Mesh.hh" // Needs ALE::Mesh
+
namespace pylith {
namespace meshio {
class MeshIO;
@@ -117,6 +119,18 @@
const int numCorners,
const int meshDim);
+ static void
+ _buildFaultMesh(const double_array& coordinates,
+ const int numVertices,
+ const int spaceDim,
+ const int_array& cells,
+ const int numCells,
+ const int numCorners,
+ const int_array& faceCells,
+ const int meshDim,
+ const Obj<Mesh>& fault,
+ Obj<ALE::Mesh>& faultBd);
+
/** Get information about vertices in mesh.
*
* Method caller is responsible for memory management.
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.cc 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.cc 2008-09-02 20:37:54 UTC (rev 12776)
@@ -43,6 +43,132 @@
{ // destructor
} // destructor
+void
+pylith::meshio::MeshIOLagrit::readFault(const std::string filename, const Obj<Mesh>& fault, Obj<ALE::Mesh>& faultBd) {
+ int faultDim = 2;
+ int fSpaceDim = 0;
+ int numFVertices = 0;
+ int numFCells = 0;
+ int numFCorners = 0;
+ double_array fCoordinates;
+ int_array fCells;
+ int_array fMaterialIds;
+ int_array faceCells;
+ int_array vertexIDs;
+
+ if (!fault->commRank()) {
+ std::ifstream fin(filename.c_str(), 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;
+
+ fSpaceDim = 3;
+ // Section 1: <num_nodes> <num_cells> <num_ndata> <num_cdata> <num_mdata>
+ fin >> numFVertices;
+ fin >> numFCells;
+ fin >> numNodeData;
+ fin >> numCellData;
+ fin >> numModelData;
+ // Section 2: <node_id 1> <x> <y> <z>
+ fCoordinates.resize(numFVertices * fSpaceDim);
+ for(int v = 0; v < numFVertices; ++v) {
+ int id;
+
+ fin >> id;
+ for(int d = 0; d < fSpaceDim; ++d) {
+ fin >> fCoordinates[v*fSpaceDim + d];
+ }
+ }
+ // Section 3: <cell_id 1> <mat_id> <cell_type> <cell_vert 1> ... <cell_vert n>
+ fMaterialIds.resize(numFCells);
+ for(int c = 0; c < numFCells; ++c) {
+ std::string cellType;
+ int cellID;
+
+ fin >> cellID;
+ fin >> fMaterialIds[c];
+ fin >> cellType;
+ if (cellType == "tri") {
+ numFCorners = 3;
+ } else if (cellType == "quad") {
+ numFCorners = 4;
+ } else {
+ std::ostringstream msg;
+ msg << "Unknown cell type " << cellType << "while reading INP file.";
+ throw std::runtime_error(msg.str());
+ }
+ if (c == 0) fCells.resize(numFCells*numFCorners);
+ for(int i = 0; i < numFCorners; ++i) fin >> fCells[c*numFCorners+i];
+ }
+ // Section 4: <num_comp for node data> <size comp 1> <size comp 2>...<size comp n>
+ int numComponents, totalSize;
+
+ fin >> numComponents;
+ totalSize = 0;
+ for(int i = 0; i < numComponents; ++i) {
+ int compSize;
+ fin >> compSize;
+ totalSize += compSize;
+ }
+ // Section 5: <node_comp_label 1> , <units_label 1>
+ for(int c = 0; c < numComponents; ++c) {
+ std::string label, typeName;
+
+ fin >> label;
+ fin >> typeName;
+ }
+ // Section 6: <node_id 1> <node_data 1> ... <node_data num_ndata>
+ vertexIDs.resize(numFVertices);
+ for(int v = 0; v < numFVertices; ++v) {
+ int id, dummy;
+
+ fin >> id;
+ fin >> vertexIDs[v];
+ fin >> dummy;
+ fin >> dummy;
+ }
+ // Section 7: <num_comp for cell's data> <size comp 1> <size comp 2>...<size comp n>
+ fin >> numComponents;
+ totalSize = 0;
+ for(int i = 0; i < numComponents; ++i) {
+ int compSize;
+ fin >> compSize;
+ totalSize += compSize;
+ }
+ // Section 8: <cell_component_label 1> , <units_label 1>
+ for(int c = 0; c < numComponents; ++c) {
+ std::string label, typeName;
+
+ fin >> label;
+ fin >> typeName;
+ }
+ // Section 9: <cell_id 1> <cell_data 1> ... <cell_data num_cdata>
+ faceCells.resize(numFCells*2);
+ for(int c = 0; c < numFCells; ++c) {
+ int id, dummy;
+
+ fin >> id;
+ fin >> dummy;
+ fin >> dummy;
+ fin >> faceCells[c*2+0];
+ fin >> faceCells[c*2+1];
+ fin >> dummy;
+ }
+
+ // Renumber vertices
+ for(int c = 0; c < numFCells; ++c) {
+ for(int corner = 0; corner < numFCorners; ++corner) {
+ fCells[c*numFCorners+corner] = vertexIDs[fCells[c*numFCorners+corner]];
+ }
+ }
+ }
+
+ _buildFaultMesh(fCoordinates, numFVertices, fSpaceDim, fCells, numFCells, numFCorners, faceCells, faultDim, fault, faultBd);
+}
+
// ----------------------------------------------------------------------
// Unpickle mesh
void
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.hh 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIOLagrit.hh 2008-09-02 20:37:54 UTC (rev 12776)
@@ -111,6 +111,8 @@
*/
bool isRecordHeader32Bit(void) const;
+ static void readFault(const std::string filename, const Obj<Mesh>& fault, Obj<ALE::Mesh>& faultBd);
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/pylith/topology/MeshImporter.py
===================================================================
--- short/3D/PyLith/trunk/pylith/topology/MeshImporter.py 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/pylith/topology/MeshImporter.py 2008-09-02 20:37:54 UTC (rev 12776)
@@ -72,11 +72,14 @@
"""
Hook for creating mesh.
"""
+ from pylith.utils.profiling import resourceUsageString
+
self._setupLogging()
logEvent = "%screate" % self._loggingPrefix
self._logger.eventBegin(logEvent)
mesh = self.importer.read(self.debug, self.interpolate)
+ self._debug.log(resourceUsageString())
self._info.log("Adjusting topology.")
self._adjustTopology(mesh, faults)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestBruneSlipFn.cc 2008-09-02 20:37:54 UTC (rev 12776)
@@ -379,9 +379,15 @@
// Create fault mesh
const bool useLagrangeConstraints = true;
- CohesiveTopology::create(faultMesh, mesh,
- mesh->getIntSection(faultLabel),
- faultId);
+ (*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),
+ faultId,
+ useLagrangeConstraints);
CPPUNIT_ASSERT(!faultMesh->isNull());
// Need to copy coordinates from mesh to fault mesh since we are not
// using create() instead of createParallel().
@@ -432,11 +438,16 @@
cs.setSpaceDim(spaceDim);
// Create fault mesh
- ALE::Obj<Mesh> faultMesh;
+ ALE::Obj<Mesh> faultMesh = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
+ ALE::Obj<ALE::Mesh> faultBd = NULL;
const bool useLagrangeConstraints = true;
- CohesiveTopology::create(&faultMesh, mesh,
- mesh->getIntSection(data.faultLabel),
- data.faultId);
+ CohesiveTopology::createFault(faultMesh, faultBd,
+ mesh,
+ mesh->getIntSection(data.faultLabel));
+ CohesiveTopology::create(faultMesh, faultBd, mesh,
+ mesh->getIntSection(data.faultLabel),
+ data.faultId,
+ useLagrangeConstraints);
CPPUNIT_ASSERT(!faultMesh.isNull());
// Need to copy coordinates from mesh to fault mesh since we are not
// using create() instead of createParallel().
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestConstRateSlipFn.cc 2008-09-02 20:37:54 UTC (rev 12776)
@@ -300,9 +300,15 @@
// Create fault mesh
const bool useLagrangeConstraints = true;
- CohesiveTopology::create(faultMesh, mesh,
- mesh->getIntSection(faultLabel),
- faultId);
+ (*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),
+ faultId,
+ useLagrangeConstraints);
CPPUNIT_ASSERT(!faultMesh->isNull());
// Need to copy coordinates from mesh to fault mesh since we are not
// using create() instead of createParallel().
@@ -347,11 +353,16 @@
cs.setSpaceDim(spaceDim);
// Create fault mesh
- ALE::Obj<Mesh> faultMesh;
+ ALE::Obj<Mesh> faultMesh = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
+ ALE::Obj<ALE::Mesh> faultBd = NULL;
const bool useLagrangeConstraints = true;
- CohesiveTopology::create(&faultMesh, mesh,
- mesh->getIntSection(data.faultLabel),
- data.faultId);
+ CohesiveTopology::createFault(faultMesh, faultBd,
+ mesh,
+ mesh->getIntSection(data.faultLabel));
+ CohesiveTopology::create(faultMesh, faultBd, mesh,
+ mesh->getIntSection(data.faultLabel),
+ data.faultId,
+ useLagrangeConstraints);
CPPUNIT_ASSERT(!faultMesh.isNull());
// Need to copy coordinates from mesh to fault mesh since we are not
// using create() instead of createParallel().
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestEqKinSrc.cc 2008-09-02 20:37:54 UTC (rev 12776)
@@ -212,9 +212,15 @@
// Create fault mesh
const bool useLagrangeConstraints = true;
- CohesiveTopology::create(faultMesh, mesh,
- mesh->getIntSection(faultLabel),
- faultId);
+ (*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),
+ faultId,
+ useLagrangeConstraints);
CPPUNIT_ASSERT(!faultMesh->isNull());
// Need to copy coordinates from mesh to fault mesh since we are not
// using create() instead of createParallel().
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc 2008-09-02 19:31:22 UTC (rev 12775)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestStepSlipFn.cc 2008-09-02 20:37:54 UTC (rev 12776)
@@ -291,9 +291,15 @@
// Create fault mesh
const bool useLagrangeConstraints = true;
- CohesiveTopology::create(faultMesh, mesh,
- mesh->getIntSection(faultLabel),
- faultId);
+ (*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),
+ faultId,
+ useLagrangeConstraints);
CPPUNIT_ASSERT(!faultMesh->isNull());
// Need to copy coordinates from mesh to fault mesh since we are not
// using create() instead of createParallel().
@@ -338,11 +344,16 @@
cs.setSpaceDim(spaceDim);
// Create fault mesh
- ALE::Obj<Mesh> faultMesh;
+ ALE::Obj<Mesh> faultMesh = new Mesh(mesh->comm(), mesh->getDimension()-1, mesh->debug());
+ ALE::Obj<ALE::Mesh> faultBd = NULL;
const bool useLagrangeConstraints = true;
- CohesiveTopology::create(&faultMesh, mesh,
- mesh->getIntSection(data.faultLabel),
- data.faultId);
+ CohesiveTopology::createFault(faultMesh, faultBd,
+ mesh,
+ mesh->getIntSection(data.faultLabel));
+ CohesiveTopology::create(faultMesh, faultBd, mesh,
+ mesh->getIntSection(data.faultLabel),
+ data.faultId,
+ useLagrangeConstraints);
CPPUNIT_ASSERT(!faultMesh.isNull());
// Need to copy coordinates from mesh to fault mesh since we are not
// using create() instead of createParallel().
More information about the cig-commits
mailing list