[cig-commits] r15544 - in short/3D/PyLith/trunk: libsrc/faults libsrc/meshio modulesrc/faults pylith/topology
brad at geodynamics.org
brad at geodynamics.org
Fri Aug 14 15:13:56 PDT 2009
Author: brad
Date: 2009-08-14 15:13:56 -0700 (Fri, 14 Aug 2009)
New Revision: 15544
Modified:
short/3D/PyLith/trunk/libsrc/faults/Fault.cc
short/3D/PyLith/trunk/libsrc/faults/Fault.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
short/3D/PyLith/trunk/libsrc/meshio/UCDFaultFile.cc
short/3D/PyLith/trunk/libsrc/meshio/UCDFaultFile.hh
short/3D/PyLith/trunk/modulesrc/faults/Fault.i
short/3D/PyLith/trunk/modulesrc/faults/FaultCohesive.i
short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py
Log:
Attempted to fix getting number of vertices for faults when using UCD files with fault meshes.
Modified: short/3D/PyLith/trunk/libsrc/faults/Fault.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/Fault.cc 2009-08-14 18:54:09 UTC (rev 15543)
+++ short/3D/PyLith/trunk/libsrc/faults/Fault.cc 2009-08-14 22:13:56 UTC (rev 15544)
@@ -48,16 +48,4 @@
return *_faultMesh;
} // faultMesh
-// ----------------------------------------------------------------------
-// Get mesh associated with fault fields.
-int
-pylith::faults::Fault::faultSize(topology::Mesh* const mesh) const
-{ // faultSize
- assert(0 != mesh);
- assert(std::string("") != label());
- const ALE::Obj<topology::Mesh::IntSection>& groupField =
- mesh->sieveMesh()->getIntSection(label());
- return groupField->size();
-} // faultSize
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/Fault.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/Fault.hh 2009-08-14 18:54:09 UTC (rev 15543)
+++ short/3D/PyLith/trunk/libsrc/faults/Fault.hh 2009-08-14 22:13:56 UTC (rev 15544)
@@ -80,9 +80,10 @@
/** Get the number of vertices on the fault.
*
* @param mesh PETSc mesh
- * @return faults size
+ * @return Number of vertices on the fault.
*/
- int faultSize(topology::Mesh* const mesh) const;
+ virtual
+ int numVertices(const topology::Mesh& mesh) const = 0;
/** Adjust mesh topology for fault implementation.
*
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2009-08-14 18:54:09 UTC (rev 15543)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.cc 2009-08-14 22:13:56 UTC (rev 15544)
@@ -63,6 +63,36 @@
} // faultMeshFilename
// ----------------------------------------------------------------------
+// Get number of vertices in fault.
+int
+pylith::faults::FaultCohesive::numVertices(const topology::Mesh& mesh) const
+{ // numVertices
+ int nvertices = 0;
+
+ if (!_useFaultMesh) {
+ // Get group of vertices associated with fault
+ const ALE::Obj<topology::Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ if (!sieveMesh->hasIntSection(label())) {
+ std::ostringstream msg;
+ msg << "Mesh missing group of vertices '" << label()
+ << "' for fault interface condition.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ assert(std::string("") != label());
+ const ALE::Obj<topology::Mesh::IntSection>& groupField =
+ mesh.sieveMesh()->getIntSection(label());
+ nvertices = groupField->size();
+ } else {
+ assert(3 == mesh.dimension());
+ nvertices = meshio::UCDFaultFile::numVertices(_faultMeshFilename.c_str());
+ } // else
+
+ return nvertices;
+} // numVertices
+
+// ----------------------------------------------------------------------
// Adjust mesh topology for fault implementation.
void
pylith::faults::FaultCohesive::adjustTopology(topology::Mesh* const mesh,
@@ -72,27 +102,37 @@
{ // adjustTopology
assert(0 != mesh);
assert(std::string("") != label());
-
+
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());
- if (_useFaultMesh) {
+ if (!_useFaultMesh) {
+ if (!sieveMesh->hasIntSection(label())) {
+ std::ostringstream msg;
+ msg << "Mesh missing group of vertices '" << label()
+ << "' for fault interface condition.";
+ throw std::runtime_error(msg.str());
+ } // if
+ const ALE::Obj<topology::Mesh::IntSection>& groupField =
+ sieveMesh->getIntSection(label());
+ assert(!groupField.isNull());
+ CohesiveTopology::createFault(&faultMesh, faultBoundary, *mesh, groupField,
+ flipFault);
+
+ CohesiveTopology::create(mesh, faultMesh, faultBoundary, groupField, id(),
+ *firstFaultVertex, *firstFaultCell, useLagrangeConstraints());
+
+ } else {
const int faultDim = 2;
assert(3 == mesh->dimension());
meshio::UCDFaultFile::read(_faultMeshFilename.c_str(),
&faultMesh, faultBoundary, *mesh);
- // TODO: This is a kludge. We need to make UCD faults report their size
- const int numFaultVertices = faultMesh.sieveMesh()->depthStratum(0)->size();
- *firstFaultCell += numFaultVertices;
- if (useLagrangeConstraints()) {
- *firstFaultCell += numFaultVertices;
- }
-
// Set coordinates in fault mesh
const ALE::Obj<topology::SubMesh::SieveMesh>& faultSieveMesh =
faultMesh.sieveMesh();
@@ -103,7 +143,7 @@
if (!sieveMesh->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
const ALE::Obj<topology::Mesh::IntSection>& groupField =
@@ -111,21 +151,6 @@
assert(!groupField.isNull());
CohesiveTopology::create(mesh, faultMesh, faultBoundary, groupField, id(),
*firstFaultVertex, *firstFaultCell, useLagrangeConstraints());
- } else {
- if (!sieveMesh->hasIntSection(label())) {
- std::ostringstream msg;
- msg << "Mesh missing group of vertices '" << label()
- << " for fault interface condition.";
- throw std::runtime_error(msg.str());
- } // if
- const ALE::Obj<topology::Mesh::IntSection>& groupField =
- sieveMesh->getIntSection(label());
- assert(!groupField.isNull());
- CohesiveTopology::createFault(&faultMesh, faultBoundary, *mesh, groupField,
- flipFault);
-
- CohesiveTopology::create(mesh, faultMesh, faultBoundary, groupField, id(),
- *firstFaultVertex, *firstFaultCell, useLagrangeConstraints());
} // if/else
} // adjustTopology
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2009-08-14 18:54:09 UTC (rev 15543)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2009-08-14 22:13:56 UTC (rev 15544)
@@ -55,6 +55,13 @@
*/
void faultMeshFilename(const char* filename);
+ /** Get the number of vertices on the fault.
+ *
+ * @param mesh PETSc mesh
+ * @return Number of vertices on the fault.
+ */
+ int numVertices(const topology::Mesh& mesh) const;
+
/** Adjust mesh topology for fault implementation.
*
* If firstFaultVertex == 0, then firstFaultVertex is set to the first point
Modified: short/3D/PyLith/trunk/libsrc/meshio/UCDFaultFile.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/UCDFaultFile.cc 2009-08-14 18:54:09 UTC (rev 15543)
+++ short/3D/PyLith/trunk/libsrc/meshio/UCDFaultFile.cc 2009-08-14 22:13:56 UTC (rev 15544)
@@ -22,9 +22,38 @@
#include <petsc.h> // USES MPI_Comm
#include <cassert> // USES assert()
+#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
// ----------------------------------------------------------------------
+int
+pylith::meshio::UCDFaultFile::numVertices(const char* filename)
+{ // read
+ int nvertices = 0;
+
+ 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
+
+ // Section 1: <num_nodes> <num_cells> <num_ndata> <num_cdata> <num_mdata>
+ fin >> nvertices;
+ fin.close();
+
+ if (nvertices <= 0) {
+ std::ostringstream msg;
+ msg << "Number of vertices (" << nvertices << ") in ASCII INP file '"
+ << filename << "' must be positive.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ return nvertices;
+} // numVertices
+
+
+// ----------------------------------------------------------------------
void
pylith::meshio::UCDFaultFile::read(const char* filename,
topology::SubMesh* faultMesh,
Modified: short/3D/PyLith/trunk/libsrc/meshio/UCDFaultFile.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/UCDFaultFile.hh 2009-08-14 18:54:09 UTC (rev 15543)
+++ short/3D/PyLith/trunk/libsrc/meshio/UCDFaultFile.hh 2009-08-14 22:13:56 UTC (rev 15544)
@@ -40,7 +40,26 @@
// PUBLIC METHODS ///////////////////////////////////////////////////////
public :
+ /** Get number of vertices in fault mesh.
+ *
+ * @param filename Name of UCD file.
+ * @returns Number of vertices.
+ */
static
+ int numVertices(const char* filename);
+
+ /** Read fault mesh from UCD file. The fault mesh must be associated
+ * with the mesh of the domain.
+ *
+ * :WARNING: Only works for 3-D domains with tetrahedral cells and
+ * UCD files with specific information.
+ *
+ * @param filename Name of UCD file.
+ * @param faultMesh Fault mesh.
+ * @param faultBoundary Boundary of fault mesh.
+ * @param mesh Domain mesh.
+ */
+ static
void read(const char* filename,
topology::SubMesh* faultMesh,
ALE::Obj<ALE::Mesh>& faultBoundary,
Modified: short/3D/PyLith/trunk/modulesrc/faults/Fault.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/Fault.i 2009-08-14 18:54:09 UTC (rev 15543)
+++ short/3D/PyLith/trunk/modulesrc/faults/Fault.i 2009-08-14 22:13:56 UTC (rev 15544)
@@ -62,10 +62,11 @@
/** Get the number of vertices on the fault.
*
* @param mesh PETSc mesh
- * @return faults size
+ * @return Number of vertices on the fault.
*/
- int faultSize(topology::Mesh* const mesh) const;
-
+ virtual
+ int numVertices(const topology::Mesh& mesh) const = 0;
+
/** Adjust mesh topology for fault implementation.
*
* @param mesh PETSc mesh
Modified: short/3D/PyLith/trunk/modulesrc/faults/FaultCohesive.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/FaultCohesive.i 2009-08-14 18:54:09 UTC (rev 15543)
+++ short/3D/PyLith/trunk/modulesrc/faults/FaultCohesive.i 2009-08-14 22:13:56 UTC (rev 15544)
@@ -49,6 +49,13 @@
*/
void faultMeshFilename(const char* filename);
+ /** Get the number of vertices on the fault.
+ *
+ * @param mesh PETSc mesh
+ * @return Number of vertices on the fault.
+ */
+ int numVertices(const topology::Mesh& mesh) const;
+
/** Adjust mesh topology for fault implementation.
*
* @param mesh PETSc mesh.
Modified: short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py
===================================================================
--- short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py 2009-08-14 18:54:09 UTC (rev 15543)
+++ short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py 2009-08-14 22:13:56 UTC (rev 15544)
@@ -95,12 +95,15 @@
firstFaultVertex = 0
firstFaultCell = 0
for interface in interfaces:
- firstFaultCell += interface.faultSize(mesh)
+ nvertices = interface.numVertices(mesh)
+ firstFaultCell += nvertices
if interface.useLagrangeConstraints():
- firstFaultCell += interface.faultSize(mesh)
+ firstFaultCell += nvertices
for interface in interfaces:
- self._info.log("Adjusting topology for fault '%s' of size %d." % (interface.label, interface.faultSize(mesh)))
- firstFaultVertex, firstFaultCell = interface.adjustTopology(mesh, firstFaultVertex, firstFaultCell)
+ self._info.log("Adjusting topology for fault '%s' with %d vertices." % \
+ (interface.label, nvertices))
+ firstFaultVertex, firstFaultCell = \
+ interface.adjustTopology(mesh, firstFaultVertex, firstFaultCell)
self._eventLogger.eventEnd(logEvent)
return
More information about the CIG-COMMITS
mailing list