[cig-commits] r14614 - in short/3D/PyLith/branches/pylith-swig: . libsrc libsrc/faults unittests/libtests/bc

brad at geodynamics.org brad at geodynamics.org
Tue Apr 7 14:45:08 PDT 2009


Author: brad
Date: 2009-04-07 14:45:08 -0700 (Tue, 07 Apr 2009)
New Revision: 14614

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/FaultCohesiveKin.cc
   short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/Makefile.am
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/TestAbsorbingDampers.cc
   short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/TestBoundaryMesh.cc
Log:
Updated FaultCohesiveKin interfaces to allow use of it in other unit tests. Updated boundary mesh tests.

Modified: short/3D/PyLith/branches/pylith-swig/TODO
===================================================================
--- short/3D/PyLith/branches/pylith-swig/TODO	2009-04-07 15:06:44 UTC (rev 14613)
+++ short/3D/PyLith/branches/pylith-swig/TODO	2009-04-07 21:45:08 UTC (rev 14614)
@@ -56,8 +56,6 @@
 
   libtests/bc/TestDirichletBoundary::testVertexField()
   libtests/bc/TestDirichletBoundary::testBoundaryMesh()
-  libtests/bc/TestBoundaryMesh
-    requires faults to be working
 
   pytests/bc/TestDirichletBoundary
   pytests/bc/TestNeumann (output)

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am	2009-04-07 15:06:44 UTC (rev 14613)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am	2009-04-07 21:45:08 UTC (rev 14614)
@@ -35,10 +35,11 @@
 	faults/BruneSlipFn.cc \
 	faults/LiuCosSlipFn.cc \
 	faults/EqKinSrc.cc \
+	faults/TopologyOps.cc \
+	faults/CohesiveTopology.cc \
 	faults/FaultCohesive.cc \
 	faults/FaultCohesiveDyn.cc \
-	faults/TopologyOps.cc \
-	faults/CohesiveTopology.cc \
+	faults/FaultCohesiveKin.cc \
 	feassemble/CellGeometry.cc \
 	feassemble/Constraint.cc \
 	feassemble/GeometryPoint1D.cc \
@@ -101,7 +102,6 @@
 	utils/EventLogger.cc \
 	utils/TestArray.cc
 
-#	faults/FaultCohesiveKin.cc \
 # 	materials/GenMaxwellIsotropic3D.cc \
 # 	topology/Distributor.cc \
 # 	topology/MeshRefiner.cc \

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc	2009-04-07 15:06:44 UTC (rev 14613)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc	2009-04-07 21:45:08 UTC (rev 14614)
@@ -19,18 +19,17 @@
 
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
-#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
-#include "pylith/topology/FieldOps.hh" // USES FieldOps
-#include "pylith/utils/array.hh" // USES double_array
-#include <petscmat.h> // USES PETSc Mat
 
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+
 #include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
-#include "spatialdata/spatialdb/SpatialDB.hh" // USES CoordSys
+#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
 #include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
 
-#include <Completion.hh> // USES completeSection
-#include <Selection.hh> // Algorithms for submeshes
-
 #include <cmath> // USES pow(), sqrt()
 #include <strings.h> // USES strcasecmp()
 #include <cstring> // USES strlen()
@@ -40,6 +39,11 @@
 #include <stdexcept> // USES std::runtime_error
 
 // ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
+
+// ----------------------------------------------------------------------
 // Default constructor.
 pylith::faults::FaultCohesiveKin::FaultCohesiveKin(void) :
   _fields(0)
@@ -90,48 +94,48 @@
     throw std::runtime_error("Normal direction for fault orientation must be "
 			     "a vector with 3 components.");
 
-  CohesiveTopology::createParallel(&_faultMesh, &_cohesiveToFault, mesh, id(),
-				   _useLagrangeConstraints());
+  CohesiveTopology::createFaultParallel(_faultMesh, &_cohesiveToFault, 
+					mesh, id(), _useLagrangeConstraints());
   //_faultMesh->getLabel("height")->view("Fault mesh height");
   //_faultMesh->view("FAULT MESH");
 
   delete _fields; 
