[cig-commits] r15533 - in short/3D/PyLith/trunk/libsrc: . topology
brad at geodynamics.org
brad at geodynamics.org
Tue Aug 11 14:56:54 PDT 2009
Author: brad
Date: 2009-08-11 14:56:53 -0700 (Tue, 11 Aug 2009)
New Revision: 15533
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/topology/Distributor.hh
short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc
short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh
short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc
short/3D/PyLith/trunk/libsrc/topology/RefineUniform.hh
Log:
Updated uniform global refinement to use current Mesh object.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2009-08-11 18:00:01 UTC (rev 15532)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2009-08-11 21:56:53 UTC (rev 15533)
@@ -107,11 +107,11 @@
topology/SubMesh.cc \
topology/SolutionFields.cc \
topology/Distributor.cc \
+ topology/MeshRefiner.cc \
+ topology/RefineUniform.cc \
utils/EventLogger.cc \
utils/TestArray.cc
-# topology/MeshRefiner.cc \
-# topology/RefineUniform.cc
libpylith_la_LDFLAGS = $(AM_LDFLAGS) $(PYTHON_LA_LDFLAGS)
libpylith_la_LIBADD = \
@@ -133,4 +133,5 @@
libpylith_la_LIBADD += -lnetcdf_c++ -lnetcdf
endif
+
# End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/Distributor.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Distributor.hh 2009-08-11 18:00:01 UTC (rev 15532)
+++ short/3D/PyLith/trunk/libsrc/topology/Distributor.hh 2009-08-11 21:56:53 UTC (rev 15533)
@@ -24,6 +24,7 @@
#include "pylith/meshio/meshiofwd.hh" // USES DataWriter<Mesh>
+// Distributor ----------------------------------------------------------
class pylith::topology::Distributor
{ // Distributor
friend class TestDistributor; // unit testing
Modified: short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc 2009-08-11 18:00:01 UTC (rev 15532)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc 2009-08-11 21:56:53 UTC (rev 15533)
@@ -14,15 +14,6 @@
#include "MeshRefiner.hh" // implementation of class methods
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
-#include "pylith/meshio/DataWriter.hh" // USES DataWriter
-
-#include <cstring> // USES strlen()
-#include <strings.h> // USES strcasecmp()
-#include <stdexcept> // USES std::runtime_error
-#include <sstream> // USES std::ostringstream
-#include <cassert> // USES assert()
-
// ----------------------------------------------------------------------
// Constructor
pylith::topology::MeshRefiner::MeshRefiner(void)
@@ -35,22 +26,5 @@
{ // destructor
} // destructor
-// ----------------------------------------------------------------------
-// Write refined mesh.
-void
-pylith::topology::MeshRefiner::write(meshio::DataWriter* const writer,
- const ALE::Obj<Mesh>& mesh,
- const spatialdata::geocoords::CoordSys* cs)
-{ // write
- assert(!mesh.isNull());
- const double t = 0.0;
- const int numTimeSteps = 0;
- writer->open(mesh, cs, numTimeSteps);
- writer->openTimeStep(t, mesh, cs);
- writer->closeTimeStep();
- writer->close();
-} // write
-
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh 2009-08-11 18:00:01 UTC (rev 15532)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh 2009-08-11 21:56:53 UTC (rev 15533)
@@ -19,25 +19,10 @@
#if !defined(pylith_topology_meshrefiner_hh)
#define pylith_topology_meshrefiner_hh
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+// Include directives ---------------------------------------------------
+#include "topologyfwd.hh" // forward declarations
-namespace pylith {
- namespace topology {
- class MeshRefiner;
- class TestMeshRefiner;
- } // topology
-
- namespace meshio {
- class DataWriter;
- } // meshio
-} // pylith
-
-namespace spatialdata {
- namespace geocoords {
- class CoordSys;
- } // geocoords
-} // spatialdata
-
+// MeshRefiner ----------------------------------------------------------
class pylith::topology::MeshRefiner
{ // MeshRefiner
friend class TestMeshRefiner; // unit testing
@@ -49,18 +34,21 @@
MeshRefiner(void);
/// Destructor
+ virtual
~MeshRefiner(void);
- /** Write refined mesh.
+ /** Refine mesh.
*
- * @param writer Data writer for partition information.
- * @param mesh Refined mesh.
- * @param cs Coordinate system for mesh.
+ * @param newMesh Refined mesh (result).
+ * @param mesh Mesh to refine.
+ * @param levels Number of levels to refine.
+ * @param fields Solution fields.
*/
- static
- void write(meshio::DataWriter* const writer,
- const ALE::Obj<Mesh>& mesh,
- const spatialdata::geocoords::CoordSys* cs);
+ virtual
+ void refine(Mesh* const newMesh,
+ const Mesh& mesh,
+ const int levels =1,
+ const SolutionFields* fields =0) = 0;
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc 2009-08-11 18:00:01 UTC (rev 15532)
+++ short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc 2009-08-11 21:56:53 UTC (rev 15533)
@@ -14,12 +14,16 @@
#include "RefineUniform.hh" // implementation of class methods
-#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh
+#include "Mesh.hh" // USES Mesh
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
#include <cassert> // USES assert()
+// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+
+// ----------------------------------------------------------------------
template<typename Point>
class Edge : public std::pair<Point, Point> {
public:
@@ -48,100 +52,139 @@
// ----------------------------------------------------------------------
// Refine mesh.
void
-pylith::topology::RefineUniform::refine(ALE::Obj<Mesh>* const newMesh,
- const ALE::Obj<Mesh>& mesh,
- const int levels)
+pylith::topology::RefineUniform::refine(Mesh* const newMesh,
+ const Mesh& mesh,
+ const int levels,
+ const SolutionFields* fields)
{ // refine
- assert(!mesh.isNull());
+ assert(0 != newMesh);
- typedef Mesh::point_type point_type;
- typedef Edge<point_type> edge_type;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
- const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->heightStratum(0);
+ assert(!cells.isNull());
+
+ newMesh->debug(mesh.debug());
+
+ // Assume number of corners per cell is the same for the entire mesh
assert(cells->size() > 0);
- // Assume number of corners per cell is the same for the entire mesh
- const int cellNumCorners = mesh->getNumCellCorners(*cells->begin());
+ const int cellNumCorners = sieveMesh->getNumCellCorners(*cells->begin());
- if (3 == mesh->getDimension() && 4 == cellNumCorners) {
- *newMesh = new Mesh(mesh->comm(), mesh->getDimension(), mesh->debug());
- Obj<Mesh::sieve_type> newSieve =
- new Mesh::sieve_type(mesh->comm(), mesh->debug());
- std::map<edge_type, point_type> edge2vertex;
+ if (3 == mesh.dimension() && 4 == cellNumCorners)
+ _refineTet4(newMesh, mesh, levels);
+
+ // TODO: Add other refinement cases here
+
+ else {
+ std::ostringstream msg;
+ msg << "Unknown case for uniform global refinement.\n"
+ << "mesh dimension: " << mesh.dimension()
+ << ", number of corners in cell: " << cellNumCorners;
+ throw std::runtime_error(msg.str());
+ } // else
+} // refine
- (*newMesh)->setSieve(newSieve);
- ALE::MeshBuilder<Mesh>::refineTetrahedra(*mesh, *(*newMesh), edge2vertex);
- // Fix material ids
- const int numCells = mesh->heightStratum(0)->size();
- const ALE::Obj<Mesh::label_type>& materials =
- mesh->getLabel("material-id");
- const ALE::Obj<Mesh::label_type>& newMaterials =
- (*newMesh)->createLabel("material-id");
+// ----------------------------------------------------------------------
+// Refine tet4 mesh.
+void
+pylith::topology::RefineUniform::_refineTet4(Mesh* const newMesh,
+ const Mesh& mesh,
+ const int levels)
+{ // _refineTet4
+ assert(0 != newMesh);
- for(int c = 0; c < numCells; ++c) {
- const int material = mesh->getValue(materials, c);
+ typedef SieveMesh::point_type point_type;
+ typedef Edge<point_type> edge_type;
- for(int i = 0; i < 8; ++i)
- (*newMesh)->setValue(newMaterials, c*8+i, material);
- } // for
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh>& newSieveMesh = newMesh->sieveMesh();
+ assert(!newSieveMesh.isNull());
- // Fix groups, assuming vertex groups
- const int numNewVertices = (*newMesh)->depthStratum(0)->size();
- const int numNewCells = (*newMesh)->heightStratum(0)->size();
- const ALE::Obj<std::set<std::string> >& sectionNames =
- mesh->getIntSections();
+ ALE::Obj<SieveMesh::sieve_type> newSieve =
+ new SieveMesh::sieve_type(mesh.comm(), mesh.debug());
- for(std::set<std::string>::const_iterator name=sectionNames->begin();
- name != sectionNames->end();
- ++name) {
- const ALE::Obj<Mesh::int_section_type>& group = mesh->getIntSection(*name);
- const ALE::Obj<Mesh::int_section_type>& newGroup =
- (*newMesh)->getIntSection(*name);
- const Mesh::int_section_type::chart_type& chart = group->getChart();
+ std::map<edge_type, point_type> edge2vertex;
+
+ newSieveMesh->setSieve(newSieve);
+ ALE::MeshBuilder<Mesh>::refineTetrahedra(*mesh.sieveMesh(), *newSieveMesh,
+ edge2vertex);
+
+ // Fix material ids
+ const int numCells = sieveMesh->heightStratum(0)->size();
+ const ALE::Obj<SieveMesh::label_type>& materials =
+ sieveMesh->getLabel("material-id");
+ const ALE::Obj<SieveMesh::label_type>& newMaterials =
+ newSieveMesh->createLabel("material-id");
+
+ const int numNewCellsPerCell = 8; // :KLUDGE: depends on levels
+ for(int icell=0; icell < numCells; ++icell) {
+ const int material = sieveMesh->getValue(materials, icell);
+
+ for(int i=0; i < numNewCellsPerCell; ++i)
+ newSieveMesh->setValue(newMaterials, icell*8+i, material);
+ } // for
+
+ // Fix groups, assuming vertex groups
+ const int numNewVertices = newSieveMesh->depthStratum(0)->size();
+ const int numNewCells = newSieveMesh->heightStratum(0)->size();
+ const ALE::Obj<std::set<std::string> >& sectionNames =
+ sieveMesh->getIntSections();
+
+ const std::set<std::string>::const_iterator namesBegin =
+ sectionNames->begin();
+ const std::set<std::string>::const_iterator namesEnd =
+ sectionNames->end();
+ for (std::set<std::string>::const_iterator name=namesBegin;
+ name != namesEnd;
+ ++name) {
+ const ALE::Obj<Mesh::IntSection>& group = sieveMesh->getIntSection(*name);
+ const ALE::Obj<Mesh::IntSection>& newGroup =
+ newSieveMesh->getIntSection(*name);
+ const Mesh::IntSection::chart_type& chart = group->getChart();
- newGroup->setChart(Mesh::int_section_type::chart_type(numNewCells,
- numNewCells + numNewVertices));
- const Mesh::int_section_type::chart_type& newChart = newGroup->getChart();
+ newGroup->setChart(Mesh::IntSection::chart_type(numNewCells,
+ numNewCells + numNewVertices));
+ const Mesh::IntSection::chart_type& newChart = newGroup->getChart();
-
- const int chartMax = chart.max();
- for(int p = chart.min(), pNew = newChart.min(); p < chartMax; ++p, ++pNew) {
- if (group->getFiberDimension(p))
- newGroup->setFiberDimension(pNew, 1);
- } // for
- const std::map<edge_type, point_type>::const_iterator edge2VertexEnd =
- edge2vertex.end();
- for(std::map<edge_type, point_type>::const_iterator e_iter=edge2vertex.begin();
- e_iter != edge2VertexEnd;
- ++e_iter) {
- const point_type vertexA = e_iter->first.first;
- const point_type vertexB = e_iter->first.second;
+ const int chartMax = chart.max();
+ for (int p = chart.min(), pNew = newChart.min(); p < chartMax; ++p, ++pNew) {
+ if (group->getFiberDimension(p))
+ newGroup->setFiberDimension(pNew, 1);
+ } // for
+ const std::map<edge_type, point_type>::const_iterator edge2VertexEnd =
+ edge2vertex.end();
+ for (std::map<edge_type, point_type>::const_iterator e_iter=edge2vertex.begin();
+ e_iter != edge2VertexEnd;
+ ++e_iter) {
+ const point_type vertexA = e_iter->first.first;
+ const point_type vertexB = e_iter->first.second;
- if (group->getFiberDimension(vertexA) && group->getFiberDimension(vertexB))
- if (group->restrictPoint(vertexA)[0] == group->restrictPoint(vertexB)[0])
- newGroup->setFiberDimension(e_iter->second, 1);
- } // for
+ if (group->getFiberDimension(vertexA) && group->getFiberDimension(vertexB))
+ if (group->restrictPoint(vertexA)[0] == group->restrictPoint(vertexB)[0])
+ newGroup->setFiberDimension(e_iter->second, 1);
+ } // for
- newGroup->allocatePoint();
- for(int p = chart.min(), pNew = newChart.min(); p < chartMax; ++p, ++pNew) {
- if (group->getFiberDimension(p))
- newGroup->updatePoint(pNew, group->restrictPoint(p));
- }
- for(std::map<edge_type, point_type>::const_iterator e_iter=edge2vertex.begin();
- e_iter != edge2VertexEnd;
- ++e_iter) {
- const point_type vertexA = e_iter->first.first;
- const point_type vertexB = e_iter->first.second;
+ newGroup->allocatePoint();
+ for (int p=chart.min(), pNew = newChart.min(); p < chartMax; ++p, ++pNew) {
+ if (group->getFiberDimension(p))
+ newGroup->updatePoint(pNew, group->restrictPoint(p));
+ } // for
+ for (std::map<edge_type, point_type>::const_iterator e_iter=edge2vertex.begin();
+ e_iter != edge2VertexEnd;
+ ++e_iter) {
+ const point_type vertexA = e_iter->first.first;
+ const point_type vertexB = e_iter->first.second;
- if (group->getFiberDimension(vertexA) && group->getFiberDimension(vertexB))
- if (group->restrictPoint(vertexA)[0] == group->restrictPoint(vertexB)[0])
- newGroup->updatePoint(e_iter->second, group->restrictPoint(vertexA));
- } // for
+ if (group->getFiberDimension(vertexA) && group->getFiberDimension(vertexB))
+ if (group->restrictPoint(vertexA)[0] == group->restrictPoint(vertexB)[0])
+ newGroup->updatePoint(e_iter->second, group->restrictPoint(vertexA));
} // for
- } else {
- *newMesh = mesh;
- } // if/else
-} // refine
+ } // for
+} // _refineTet4
// End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/RefineUniform.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/RefineUniform.hh 2009-08-11 18:00:01 UTC (rev 15532)
+++ short/3D/PyLith/trunk/libsrc/topology/RefineUniform.hh 2009-08-11 21:56:53 UTC (rev 15533)
@@ -19,15 +19,10 @@
#if !defined(pylith_topology_refineuniform_hh)
#define pylith_topology_refineuniform_hh
+// Include directives ---------------------------------------------------
#include "MeshRefiner.hh" // ISA MeshRefiner
-namespace pylith {
- namespace topology {
- class RefineUniform;
- class TestRefineUniform;
- } // topology
-} // pylith
-
+// RefineUniform --------------------------------------------------------
class pylith::topology::RefineUniform : public MeshRefiner
{ // RefineUniform
friend class TestRefineUniform; // unit testing
@@ -46,11 +41,26 @@
* @param newMesh Refined mesh (result).
* @param mesh Mesh to refine.
* @param levels Number of levels to refine.
+ * @param fields Solution fields.
*/
- void refine(ALE::Obj<Mesh>* const newMesh,
- const ALE::Obj<Mesh>& mesh,
- const int levels);
+ void refine(Mesh* const newMesh,
+ const Mesh& mesh,
+ const int levels =1,
+ const SolutionFields* fields =0);
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ /** Refine tet4 mesh.
+ *
+ * @param newMesh Refined mesh (result).
+ * @param mesh Mesh to refine.
+ * @param levels Number of levels to refine.
+ */
+ void _refineTet4(Mesh* const newMesh,
+ const Mesh& mesh,
+ const int levels =1);
+
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
More information about the CIG-COMMITS
mailing list