[cig-commits] r11636 - in short/3D/PyLith/trunk: . libsrc libsrc/topology modulesrc/topology pylith/problems pylith/topology tests/2d/bar_quad4 unittests/libtests/topology unittests/pytests/topology
brad at geodynamics.org
brad at geodynamics.org
Sat Mar 29 20:25:44 PDT 2008
Author: brad
Date: 2008-03-29 20:25:44 -0700 (Sat, 29 Mar 2008)
New Revision: 11636
Added:
short/3D/PyLith/trunk/libsrc/topology/MeshOps.cc
short/3D/PyLith/trunk/libsrc/topology/MeshOps.hh
short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.hh
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/topology/Makefile.am
short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src
short/3D/PyLith/trunk/pylith/problems/Problem.py
short/3D/PyLith/trunk/pylith/topology/Mesh.py
short/3D/PyLith/trunk/tests/2d/bar_quad4/dislocation_dyn.cfg
short/3D/PyLith/trunk/unittests/libtests/topology/Makefile.am
short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py
Log:
Added check of cell material ids in mesh (and added corresponding unit tests).
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2008-03-30 02:54:05 UTC (rev 11635)
+++ short/3D/PyLith/trunk/TODO 2008-03-30 03:25:44 UTC (rev 11636)
@@ -4,16 +4,8 @@
Release 1.1
- 1. Fix parallel output for VTK files.
- a. Cell fields where proc 0 doesn't have cells.
- b. Faults are missing cells.
-
2. Fix calcTractionsChange.
- 3. Add more error checking.
-
- a. Add check to make sure every material in mesh has a material model.
-
Benchmarks
1. Strike-slip benchmark.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2008-03-30 02:54:05 UTC (rev 11635)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2008-03-30 03:25:44 UTC (rev 11636)
@@ -91,6 +91,7 @@
meshio/VertexFilterVecNorm.cc \
topology/FieldsManager.cc \
topology/Distributor.cc \
+ topology/MeshOps.cc \
utils/EventLogger.cc
libpylith_la_LDFLAGS = $(AM_LDFLAGS) $(PYTHON_LA_LDFLAGS)
Modified: short/3D/PyLith/trunk/libsrc/topology/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2008-03-30 02:54:05 UTC (rev 11635)
+++ short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2008-03-30 03:25:44 UTC (rev 11636)
@@ -14,8 +14,9 @@
include $(top_srcdir)/subpackage.am
subpkginclude_HEADERS = \
+ Distributor.hh \
FieldsManager.hh \
- Distributor.hh
+ MeshOps.hh
noinst_HEADERS =
Added: short/3D/PyLith/trunk/libsrc/topology/MeshOps.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshOps.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshOps.cc 2008-03-30 03:25:44 UTC (rev 11636)
@@ -0,0 +1,56 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "MeshOps.hh" // implementation of class methods
+
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
+#include <assert.h> // USES assert()
+
+#include <algorithm> // USES std::sort, std::find
+
+// ----------------------------------------------------------------------
+void
+pylith::topology::MeshOps::checkMaterialIds(const ALE::Obj<Mesh>& mesh,
+ int* materialIds,
+ const int numMaterials)
+{ // checkMaterialIds
+ assert( (0 == numMaterials && 0 == materialIds) ||
+ (0 < numMaterials && 0 != materialIds) );
+
+ const ALE::Obj<ALE::Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const Mesh::label_sequence::iterator cellsEnd = cells->end();
+ const ALE::Obj<Mesh::label_type>& materialsLabel = mesh->getLabel("material-id");
+
+ int* matBegin = materialIds;
+ int* matEnd = materialIds + numMaterials;
+ std::sort(matBegin, matEnd);
+
+ for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ c_iter != cellsEnd;
+ ++c_iter) {
+ const int cellId = mesh->getValue(materialsLabel, *c_iter);
+ const int* result = std::find(matBegin, matEnd, cellId);
+ if (result == matEnd) {
+ std::ostringstream msg;
+ msg << "Material id '" << cellId << "' for cell '" << *c_iter
+ << "' does not match the id of any available materials or interfaces.";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // for
+} // checkMaterialIds
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/MeshOps.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshOps.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshOps.hh 2008-03-30 03:25:44 UTC (rev 11636)
@@ -0,0 +1,68 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/topology/MeshOps.hh
+ *
+ * @brief Temporary object for doing operations on a PETSc
+ * mesh. Object will be replaced by a PyLith Mesh object that inherits
+ * or templates over the PETSc mesh object.
+ */
+
+#if !defined(pylith_topology_meshops_hh)
+#define pylith_topology_meshops_hh
+
+#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+
+namespace pylith {
+ namespace topology {
+ class MeshOps;
+ class TestMeshOps;
+ } // topology
+} // pylith
+
+class pylith::topology::MeshOps
+{ // MeshOps
+ friend class TestMeshOps; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Check to make sure material id of every cell matches the id of
+ * one of the materials.
+ *
+ * @param mesh PETSc mesh.
+ * @param materialIds Array of ids for all materials and cohesive cell interfaces.
+ * @param numMaterials Size of array.
+ */
+ static
+ void checkMaterialIds(const ALE::Obj<Mesh>& mesh,
+ int* materialIds,
+ const int numMaterials);
+
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ MeshOps(const MeshOps&);
+
+ /// Not implemented
+ const MeshOps& operator=(const MeshOps&);
+
+
+}; // MeshOps
+
+#endif // pylith_topology_meshops_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src 2008-03-30 02:54:05 UTC (rev 11635)
+++ short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src 2008-03-30 03:25:44 UTC (rev 11636)
@@ -13,6 +13,7 @@
#header{
#include "pylith/topology/FieldsManager.hh"
#include "pylith/topology/Distributor.hh"
+#include "pylith/topology/MeshOps.hh"
#include "pylith/utils/sievetypes.hh"
#include "pylith/utils/petscfwd.h"
#include <petscmesh.hh>
@@ -286,7 +287,41 @@
Mesh_view(self.thisptr)
return
+
+ def checkMaterialIds(self, materialIds):
+ """
+ Make sure material id for each cell matches id of a material.
+ """
+ # create shim for checkMaterialIds
+ #embed{ void* Mesh_checkMaterialIds(void* objVptr, int* materialIds, int numMaterials)
+ try {
+ ALE::Obj<ALE::Mesh>* mesh = (ALE::Obj<ALE::Mesh>*) objVptr;
+ assert(0 != mesh);
+ assert(!mesh->isNull());
+ pylith::topology::MeshOps::checkMaterialIds(*mesh, materialIds, numMaterials);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ cdef int* materialIdVals
+ materialIdVals = NULL
+ numMaterials = len(materialIds);
+ if numMaterials > 0:
+ materialIdVals = <int*> malloc(numMaterials*sizeof(int));
+ for i from 0 <= i < numMaterials:
+ materialIdVals[i] = materialIds[i]
+ Mesh_checkMaterialIds(self.thisptr, materialIdVals, numMaterials)
+ free(materialIdVals)
+ return
+
def _createHandle(self):
"""
Wrap pointer to C++ object in PyCObject.
Modified: short/3D/PyLith/trunk/pylith/problems/Problem.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Problem.py 2008-03-30 02:54:05 UTC (rev 11635)
+++ short/3D/PyLith/trunk/pylith/problems/Problem.py 2008-03-30 03:25:44 UTC (rev 11636)
@@ -133,6 +133,7 @@
"for spatial dimension '%d'." % \
(self.dimension, mesh.dimension)
+ # Check to make sure ids of materials and interfaces are unique
materialIds = {}
for material in self.materials.components():
if material.quadrature.spaceDim != self.dimension:
@@ -155,6 +156,9 @@
(materialIds[interface.id], interface.label, interface.id)
materialIds[interface.id] = interface.label
+ # Check to make sure material-id for each cell matches the id of a material
+ self.mesh.checkMaterialIds(materialIds.keys())
+
self._logger.eventEnd(logEvent)
return
Modified: short/3D/PyLith/trunk/pylith/topology/Mesh.py
===================================================================
--- short/3D/PyLith/trunk/pylith/topology/Mesh.py 2008-03-30 02:54:05 UTC (rev 11635)
+++ short/3D/PyLith/trunk/pylith/topology/Mesh.py 2008-03-30 03:25:44 UTC (rev 11636)
@@ -119,6 +119,15 @@
return self.cppHandle.createMatrix(field)
+ def checkMaterialIds(self, materialIds):
+ """
+ Make sure material id for each cell matches id of a material.
+ """
+ assert(None != self.cppHandle)
+ self.cppHandle.checkMaterialIds(materialIds)
+ return
+
+
# PRIVATE METHODS ////////////////////////////////////////////////////
def _configure(self):
Modified: short/3D/PyLith/trunk/tests/2d/bar_quad4/dislocation_dyn.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/bar_quad4/dislocation_dyn.cfg 2008-03-30 02:54:05 UTC (rev 11635)
+++ short/3D/PyLith/trunk/tests/2d/bar_quad4/dislocation_dyn.cfg 2008-03-30 03:25:44 UTC (rev 11636)
@@ -4,17 +4,16 @@
# ----------------------------------------------------------------------
# journal
# ----------------------------------------------------------------------
-#[pylithapp.journal.info]
-#timedependent = 1
-#explicit = 1
-#petsc = 1
-#solverlinear = 1
-#meshioascii = 1
-#homogeneous = 1
-#explicitelasticity = 1
-#quadrature2d = 1
-#fiatlagrange = 1
-#faultcohesivekin = 1
+[pylithapp.journal.info]
+timedependent = 1
+explicit = 1
+petsc = 1
+solverlinear = 1
+homogeneous = 1
+explicitelasticity = 1
+quadrature2d = 1
+fiatlagrange = 1
+faultcohesivekin = 1
# ----------------------------------------------------------------------
# mesh_generator
@@ -35,9 +34,9 @@
dimension = 2
formulation = pylith.problems.Explicit
#formulation.solver.initial_guess_zero = False
-bc = pylith.bc.BCTwoSides
+bc = [pos,neg]
bc.pos = pylith.bc.AbsorbingDampers
-interfaces = pylith.faults.SingleFault
+interfaces = [fault]
# ----------------------------------------------------------------------
# materials
@@ -58,7 +57,6 @@
# boundary conditions
# ----------------------------------------------------------------------
[pylithapp.timedependent.bc.pos]
-id = 10
label = end points
db.label = Absorbing BC
db.iohandler.filename = matprops.spatialdb
@@ -69,7 +67,6 @@
[pylithapp.timedependent.bc.neg]
fixed_dof = [0]
-id = 11
label = sides
db.label = Dirichlet BC
@@ -116,10 +113,15 @@
# output
# ----------------------------------------------------------------------
[pylithapp.problem.formulation.output.output.writer]
-filename = output/dislocation-dyn.vtk
+filename = output/dislocation-soln.vtk
time_format = %05.3f
# Give basename for vtk fault output.
[pylithapp.timedependent.interfaces.fault.output.writer]
filename = output/dislocation-fault.vtk
time_format = %05.3f
+
+[pylithapp.timedependent.materials.material.output]
+cell_filter = pylith.meshio.CellFilterAvg
+writer.filename = dislocation-statevars.vtk
+writer.time_format = %05.3f
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/Makefile.am 2008-03-30 02:54:05 UTC (rev 11635)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/Makefile.am 2008-03-30 03:25:44 UTC (rev 11636)
@@ -22,6 +22,7 @@
# Primary source files
testtopology_SOURCES = \
TestFieldsManager.cc \
+ TestMeshOps.cc \
test_topology.cc
noinst_HEADERS =
Added: short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.cc (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.cc 2008-03-30 03:25:44 UTC (rev 11636)
@@ -0,0 +1,56 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestMeshOps.hh" // Implementation of class methods
+
+#include "pylith/topology/MeshOps.hh" // USES MeshOps
+
+#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
+#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::topology::TestMeshOps );
+
+// ----------------------------------------------------------------------
+// Test checkMaterialIds().
+void
+pylith::topology::TestMeshOps::testCheckMaterialIds(void)
+{ // testCheckMaterialIds
+ ALE::Obj<Mesh> mesh;
+
+ meshio::MeshIOAscii iohandler;
+ iohandler.filename("data/tri3.mesh");
+ iohandler.read(&mesh);
+ CPPUNIT_ASSERT(!mesh.isNull());
+
+ const int numMaterials = 2;
+ int materialIds[numMaterials];
+ materialIds[0] = 4;
+ materialIds[1] = 3;
+
+ MeshOps::checkMaterialIds(mesh, materialIds, numMaterials);
+
+ bool caughtError = false;
+ try {
+ materialIds[0] = 99;
+
+ MeshOps::checkMaterialIds(mesh, materialIds, numMaterials);
+ } catch (const std::runtime_error& err) {
+ caughtError = true;
+ } // try/catch
+ CPPUNIT_ASSERT(caughtError);
+} // testCheckMaterialIds
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.hh (rev 0)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestMeshOps.hh 2008-03-30 03:25:44 UTC (rev 11636)
@@ -0,0 +1,55 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/topology/TestMeshOps.hh
+ *
+ * @brief C++ TestMeshOps object.
+ *
+ * C++ unit testing for MeshOps.
+ */
+
+#if !defined(pylith_topology_testmeshops_hh)
+#define pylith_topology_testmeshops_hh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace topology {
+ class TestMeshOps;
+ } // topology
+} // pylith
+
+/// C++ unit testing for MeshOps.
+class pylith::topology::TestMeshOps : public CppUnit::TestFixture
+{ // class TestMeshOps
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestMeshOps );
+
+ CPPUNIT_TEST( testCheckMaterialIds );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Test checkMaterialIds().
+ void testCheckMaterialIds(void);
+
+}; // class TestMeshOps
+
+#endif // pylith_topology_meshops_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py 2008-03-30 02:54:05 UTC (rev 11635)
+++ short/3D/PyLith/trunk/unittests/pytests/topology/TestMesh.py 2008-03-30 03:25:44 UTC (rev 11636)
@@ -120,6 +120,23 @@
return
+ def test_checkMaterialIds(self):
+ """
+ Test createRealSection().
+ """
+ mesh = self._getMesh()
+ materialIds = [4, 3]
+ mesh.checkMaterialIds(materialIds)
+
+ materialIds[0] = -2
+ caughtError = False
+ try:
+ mesh.checkMaterialIds(materialIds)
+ except RuntimeError:
+ caughtError = True
+ self.assertTrue(caughtError)
+ return
+
# PRIVATE METHODS ////////////////////////////////////////////////////
def _getMesh(self):
More information about the cig-commits
mailing list