-  _fields = new topology::Fields<topology::Field<topology::SubMesh> >;
+  _fields = new topology::Fields<topology::Field<topology::SubMesh> >(*_faultMesh);
 
-  // Setup pseudo-stiffness of cohesive cells to improve conditioning
-  // of Jacobian matrix
-  _calcConditioning(cs, matDB);
-
-  // Compute orientation at vertices in fault mesh.
-  _calcOrientation(upDir, normalDir);
-
-  // Compute tributary area for each vertex in fault mesh.
-  _calcArea();
-
   const srcs_type::const_iterator srcsEnd = _eqSrcs.end();
   for (srcs_type::iterator s_iter=_eqSrcs.begin(); 
        s_iter != srcsEnd; 
        ++s_iter) {
     EqKinSrc* src = s_iter->second;
     assert(0 != src);
-    src->initialize(_faultMesh, *_normalizer);
+    src->initialize(*_faultMesh, *_normalizer);
   } // for
 
   // Allocate slip field
-  const ALE::Obj<SubMesh::SieveMesh>& faultSieveMesh = _faultMesh->sieveMesh();
-  const ALE::Obj<SubMesh::label_sequence>& vertices =
-    _faultSieveMesh->depthStratum(0);
+  const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+  const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
+    faultSieveMesh->depthStratum(0);
   assert(!vertices.isNull());
-  _fields->add("slip");
+  _fields->add("slip", "slip");
   topology::Field<topology::SubMesh>& slip = _fields->get("slip");
   slip.newSection(vertices, cs->spaceDim());
   slip.allocate();
 
   // Allocate cumulative slip field
-  _fields->add("cumulative slip");
+  _fields->add("cumulative slip", "cumulative_slip");
   topology::Field<topology::SubMesh>& cumSlip = _fields->get("cumulative slip");
   cumSlip.newSection(slip);
   cumSlip.allocate();
+
+  // Setup pseudo-stiffness of cohesive cells to improve conditioning
+  // of Jacobian matrix
+  _calcConditioning(cs, matDB);
+
+  // Compute orientation at vertices in fault mesh.
+  _calcOrientation(upDir, normalDir);
+
+  // Compute tributary area for each vertex in fault mesh.
+  _calcArea();
 } // initialize
 
 // ----------------------------------------------------------------------
@@ -143,6 +147,7 @@
 			     const double t,
 			     topology::SolutionFields* const fields)
 { // integrateResidual
+#if 0
   assert(0 != fields);
 
   // Cohesive cells with normal vertices i and j, and constraint
@@ -276,6 +281,7 @@
 
   // FIX THIS
   PetscLogFlops(cellsCohesiveSize*numConstraintVert*spaceDim*spaceDim*7);
+#endif
 } // integrateResidual
 
 // ----------------------------------------------------------------------
@@ -283,12 +289,11 @@
 // not require assembly across cells, vertices, or processors.
 void
 pylith::faults::FaultCohesiveKin::integrateResidualAssembled(
-				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)
 { // integrateResidualAssembled
+#if 0
   assert(!residual.isNull());
   assert(0 != fields);
   assert(!mesh.isNull());
@@ -344,6 +349,7 @@
       assert(0 != slip);
       residual->updatePoint(vertexMesh, slip);
     } // if
+#endif
 } // integrateResidualAssembled
 
 // ----------------------------------------------------------------------
@@ -351,11 +357,11 @@
 // require assembly across cells, vertices, or processors.
 void
 pylith::faults::FaultCohesiveKin::integrateJacobianAssembled(
-				    PetscMat* mat,
-				    const double t,
-				    topology::FieldsManager* const fields,
-				    const ALE::Obj<Mesh>& mesh)
+				       topology::Jacobian* jacobian,
+				       const double t,
+				       topology::SolutionFields* const fields)
 { // integrateJacobianAssembled
+#if 0
   assert(0 != mat);
   assert(0 != fields);
   assert(!mesh.isNull());
@@ -485,6 +491,7 @@
   } // for
   PetscLogFlops(cellsCohesiveSize*numConstraintVert*spaceDim*spaceDim*4);
   _needNewJacobian = false;
+#endif
 } // integrateJacobianAssembled
   
 // ----------------------------------------------------------------------
@@ -497,8 +504,8 @@
   assert(0 != _fields);
 
   // Update cumulative slip
-  topology::Field<topology::SubMesh>& cumSlip = fields->get("cumulative slip");
-  topology::Field<topology::SubMesh>& slip = fields->get("slip");
+  topology::Field<topology::SubMesh>& cumSlip = _fields->get("cumulative slip");
+  topology::Field<topology::SubMesh>& slip = _fields->get("slip");
   if (!_useSolnIncr)
     cumSlip.zero();
   cumSlip += slip;
@@ -535,15 +542,14 @@
   } // if
 
   const int numCorners = _quadrature->refGeometry().numCorners();
-  const ALE::Obj<Mesh::SieveMesh::label_sequence>& cells = 
+  const ALE::Obj<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();
-  for (Mesh::label_sequence::iterator c_iter=cellsBegin;
+  const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+  for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
        c_iter != cellsEnd;
        ++c_iter) {
-    const int cellNumCorners = mesh->getNumCellCorners(*c_iter);
+    const int cellNumCorners = sieveMesh->getNumCellCorners(*c_iter);
     if (3*numCorners != cellNumCorners) {
       std::ostringstream msg;
       msg << "Number of vertices in reference cell (" << numCorners 
@@ -557,13 +563,12 @@
 
 // ----------------------------------------------------------------------
 // Get vertex field associated with integrator.
-const ALE::Obj<pylith::real_section_type>&
+const pylith::topology::Field<pylith::topology::SubMesh>&
 pylith::faults::FaultCohesiveKin::vertexField(
-				    VectorFieldEnum* fieldType,
-				    const char* name,
-				    const ALE::Obj<Mesh>& mesh,
-				    topology::FieldsManager* fields)
+				  const char* name,
+				  const topology::SolutionFields& fields)
 { // vertexField
+#if 0
   assert(!_faultMesh.isNull());
   assert(!_orientation.isNull());
   assert(0 != _normalizer);
@@ -662,16 +667,17 @@
   } // if
 
   return _bufferTmp;
+#endif
 } // vertexField
 
 // ----------------------------------------------------------------------
 // Get cell field associated with integrator.
-const ALE::Obj<pylith::real_section_type>&
-pylith::faults::FaultCohesiveKin::cellField(VectorFieldEnum* fieldType,
-					    const char* name,
-					    const ALE::Obj<Mesh>& mesh,
-					    topology::FieldsManager* fields)
+const pylith::topology::Field<pylith::topology::SubMesh>&
+pylith::faults::FaultCohesiveKin::cellField(
+				      const char* name,
+				      const topology::SolutionFields& fields)
 { // cellField
+#if 0
   // Should not reach this point if requested field was found
   std::ostringstream msg;
   msg << "Request for unknown cell field '" << name
@@ -680,6 +686,7 @@
 
   // Return generic section to satisfy member function definition.
   //return _outputCellVector;
+#endif
 } // cellField
 
 // ----------------------------------------------------------------------
@@ -688,6 +695,7 @@
 pylith::faults::FaultCohesiveKin::_calcOrientation(const double_array& upDir,
 						   const double_array& normalDir)
 { // _calcOrientation
+#if 0
   assert(!_faultMesh.isNull());
 
   // Get vertices in fault mesh
@@ -844,6 +852,7 @@
   } // if
 
   //_orientation->view("ORIENTATION");
+#endif
 } // _calcOrientation
 
 // ----------------------------------------------------------------------
@@ -852,6 +861,7 @@
 				 const spatialdata::geocoords::CoordSys* cs,
 				 spatialdata::spatialdb::SpatialDB* matDB)
 { // _calcConditioning
+#if 0
   assert(0 != cs);
   assert(0 != matDB);
   assert(!_faultMesh.isNull());
@@ -916,12 +926,14 @@
   PetscLogFlops(count * 2);
 
   matDB->close();
+#endif
 } // _calcConditioning
 
 // ----------------------------------------------------------------------
 void
 pylith::faults::FaultCohesiveKin::_calcArea(void)
 { // _calcArea
+#if 0
   assert(!_faultMesh.isNull());
 
   // Get vertices in fault mesh
@@ -995,6 +1007,7 @@
   _faultMesh->getSendOverlap()->view("Send fault overlap");
   _faultMesh->getRecvOverlap()->view("Receive fault overlap");
 #endif
+#endif
 } // _calcArea
 
 // ----------------------------------------------------------------------
@@ -1002,10 +1015,10 @@
 // NOTE: We must convert vertex labels to fault vertex labels
 void
 pylith::faults::FaultCohesiveKin::_calcTractionsChange(
-				 ALE::Obj<real_section_type>* tractions,
-				 const ALE::Obj<Mesh>& mesh,
-				 const ALE::Obj<real_section_type>& solution)
+			     topology::Field<topology::SubMesh>* tractions,
+			     const topology::Field<topology::Mesh>& solution)
 { // _calcTractionsChange
+#if 0
   assert(0 != tractions);
   assert(!mesh.isNull());
   assert(!solution.isNull());
@@ -1097,6 +1110,7 @@
   _pseudoStiffness->view("CONDITIONING");
   (*tractions)->view("TRACTIONS");
 #endif
+#endif
 } // _calcTractionsChange
 
 

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh	2009-04-07 15:06:44 UTC (rev 14613)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh	2009-04-07 21:45:08 UTC (rev 14614)
@@ -159,7 +159,7 @@
    */
   const topology::Field<topology::SubMesh>&
   cellField(const char* name,
-	      const topology::SolutionFields& fields);
+	    const topology::SolutionFields& fields);
 
   // PROTECTED METHODS //////////////////////////////////////////////////
 protected :
@@ -229,7 +229,8 @@
   topology::Fields<topology::Field<topology::SubMesh> >* _fields;
 
   /// Map label of cohesive cell to label of cells in fault mesh.
-  std::map<Mesh::point_type, Mesh::point_type> _cohesiveToFault;
+  std::map<topology::Mesh::SieveMesh::point_type, 
+	   topology::SubMesh::SieveMesh::point_type> _cohesiveToFault;
 
 }; // class FaultCohesiveKin
 

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/Makefile.am	2009-04-07 15:06:44 UTC (rev 14613)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/Makefile.am	2009-04-07 21:45:08 UTC (rev 14614)
@@ -49,14 +49,14 @@
 	TestAbsorbingDampersQuad4.cc \
 	TestAbsorbingDampersTet4.cc \
 	TestAbsorbingDampersHex8.cc \
+	TestBoundaryMesh.cc \
+	TestBoundaryMeshTri3.cc \
+	TestBoundaryMeshQuad4.cc \
+	TestBoundaryMeshTet4.cc \
+	TestBoundaryMeshHex8.cc \
 	test_bc.cc
 
 
-#	TestBoundaryMesh.cc \
-#	TestBoundaryMeshTri3.cc \
-#	TestBoundaryMeshQuad4.cc \
-#	TestBoundaryMeshTet4.cc \
-#	TestBoundaryMeshHex8.cc
 
 noinst_HEADERS = \
 	TestAbsorbingDampers.hh \

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/TestAbsorbingDampers.cc	2009-04-07 15:06:44 UTC (rev 14613)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/TestAbsorbingDampers.cc	2009-04-07 21:45:08 UTC (rev 14614)
@@ -278,7 +278,7 @@
     iohandler.filename(_data->meshFilename);
     iohandler.read(mesh);
 
-    // Set up coordinates
+    // Set coordinate system
     spatialdata::geocoords::CSCart cs;
     cs.setSpaceDim(mesh->dimension());
     cs.initialize();

Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/TestBoundaryMesh.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/TestBoundaryMesh.cc	2009-04-07 15:06:44 UTC (rev 14613)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/bc/TestBoundaryMesh.cc	2009-04-07 21:45:08 UTC (rev 14614)
@@ -17,12 +17,17 @@
 #include "data/BoundaryMeshData.hh" // USES BoundaryMeshData
 
 #include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
 #include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
 #include "pylith/faults/FaultCohesiveKin.hh" // USES FaultsCohesiveKin
 
-#include <Selection.hh> // USES submesh algorithms
+#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
 
 // ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
+typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
+
+// ----------------------------------------------------------------------
 // Setup testing data.
 void
 pylith::bc::TestBoundaryMesh::setUp(void)
@@ -46,25 +51,24 @@
   CPPUNIT_ASSERT(0 != _data);
 
   topology::Mesh mesh;
-
   meshio::MeshIOAscii iohandler;
   iohandler.filename(_data->filename);
   iohandler.read(&mesh);
 
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
+  // Set coordinate system
+  spatialdata::geocoords::CSCart cs;
+  cs.setSpaceDim(mesh.dimension());
+  cs.initialize();
+  mesh.coordsys(&cs);
 
-  const char* label = _data->bcLabel;
-  const ALE::Obj<SieveSubMesh>& subMesh = 
-    ALE::Selection<Mesh>::submeshV<SieveSubMesh>(sieveMesh, 
-					    sieveMesh->getIntSection(label));
-  CPPUNIT_ASSERT(!subMesh.isNull());
+  // Create submesh
+  topology::SubMesh submesh(mesh, _data->bcLabel);
+  CPPUNIT_ASSERT(!submesh.sieveMesh().isNull());
 
-  //subMesh->view("SUBMESH WITHOUT FAULT");
-
-  const ALE::Obj<SieveSubMesh::label_sequence>& vertices = subMesh->depthStratum(0);
+  // Check vertices
+  const ALE::Obj<SieveSubMesh::label_sequence>& vertices = 
+    submesh.sieveMesh()->depthStratum(0);
   const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
-
   CPPUNIT_ASSERT_EQUAL(_data->numVerticesNoFault, int(vertices->size()));
 
   int ipt = 0;
@@ -73,24 +77,25 @@
        ++v_iter, ++ipt)
     CPPUNIT_ASSERT_EQUAL(_data->verticesNoFault[ipt], *v_iter);
 
-  const ALE::Obj<SieveSubMesh::label_sequence>& cells = subMesh->heightStratum(1);
+  // Check cells
+  const ALE::Obj<SieveSubMesh::label_sequence>& cells = 
+    submesh.sieveMesh()->heightStratum(1);
   const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
-  const ALE::Obj<sieve_type>& sieve = subMesh->getSieve();
+  const ALE::Obj<SieveSubMesh::sieve_type>& sieve = 
+    submesh.sieveMesh()->getSieve();
   assert(!sieve.isNull());
+  CPPUNIT_ASSERT_EQUAL(_data->numCells, int(cells->size()));
 
-  CPPUNIT_ASSERT_EQUAL(_data->numCells, (int) cells->size());
+  ALE::ISieveVisitor::NConeRetriever<SieveSubMesh::sieve_type> ncV(*sieve, (int) pow(sieve->getMaxConeSize(), submesh.sieveMesh()->depth()));
 
-  ALE::ISieveVisitor::NConeRetriever<sieve_type> ncV(*sieve, (int) pow(sieve->getMaxConeSize(), subMesh->depth()));
-
   int icell = 0;
   int index = 0;
   for (SieveSubMesh::label_sequence::iterator c_iter=cells->begin();
        c_iter != cellsEnd;
        ++c_iter, ++icell) {
-    ALE::ISieveTraversal<sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
+    ALE::ISieveTraversal<SieveSubMesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
     const int coneSize = ncV.getSize();
-    const SieveMesh::point_type *cone = ncV.getPoints();
-
+    const SieveSubMesh::point_type *cone = ncV.getPoints();
     CPPUNIT_ASSERT_EQUAL(_data->numCorners, coneSize);
 
     for(int v = 0; v < coneSize; ++v, ++index)
@@ -107,29 +112,30 @@
   CPPUNIT_ASSERT(0 != _data);
 
   topology::Mesh mesh;
-
   meshio::MeshIOAscii iohandler;
   iohandler.filename(_data->filename);
   iohandler.read(&mesh);
 
-  const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
-  CPPUNIT_ASSERT(!sieveMesh.isNull());
+  // Set coordinate system
+  spatialdata::geocoords::CSCart cs;
+  cs.setSpaceDim(mesh.dimension());
+  cs.initialize();
+  mesh.coordsys(&cs);
 
+  // Adjust topology
   faults::FaultCohesiveKin fault;
   fault.label(_data->faultLabel);
   fault.id(_data->faultId);
-  fault.adjustTopology(sieveMesh, _flipFault);
+  fault.adjustTopology(&mesh, _flipFault);
 
-  const char* label = _data->bcLabel;
-  const ALE::Obj<SieveSubMesh>& subMesh = 
-    ALE::Selection<SieveMesh>::submeshV<SieveSubMesh>(sievMesh,
-						      sievMesh->getIntSection(label));
-  CPPUNIT_ASSERT(!subMesh.isNull());
+  // Create submesh
+  topology::SubMesh submesh(mesh, _data->bcLabel);
+  CPPUNIT_ASSERT(!submesh.sieveMesh().isNull());
 
-  //subMesh->view("Submesh for mesh w/fault");
-  const ALE::Obj<SieveSubMesh::label_sequence>& vertices = subMesh->depthStratum(0);
+  // Check vertices
+  const ALE::Obj<SieveSubMesh::label_sequence>& vertices = 
+    submesh.sieveMesh()->depthStratum(0);
   const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
-
   CPPUNIT_ASSERT_EQUAL(_data->numVerticesFault, int(vertices->size()));
 
   int ipt = 0;
@@ -137,27 +143,26 @@
        v_iter != verticesEnd;
        ++v_iter, ++ipt)
     CPPUNIT_ASSERT_EQUAL(_data->verticesFault[ipt], *v_iter);
-    
-  const ALE::Obj<SieveSubMesh::label_sequence>& cells = subMesh->depthStratum(1);
+
+  // Check cells
+  const ALE::Obj<SieveSubMesh::label_sequence>& cells = 
+    submesh.sieveMesh()->heightStratum(1);
   const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
-  const ALE::Obj<sieve_type>& sieve = subMesh->getSieve();
+  const ALE::Obj<SieveSubMesh::sieve_type>& sieve = 
+    submesh.sieveMesh()->getSieve();
   assert(!sieve.isNull());
-  const int depth = 1;
-  typedef ALE::SieveAlg<SieveMesh> SieveAlg;
-
   CPPUNIT_ASSERT_EQUAL(_data->numCells, int(cells->size()));
 
-  ALE::ISieveVisitor::NConeRetriever<sieve_type> ncV(*sieve, (int) pow(sieve->getMaxConeSize(), subMesh->depth()));
+  ALE::ISieveVisitor::NConeRetriever<SieveSubMesh::sieve_type> ncV(*sieve, (int) pow(sieve->getMaxConeSize(), submesh.sieveMesh()->depth()));
 
   int icell = 0;
   int index = 0;
   for (SieveSubMesh::label_sequence::iterator c_iter=cells->begin();
        c_iter != cellsEnd;
        ++c_iter, ++icell) {
-    ALE::ISieveTraversal<sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
+    ALE::ISieveTraversal<SieveSubMesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
     const int coneSize = ncV.getSize();
-    const SieveMesh::point_type *cone = ncV.getPoints();
-
+    const SieveSubMesh::point_type *cone = ncV.getPoints();
     CPPUNIT_ASSERT_EQUAL(_data->numCorners, coneSize);
 
     for(int v = 0; v < coneSize; ++v, ++index)



More information about the CIG-COMMITS mailing list