[cig-commits] r17213 - in short/3D/PyLith/trunk: libsrc libsrc/topology unittests/libtests/topology
brad at geodynamics.org
brad at geodynamics.org
Fri Sep 24 14:33:18 PDT 2010
Author: brad
Date: 2010-09-24 14:33:18 -0700 (Fri, 24 Sep 2010)
New Revision: 17213
Added:
short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc
short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh
short/3D/PyLith/trunk/libsrc/topology/MeshOrder.cc
short/3D/PyLith/trunk/libsrc/topology/MeshOrder.hh
short/3D/PyLith/trunk/libsrc/topology/MeshOrder.icc
short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc
short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/topology/Makefile.am
short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc
short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh
short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.hh
Log:
Work on better implementation of global uniform refinements.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2010-09-24 03:35:36 UTC (rev 17212)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2010-09-24 21:33:18 UTC (rev 17213)
@@ -139,6 +139,13 @@
utils/TestArray.cc
+# TEMPORARY
+libpylith_la_SOURCES += \
+ topology/MeshRefiner.cc \
+ topology/CellRefinerTri3.cc \
+ topology/MeshOrder.cc
+
+
libpylith_la_LDFLAGS = $(AM_LDFLAGS) $(PYTHON_LA_LDFLAGS)
libpylith_la_LIBADD = \
-lspatialdata \
Added: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc 2010-09-24 21:33:18 UTC (rev 17213)
@@ -0,0 +1,320 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "CellRefinerTri3.hh" // implementation of class methods
+
+#include "MeshOrder.hh" // USES MeshOrder
+
+#include <cassert> // USES assert()
+
+#include <iostream> // TEMPORARY
+// ----------------------------------------------------------------------
+// Constructor
+ALE::CellRefinerTri3::CellRefinerTri3(const mesh_type& mesh) :
+ _mesh(mesh)
+{ // constructor
+ assert(2 == mesh.getDimension());
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+ALE::CellRefinerTri3::~CellRefinerTri3(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Get number of refined cells for each original cell.
+int
+ALE::CellRefinerTri3::numNewCells(const point_type cell)
+{ // numNewCells
+ switch (_cellType(cell)) {
+ case TRIANGLE:
+ return 4;
+ case LINE_COHESIVE_LAGRANGE:
+ return 2;
+ default:
+ assert(0);
+ throw ALE::Exception("Unknown cell type.");
+ } // switch
+} // numNewCells
+
+// ----------------------------------------------------------------------
+// Split cell into smaller cells of same type.
+void
+ALE::CellRefinerTri3::splitCell(const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ point_type* curNewVertex)
+{ // splitCell
+ assert(curNewVertex);
+
+ int numEdges = 0;
+ const EdgeType* edges;
+
+ switch (_cellType(cell)) {
+ case TRIANGLE:
+ _edges_TRIANGLE(&edges, &numEdges, cone, coneSize);
+ break;
+ case LINE_COHESIVE_LAGRANGE:
+ _edges_LINE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize);
+ break;
+ default:
+ throw ALE::Exception("Unknown cell type.");
+ } // switch
+
+ for(int iEdge=0; iEdge < numEdges; ++iEdge) {
+ if (_edgeToVertex.find(edges[iEdge]) == _edgeToVertex.end()) {
+ // if vertex does not exist
+ std::cout << "Edge: " << edges[iEdge] << ", new vertex: " << *curNewVertex << std::endl;
+ _edgeToVertex[edges[iEdge]] = *curNewVertex;
+ ++(*curNewVertex);
+ } // if
+ } // for
+} // splitCell
+
+// ----------------------------------------------------------------------
+// Get refined cells.
+void
+ALE::CellRefinerTri3::getNewCells(const point_type** cells,
+ int* numCells,
+ const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ const MeshOrder& orderOldMesh,
+ const MeshOrder& orderNewMesh)
+{ // getNewCells
+ assert(cells);
+ assert(numCells);
+
+ switch (_cellType(cell)) {
+ case TRIANGLE: {
+ const int coneVertexOffset = orderNewMesh.verticesNormal().min() - orderOldMesh.verticesNormal().min();
+ _newCells_TRIANGLE(cells, numCells, cone, coneSize, coneVertexOffset);
+ break;
+ } // TRIANGLE
+ case LINE_COHESIVE_LAGRANGE: {
+ const int coneVertexOffset = orderNewMesh.verticesCensored().min() - orderOldMesh.verticesCensored().min();
+ _newCells_LINE_COHESIVE_LAGRANGE(cells, numCells, cone, coneSize, coneVertexOffset);
+ break;
+ } // LINE_COHESIVE_LAGRANGE
+ default:
+ throw ALE::Exception("Unknown cell type.");
+ } // switch
+} // getNewCells
+
+// ----------------------------------------------------------------------
+// Set coordinates of new vertices.
+void
+ALE::CellRefinerTri3::setCoordsNewVertices(const ALE::Obj<mesh_type::real_section_type>& newCoordsSection,
+ const ALE::Obj<mesh_type::real_section_type>& oldCoordsSection)
+{ // setCoordsNewVertices
+ assert(!newCoordsSection.isNull());
+ assert(!oldCoordsSection.isNull());
+
+ double coordinatesVertex[3];
+
+ assert(_edgeToVertex.size() > 0);
+ const int spaceDim = newCoordsSection->getFiberDimension(_edgeToVertex.begin()->second);
+ assert(spaceDim > 0 && spaceDim <= 3);
+
+ const edge_map_type::const_iterator edgesEnd = _edgeToVertex.end();
+ for (edge_map_type::const_iterator e_iter = _edgeToVertex.begin(); e_iter != edgesEnd; ++e_iter) {
+ const point_type newVertex = e_iter->second;
+ const point_type edgeVertexA = e_iter->first.first;
+ const point_type edgeVertexB = e_iter->first.second;
+
+ assert(spaceDim == oldCoordsSection->getFiberDimension(edgeVertexA));
+ assert(spaceDim == oldCoordsSection->getFiberDimension(edgeVertexB));
+ assert(spaceDim == newCoordsSection->getFiberDimension(newVertex));
+
+ const mesh_type::real_section_type::value_type* coordsA = oldCoordsSection->restrictPoint(edgeVertexA);
+ const mesh_type::real_section_type::value_type* coordsB = oldCoordsSection->restrictPoint(edgeVertexB);
+ for (int i=0; i < spaceDim; ++i)
+ coordinatesVertex[i] = 0.5*(coordsA[i] + coordsB[i]);
+
+ newCoordsSection->updatePoint(newVertex, coordinatesVertex);
+ } // for
+} // setCoordsNewVertices
+
+// ----------------------------------------------------------------------
+// Get cell type.
+ALE::CellRefinerTri3::CellEnum
+ALE::CellRefinerTri3::_cellType(const point_type cell)
+{ // _cellType
+ assert(!_mesh.getSieve().isNull());
+
+ switch (_mesh.getSieve()->getConeSize(cell)) {
+ case 3:
+ return TRIANGLE;
+ case 6:
+ return LINE_COHESIVE_LAGRANGE;
+ case 0: {
+ std::ostringstream msg;
+ std::cerr << "Internal error. Cone size for mesh point " << cell << " is zero. May be a vertex.";
+ assert(0);
+ throw ALE::Exception("Could not determine cell type during uniform global refinement.");
+ } // case 0
+ default : {
+ std::ostringstream msg;
+ std::cerr << "Internal error. Unknown cone size for mesh point " << cell << ". Unknown cell type.";
+ assert(0);
+ throw ALE::Exception("Could not determine cell type during uniform global refinement.");
+ } // default
+ } // switch
+} // _cellType
+
+// ----------------------------------------------------------------------
+// Get edges of triangular cell.
+void
+ALE::CellRefinerTri3::_edges_TRIANGLE(const EdgeType** edges,
+ int* numEdges,
+ const point_type cone[],
+ const int coneSize)
+{ // _edges_TRIANGLE
+ static EdgeType triEdges[3];
+
+ assert(coneSize == 3);
+ triEdges[0] = EdgeType(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
+ triEdges[1] = EdgeType(std::min(cone[1], cone[2]), std::max(cone[1], cone[2]));
+ triEdges[2] = EdgeType(std::min(cone[2], cone[0]), std::max(cone[2], cone[0]));
+ *numEdges = 3;
+ *edges = triEdges;
+} // _edges_TRIANGLE
+
+// ----------------------------------------------------------------------
+// Get edges of line cohesive cell with Lagrange multipler vertices.
+void
+ALE::CellRefinerTri3::_edges_LINE_COHESIVE_LAGRANGE(const EdgeType** edges,
+ int* numEdges,
+ const point_type cone[],
+ const int coneSize)
+{ // _edges_LINE_COHESIVE_LAGRANGE
+ static EdgeType lineEdges[6];
+
+ assert(coneSize == 6);
+ lineEdges[0] = EdgeType(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
+ lineEdges[1] = EdgeType(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]));
+ lineEdges[2] = EdgeType(std::min(cone[4], cone[5]), std::max(cone[4], cone[5]));
+ *numEdges = 3;
+ *edges = lineEdges;
+} // _edges_LINE_COHESIVE_LAGRANGE
+
+// ----------------------------------------------------------------------
+// Get new cells from refinement of a triangular cell.
+void
+ALE::CellRefinerTri3::_newCells_TRIANGLE(const point_type** cells,
+ int *numCells,
+ const point_type cone[],
+ const int coneSize,
+ const int coneVertexOffset)
+{ // _newCells_TRIANGLE
+ const int coneSizeTri3 = 3;
+ const int numEdgesTri3 = 3;
+ const int numNewCells = 4;
+ const int numNewVertices = 3;
+
+ int numEdges = 0;
+ const EdgeType *edges;
+ _edges_TRIANGLE(&edges, &numEdges, cone, coneSize);
+ assert(numEdgesTri3 == numEdges);
+
+ static point_type triCells[numNewCells*coneSizeTri3];
+ point_type newVertices[numNewVertices];
+ for(int iEdge=0, iNewVertex=0; iEdge < numEdgesTri3; ++iEdge) {
+ if (_edgeToVertex.find(edges[iEdge]) == _edgeToVertex.end()) {
+ throw ALE::Exception("Missing edge in refined mesh");
+ } // if
+ newVertices[iNewVertex++] = _edgeToVertex[edges[iEdge]];
+ } // for
+
+ // new cell 0
+ triCells[0*3+0] = cone[0] + coneVertexOffset;
+ triCells[0*3+1] = newVertices[0];
+ triCells[0*3+2] = newVertices[2];
+
+ // new cell 1
+ triCells[1*3+0] = newVertices[0];
+ triCells[1*3+1] = newVertices[1];
+ triCells[1*3+2] = newVertices[2];
+
+ // new cell 2
+ triCells[2*3+0] = cone[1] + coneVertexOffset;
+ triCells[2*3+1] = newVertices[1];
+ triCells[2*3+2] = newVertices[0];
+
+ // new cell 3
+ triCells[3*3+0] = cone[2] + coneVertexOffset;
+ triCells[3*3+1] = newVertices[2];
+ triCells[3*3+2] = newVertices[1];
+
+ *numCells = numNewCells;
+ *cells = triCells;
+} // _newCells_TRIANGLE
+
+// ----------------------------------------------------------------------
+// Get new cells from refinement of a line cohseive cell with Lagrange
+// multiplier vertices.
+void
+ALE::CellRefinerTri3::_newCells_LINE_COHESIVE_LAGRANGE(const point_type** cells,
+ int *numCells,
+ const point_type cone[],
+ const int coneSize,
+ const int coneVertexOffset)
+{ // _newCells_LINE_COHESIVE_LAGRANGE
+ const int coneSizeLine6 = 6;
+ const int numEdgesLine6 = 3;
+ const int numNewCells = 2;
+ const int numNewVertices = 3;
+
+ int numEdges = 0;
+ const EdgeType *edges;
+ _edges_LINE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize);
+ assert(numEdgesLine6 == numEdges);
+
+ static point_type lineCells[numNewCells*coneSizeLine6];
+ point_type newVertices[numNewVertices];
+ for(int iEdge=0, iNewVertex=0; iEdge < numEdgesLine6; ++iEdge) {
+ if (_edgeToVertex.find(edges[iEdge]) == _edgeToVertex.end()) {
+ throw ALE::Exception("Missing edge in refined mesh");
+ } // if
+ newVertices[iNewVertex++] = _edgeToVertex[edges[iEdge]];
+ } // for
+
+ // new cell 0
+ lineCells[0*6+0] = cone[0] + coneVertexOffset;
+ lineCells[0*6+1] = newVertices[0];
+ lineCells[0*6+2] = cone[2] + coneVertexOffset;
+ lineCells[0*6+3] = newVertices[1];
+ lineCells[0*6+4] = cone[4] + coneVertexOffset;
+ lineCells[0*6+5] = newVertices[2];
+
+ // new cell 1
+ lineCells[1*6+0] = newVertices[0];
+ lineCells[1*6+1] = cone[1] + coneVertexOffset;
+ lineCells[1*6+2] = newVertices[1];
+ lineCells[1*6+3] = cone[3] + coneVertexOffset;
+ lineCells[1*6+4] = newVertices[2];
+ lineCells[1*6+5] = cone[5] + coneVertexOffset;
+
+ *numCells = 2;
+ *cells = lineCells;
+} // _newCells_LINE_COHESIVE_LAGRANGE
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh 2010-09-24 21:33:18 UTC (rev 17213)
@@ -0,0 +1,204 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/topology/CellRefinerTri3.hh
+ *
+ * @brief Object for tri3 refinement of cells.
+ */
+
+#if !defined(pylith_topology_cellrefinertri3_hh)
+#define pylith_topology_cellrefinertri3_hh
+
+// Include directives ---------------------------------------------------
+#include "topologyfwd.hh" // forward declarations
+
+#include <list> // USES std::pair
+
+// CellRefinerTri3 ------------------------------------------------------
+/// Object for tri3 refinement of cells.
+class ALE::CellRefinerTri3
+{ // CellRefinerTri3
+ typedef IMesh<> mesh_type;
+ typedef mesh_type::point_type point_type;
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Constructor
+ *
+ * @param mesh Finite-element mesh.
+ */
+ CellRefinerTri3(const mesh_type& mesh);
+
+ /// Destructor
+ ~CellRefinerTri3(void);
+
+ /** Get number of refined cells for each original cell.
+ *
+ * @param cell Original cell.
+ *
+ * @returns Number of refined cells.
+ */
+ int numNewCells(const point_type cell);
+
+ /** Split cell into smaller cells of same type.
+ *
+ * @param cell Original cell.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ * @param curNewVertex Value for next new vertex.
+ */
+ void splitCell(const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ point_type* curNewVertex);
+
+ /** Get refined cells.
+ *
+ * @param cells Vertices in refined cells (refined mesh).
+ * @param numCells Number of refined cells.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ * @param orderOldMesh Order in old mesh.
+ * @param orderNewMesh Order in new mesh.
+ */
+ void getNewCells(const point_type** cells,
+ int* numCells,
+ const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ const MeshOrder& orderOldMesh,
+ const MeshOrder& orderNewMesh);
+
+ /** Set coordinates of new vertices.
+ *
+ * @param newCoordsSection Coordinates of vertices in new mesh.
+ * @param oldCoordsSection Coordinates of vertices in original mesh.
+ */
+ void setCoordsNewVertices(const ALE::Obj<mesh_type::real_section_type>& newCoordsSection,
+ const ALE::Obj<mesh_type::real_section_type>& oldCoordsSection);
+
+// PRIVATE TYPEDEFS /////////////////////////////////////////////////////
+private :
+
+ template<typename Point>
+ class Edge : public std::pair<Point, Point> {
+ public:
+ Edge() : std::pair<Point, Point>() {};
+ Edge(const Point l) : std::pair<Point, Point>(l, l) {};
+ Edge(const Point l, const Point r) : std::pair<Point, Point>(l, r) {};
+ ~Edge() {};
+ friend std::ostream& operator<<(std::ostream& stream, const Edge& edge) {
+ stream << "(" << edge.first << ", " << edge.second << ")";
+ return stream;
+ };
+ };
+
+ typedef Edge<point_type> EdgeType;
+ typedef std::map<EdgeType, point_type> edge_map_type;
+
+// PRIVATE ENUMS ////////////////////////////////////////////////////////
+private :
+
+ enum CellEnum {
+ TRIANGLE, // Normal triangular cell
+ LINE_COHESIVE_LAGRANGE, // Cohesive cell with Lagrange multiplier vertices
+ };
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /** Get cell type.
+ *
+ * @param cell Cell in original mesh.
+ * @returns Cell type.
+ */
+ CellEnum _cellType(const point_type cell);
+
+ /** Get edges of triangular cell.
+ *
+ * @param edges Edges of cell.
+ * @param numEdges Number of edges.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ */
+ void _edges_TRIANGLE(const EdgeType** edges,
+ int* numEdges,
+ const point_type cone[],
+ const int coneSize);
+
+ /** Get edges of line cohesive cell with Lagrange multipler vertices.
+ *
+ * @param edges Edges of cell.
+ * @param numEdges Number of edges.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ */
+ void _edges_LINE_COHESIVE_LAGRANGE(const EdgeType** edges,
+ int* numEdges,
+ const point_type cone[],
+ const int coneSize);
+
+ /** Get new cells from refinement of a triangular cell.
+ *
+ * @param cells Vertices in refined cells (refined mesh).
+ * @param numCells Number of refined cells.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ * @param coneVertexOffset Offset for cone vertices.
+ */
+ void _newCells_TRIANGLE(const point_type** cells,
+ int *numCells,
+ const point_type cone[],
+ const int coneSize,
+ const int coneVertexOffset);
+
+ /** Get new cells from refinement of a line cohseive cell with
+ * Lagrange multiplier vertices.
+ *
+ * @param cells Vertices in refined cells (refined mesh).
+ * @param numCells Number of refined cells.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ * @param coneVertexOffset Offset for cone vertices.
+ */
+ void _newCells_LINE_COHESIVE_LAGRANGE(const point_type** cells,
+ int *numCells,
+ const point_type cone[],
+ const int coneSize,
+ const int coneVertexOffset);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ const mesh_type& _mesh;
+ edge_map_type _edgeToVertex;
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ CellRefinerTri3(const CellRefinerTri3&); ///< Not implemented
+ const CellRefinerTri3& operator=(const CellRefinerTri3&); ///< Not implemented
+
+}; // CellRefinerTri3
+
+#endif // pylith_topology_cellrefinertri3_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2010-09-24 03:35:36 UTC (rev 17212)
+++ short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2010-09-24 21:33:18 UTC (rev 17213)
@@ -46,6 +46,13 @@
noinst_HEADERS =
+# TEMPORARY
+noinst_HEADERS += \
+ MeshRefiner.hh \
+ CellRefinerTri3.hh \
+ MeshOrder.hh
+
+
# export
clean-local: clean-subpkgincludeHEADERS
BUILT_SOURCES = export-subpkgincludeHEADERS
Added: short/3D/PyLith/trunk/libsrc/topology/MeshOrder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshOrder.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshOrder.cc 2010-09-24 21:33:18 UTC (rev 17213)
@@ -0,0 +1,122 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "MeshOrder.hh" // implementation of class methods
+
+#include <cassert> // USES assert()
+
+// ----------------------------------------------------------------------
+// Constructor
+ALE::MeshOrder::MeshOrder(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+ALE::MeshOrder::~MeshOrder(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Determine order from pre-existing mesh.
+void
+ALE::MeshOrder::initialize(const ALE::Obj<mesh_type>& mesh)
+{ // initialize
+ const ALE::Obj<mesh_type::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const int numCells = cells->size();
+ _cellsNormal = ALE::Interval<point_type>(0, numCells);
+
+ const Obj<mesh_type::label_sequence>& vertices = mesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const int numVertices = vertices->size();
+ _verticesNormal = ALE::Interval<point_type>(numCells, numCells+numVertices);
+
+ const int numEntities = numCells + numVertices;
+ _cellsCensored = ALE::Interval<point_type>(numEntities, numEntities);
+ _verticesCensored = ALE::Interval<point_type>(numEntities, numEntities);
+} // initialize
+
+// ----------------------------------------------------------------------
+// Set range for normal cells.
+void
+ALE::MeshOrder::cellsNormal(const point_type min, const point_type max)
+{ // cellsNormal
+ _cellsNormal = ALE::Interval<point_type>(min, max);
+} // cellsNormal
+
+// ----------------------------------------------------------------------
+// Get range for normal cells.
+const ALE::Interval<ALE::MeshOrder::point_type>&
+ALE::MeshOrder::cellsNormal(void) const
+{ // cellsNormal
+ return _cellsNormal;
+} // cellsNormal
+
+// ----------------------------------------------------------------------
+// Set range for normal vertices.
+void
+ALE::MeshOrder::verticesNormal(const point_type min, const point_type max)
+{ // verticesNormal
+ _verticesNormal = ALE::Interval<point_type>(min, max);
+} // verticesNormal
+
+// ----------------------------------------------------------------------
+// Get range for normal vertices.
+const ALE::Interval<ALE::MeshOrder::point_type>&
+ALE::MeshOrder::verticesNormal(void) const
+{ // verticesNormal
+ return _verticesNormal;
+} // verticesNormal
+
+// ----------------------------------------------------------------------
+// Set range for censored cells.
+void
+ALE::MeshOrder::cellsCensored(const point_type min, const point_type max)
+{ // cellsCensored
+ _cellsCensored = ALE::Interval<point_type>(min, max);
+} // cellsCensored
+
+// ----------------------------------------------------------------------
+// Get range for censored cells.
+const ALE::Interval<ALE::MeshOrder::point_type>&
+ALE::MeshOrder::cellsCensored(void) const
+{ // cellsCensored
+ return _cellsCensored;
+} // cellsCensored
+
+// ----------------------------------------------------------------------
+// Set range for censored vertices.
+void
+ALE::MeshOrder::verticesCensored(const point_type min, const point_type max)
+{ // verticesCensored
+ _verticesCensored = ALE::Interval<point_type>(min, max);
+} // verticesCensored
+
+// ----------------------------------------------------------------------
+// Get range for censored vertices.
+const ALE::Interval<ALE::MeshOrder::point_type>&
+ALE::MeshOrder::verticesCensored(void) const
+{ // verticesCensored
+ return _verticesCensored;
+} // verticesCensored
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/MeshOrder.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshOrder.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshOrder.hh 2010-09-24 21:33:18 UTC (rev 17213)
@@ -0,0 +1,129 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/topology/MeshOrder.hh
+ *
+ * @brief Object for managing order of mesh entities.
+ *
+ * Entities are stored in contiguous ranges.
+ */
+
+#if !defined(pylith_topology_meshorder_hh)
+#define pylith_topology_meshorder_hh
+
+// Include directives ---------------------------------------------------
+#include "topologyfwd.hh" // forward declarations
+
+#include <petscmesh.hh> // HASA ALE::IMesh
+
+// MeshOrder ------------------------------------------------------------
+/// Object for managing order of mesh entities.
+class ALE::MeshOrder
+{ // MeshOrder
+ typedef int point_type;
+ typedef ALE::IMesh<> mesh_type;
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ MeshOrder(void);
+
+ /// Destructor
+ ~MeshOrder(void);
+
+ /** Determine order from pre-existing mesh.
+ *
+ * @param mesh Finite-element mesh.
+ */
+ void initialize(const ALE::Obj<mesh_type>& mesh);
+
+ /** Set range for normal cells.
+ *
+ * @param min Minimum cell label.
+ * @param max Maximum cell label.
+ */
+ void cellsNormal(const point_type min, const point_type max);
+
+ /** Get range for normal cells.
+ *
+ * @returns Interval with range of normal cells.
+ */
+ const ALE::Interval<point_type>& cellsNormal(void) const;
+
+ /** Set range for normal vertices.
+ *
+ * @param min Minimum vertex label.
+ * @param max Maximum vertex label.
+ */
+ void verticesNormal(const point_type min, const point_type max);
+
+ /** Get range for normal vertices.
+ *
+ * @returns Interval with range of normal vertices.
+ */
+ const ALE::Interval<point_type>& verticesNormal(void) const;
+
+ /** Set range for censored cells.
+ *
+ * @param min Minimum cell label.
+ * @param max Maximum cell label.
+ */
+ void cellsCensored(const point_type min, const point_type max);
+
+ /** Get range for censored cells.
+ *
+ * @returns Interval with range of censored cells.
+ */
+ const ALE::Interval<point_type>& cellsCensored(void) const;
+
+ /** Set range for censored vertices.
+ *
+ * @param min Minimum vertex label.
+ * @param max Maximum vertex label.
+ */
+ void verticesCensored(const point_type min, const point_type max);
+
+ /** Get range for censored vertices.
+ *
+ * @returns Interval with range of censored vertices.
+ */
+ const ALE::Interval<point_type>& verticesCensored(void) const;
+
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ ALE::Interval<point_type> _cellsNormal;
+ ALE::Interval<point_type> _verticesNormal;
+ ALE::Interval<point_type> _verticesCensored;
+ ALE::Interval<point_type> _cellsCensored;
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ MeshOrder(const MeshOrder&); ///< Not implemented
+ const MeshOrder& operator=(const MeshOrder&); ///< Not implemented
+
+}; // MeshOrder
+
+#endif // pylith_topology_meshorder_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/MeshOrder.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshOrder.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshOrder.icc 2010-09-24 21:33:18 UTC (rev 17213)
@@ -0,0 +1,26 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+#if !defined(pylith_topology_meshorder_hh)
+#error "MeshOrder.icc must be included only from MeshOrder.hh"
+#else
+
+#endif
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc 2010-09-24 21:33:18 UTC (rev 17213)
@@ -0,0 +1,387 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "MeshRefiner.hh" // implementation of class methods
+
+#include "CellRefinerTri3.hh" // USES CellRefinerTri3
+#include "MeshOrder.hh" // USES MeshOrder
+
+#include <cassert> // USES assert()
+
+// ----------------------------------------------------------------------
+// Constructor
+ALE::MeshRefiner::MeshRefiner(void) :
+ _orderOldMesh(new MeshOrder),
+ _orderNewMesh(new MeshOrder)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+ALE::MeshRefiner::~MeshRefiner(void)
+{ // destructor
+ delete _orderOldMesh; _orderOldMesh = PETSC_NULL;
+ delete _orderNewMesh; _orderNewMesh = PETSC_NULL;
+} // destructor
+
+// ----------------------------------------------------------------------
+// Refine mesh.
+void
+ALE::MeshRefiner::refine(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ CellRefinerTri3& refiner)
+{ // refine
+ assert(!mesh.isNull());
+ if (mesh->hasLabel("censored depth")) {
+ _refineCensored(newMesh, mesh, refiner);
+ } else {
+ _refine(newMesh, mesh, refiner);
+ } // if/else
+} // refine
+
+// ----------------------------------------------------------------------
+// Refine mesh.
+void
+ALE::MeshRefiner::_refine(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ CellRefinerTri3& refiner)
+{ // _refine
+ typedef Interval<point_type>::const_iterator interval_type;
+
+ assert(_orderOldMesh);
+ assert(_orderNewMesh);
+
+ // Calculate order in old mesh.
+ _orderOldMesh->initialize(mesh);
+
+ assert(!mesh.isNull());
+ assert(!newMesh.isNull());
+
+ // Get original mesh stuff.
+ const Obj<mesh_type::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const mesh_type::label_sequence::iterator cellsEnd = cells->end();
+
+ const Obj<mesh_type::label_sequence>& vertices = mesh->depthStratum(0);
+ assert(!vertices.isNull());
+
+ const Obj<mesh_type::sieve_type>& sieve = mesh->getSieve();
+ assert(!sieve.isNull());
+ ALE::ISieveVisitor::PointRetriever<mesh_type::sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
+
+ // Count number of cells.
+ const int oldNumCells = cells->size();
+ int newNumCells = 0;
+ for (mesh_type::label_sequence::iterator c_iter = cells->begin(); c_iter != cellsEnd; ++c_iter) {
+ newNumCells += refiner.numNewCells(*c_iter);
+ } // for
+
+ // Count number of vertices.
+ const int oldNumVertices = vertices->size();
+ int counterBegin = newNumCells + oldNumVertices;
+ point_type curNewVertex = counterBegin;
+ for(mesh_type::label_sequence::iterator c_iter = cells->begin(); c_iter != cellsEnd; ++c_iter) {
+ cV.clear();
+ sieve->cone(*c_iter, cV);
+ refiner.splitCell(*c_iter, cV.getPoints(), cV.getSize(), &curNewVertex);
+ } // for
+ const int newNumVertices = oldNumVertices + curNewVertex - counterBegin;
+
+ _orderNewMesh->cellsNormal(0, newNumCells);
+ _orderNewMesh->verticesNormal(newNumCells, newNumCells+newNumVertices);
+ _orderNewMesh->verticesCensored(newNumCells+newNumVertices, newNumCells+newNumVertices);
+ _orderNewMesh->cellsCensored(newNumCells+newNumVertices, newNumCells+newNumVertices);
+
+ // Allocate chart for new sieve.
+ const Obj<mesh_type::sieve_type>& newSieve = newMesh->getSieve();
+ assert(!newSieve.isNull());
+ newSieve->setChart(mesh_type::sieve_type::chart_type(0, _orderNewMesh->cellsCensored().max()));
+
+ // Create new sieve with correct sizes for refined cells
+ point_type curNewCell = _orderNewMesh->cellsNormal().min();
+ const interval_type::const_iterator oldCellsEnd = _orderOldMesh->cellsNormal().end();
+ for (interval_type::const_iterator c_iter=_orderOldMesh->cellsNormal().begin(); c_iter != oldCellsEnd; ++c_iter) {
+ // Set new cone and support sizes
+ cV.clear();
+ sieve->cone(*c_iter, cV);
+ const point_type* cone = cV.getPoints();
+ const int coneSize = cV.getSize();
+
+ const point_type* newCells;
+ int numNewCells = 0;
+ refiner.getNewCells(&newCells, &numNewCells, *c_iter, cone, coneSize, *_orderOldMesh, *_orderNewMesh);
+
+ for(int iCell=0; iCell < numNewCells; ++iCell, ++curNewCell) {
+ newSieve->setConeSize(curNewCell, coneSize);
+ for(int iVertex=0; iVertex < coneSize; ++iVertex) {
+ newSieve->addSupportSize(newCells[iCell*coneSize+iVertex], 1);
+ } // for
+ } // for
+
+ } // for
+ newSieve->allocate();
+
+ // Create refined cells in new sieve.
+ curNewCell = _orderNewMesh->cellsNormal().min();
+ for (interval_type::const_iterator c_iter=_orderOldMesh->cellsNormal().begin(); c_iter != oldCellsEnd; ++c_iter) {
+ cV.clear();
+ sieve->cone(*c_iter, cV);
+ const point_type *cone = cV.getPoints();
+ const int coneSize = cV.getSize();
+
+ const point_type* newCells;
+ int numNewCells = 0;
+ refiner.getNewCells(&newCells, &numNewCells, *c_iter, cone, coneSize, *_orderOldMesh, *_orderNewMesh);
+
+ for(int iCell=0; iCell < numNewCells; ++iCell, ++curNewCell) {
+ newSieve->setCone(&newCells[iCell*coneSize], curNewCell);
+ } // for
+ } // for
+ newSieve->symmetrize();
+
+ // Set coordinates in refined mesh.
+ const Obj<mesh_type::real_section_type>& coordinates = mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ const Obj<mesh_type::real_section_type>& newCoordinates = newMesh->getRealSection("coordinates");
+ assert(!newCoordinates.isNull());
+
+ const mesh_type::label_sequence::const_iterator verticesEnd = vertices->end();
+ assert(vertices->size() > 0);
+ const int spaceDim = coordinates->getFiberDimension(*vertices->begin());
+ assert(spaceDim > 0);
+ newCoordinates->setChart(mesh_type::sieve_type::chart_type(_orderNewMesh->verticesNormal().min(), _orderNewMesh->verticesCensored().max()));
+
+ interval_type::const_iterator newVerticesEnd = _orderNewMesh->verticesCensored().end();
+ for (interval_type::const_iterator v_iter=_orderNewMesh->verticesNormal().begin(); v_iter != newVerticesEnd; ++v_iter) {
+ newCoordinates->setFiberDimension(*v_iter, spaceDim);
+ } // for
+ newCoordinates->allocatePoint();
+
+ const interval_type::const_iterator oldVerticesEnd = _orderOldMesh->verticesCensored().end();
+ for (interval_type::const_iterator vOld_iter=_orderOldMesh->verticesNormal().begin(), vNew_iter=_orderNewMesh->verticesNormal().begin(); vOld_iter != oldVerticesEnd; ++vOld_iter, ++vNew_iter) {
+ std::cout << "Copy coordinates from old vertex " << *vOld_iter << " to new vertex " << *vNew_iter << std::endl;
+ newCoordinates->updatePoint(*vNew_iter, coordinates->restrictPoint(*vOld_iter));
+ } // for
+
+ refiner.setCoordsNewVertices(newCoordinates, coordinates);
+
+ _stratify(newMesh);
+ _calcNewOverlap(newMesh, mesh);
+} // _refine
+
+// ----------------------------------------------------------------------
+// Refine mesh with a censored depth.
+ void
+ ALE::MeshRefiner::_refineCensored(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ CellRefinerTri3& refiner)
+{ // _refineCensored
+} // _refineCensored
+
+// ----------------------------------------------------------------------
+// Stratify mesh.
+void
+ALE::MeshRefiner::_stratify(const Obj<mesh_type>& mesh)
+{ // _stratify
+ typedef Interval<point_type>::const_iterator interval_type;
+
+ assert(_orderNewMesh);
+ assert(_orderOldMesh);
+
+ // Fast stratification
+ const Obj<mesh_type::label_type>& height = mesh->createLabel("height");
+ const Obj<mesh_type::label_type>& depth = mesh->createLabel("depth");
+
+ // Set height/depth of cells
+ interval_type::const_iterator cellsEnd = _orderNewMesh->cellsNormal().end();
+ for (interval_type::const_iterator c_iter = _orderNewMesh->cellsNormal().begin(); c_iter != cellsEnd; ++c_iter) {
+ height->setCone(0, *c_iter);
+ depth->setCone(1, *c_iter);
+ } // for
+ cellsEnd = _orderNewMesh->cellsCensored().end();
+ for (interval_type::const_iterator c_iter = _orderNewMesh->cellsCensored().begin(); c_iter != cellsEnd; ++c_iter) {
+ height->setCone(0, *c_iter);
+ depth->setCone(1, *c_iter);
+ } // for
+
+ // Set height/depth of vertices
+ interval_type::const_iterator verticesEnd = _orderNewMesh->verticesNormal().end();
+ for (interval_type::const_iterator v_iter = _orderNewMesh->verticesNormal().begin(); v_iter != verticesEnd; ++v_iter) {
+ height->setCone(1, *v_iter);
+ depth->setCone(0, *v_iter);
+ } // for
+ verticesEnd = _orderNewMesh->verticesCensored().end();
+ for (interval_type::const_iterator v_iter = _orderNewMesh->verticesCensored().begin(); v_iter != verticesEnd; ++v_iter) {
+ height->setCone(1, *v_iter);
+ depth->setCone(0, *v_iter);
+ } // for
+
+ mesh->setHeight(1);
+ mesh->setDepth(1);
+} // _stratify
+
+// ----------------------------------------------------------------------
+// Calculate new overlap.
+void
+ALE::MeshRefiner::_calcNewOverlap(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh)
+{ // _calcNewOverlap
+ // Exchange new boundary vertices
+ // We can convert endpoints, and then just match to new vertex on this side
+ // 1) Create the overlap of edges which are vertex pairs (do not need for interpolated meshes)
+ // 2) Create a section of overlap edge --> new vertex (this will generalize to other split points in interpolated meshes)
+ // 3) Copy across new overlap
+ // 4) Fuse matches new vertex pairs and inserts them into the old overlap
+
+ // Create the parallel overlap
+#if 0
+ int *oldNumCellsP = new int[mesh->commSize()];
+ int *newNumCellsP = new int[newMesh->commSize()];
+ int ierr;
+
+ ierr = MPI_Allgather((void *) &oldNumCells, 1, MPI_INT, oldNumCellsP, 1, MPI_INT, mesh->comm());CHKERRXX(ierr);
+ ierr = MPI_Allgather((void *) &newNumCells, 1, MPI_INT, newNumCellsP, 1, MPI_INT, newMesh->comm());CHKERRXX(ierr);
+ Obj<mesh_type::send_overlap_type> newSendOverlap = newMesh->getSendOverlap();
+ Obj<mesh_type::recv_overlap_type> newRecvOverlap = newMesh->getRecvOverlap();
+ const Obj<mesh_type::send_overlap_type>& sendOverlap = mesh->getSendOverlap();
+ const Obj<mesh_type::recv_overlap_type>& recvOverlap = mesh->getRecvOverlap();
+ Obj<mesh_type::send_overlap_type::traits::capSequence> sendPoints = sendOverlap->cap();
+ const mesh_type::send_overlap_type::source_type localOffset = newNumCellsP[newMesh->commRank()] - oldNumCellsP[mesh->commRank()];
+
+ for(mesh_type::send_overlap_type::traits::capSequence::iterator p_iter = sendPoints->begin(); p_iter != sendPoints->end(); ++p_iter) {
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& ranks = sendOverlap->support(*p_iter);
+ const mesh_type::send_overlap_type::source_type& localPoint = *p_iter;
+
+ for(mesh_type::send_overlap_type::traits::supportSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
+ const int rank = *r_iter;
+ const mesh_type::send_overlap_type::source_type& remotePoint = r_iter.color();
+ const mesh_type::send_overlap_type::source_type remoteOffset = newNumCellsP[rank] - oldNumCellsP[rank];
+
+ newSendOverlap->addArrow(localPoint+localOffset, rank, remotePoint+remoteOffset);
+ }
+ }
+ Obj<mesh_type::recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();
+
+ for(mesh_type::recv_overlap_type::traits::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
+ const Obj<mesh_type::recv_overlap_type::traits::coneSequence>& ranks = recvOverlap->cone(*p_iter);
+ const mesh_type::recv_overlap_type::target_type& localPoint = *p_iter;
+
+ for(mesh_type::recv_overlap_type::traits::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
+ const int rank = *r_iter;
+ const mesh_type::recv_overlap_type::target_type& remotePoint = r_iter.color();
+ const mesh_type::recv_overlap_type::target_type remoteOffset = newNumCellsP[rank] - oldNumCellsP[rank];
+
+ newRecvOverlap->addArrow(rank, localPoint+localOffset, remotePoint+remoteOffset);
+ }
+ }
+ newMesh->setCalculatedOverlap(true);
+ delete [] oldNumCellsP;
+ delete [] newNumCellsP;
+ // Check edges in edge2vertex for both endpoints sent to same process
+ // Put it in section with point being the lowest numbered vertex and value (other endpoint, new vertex)
+ Obj<ALE::Section<point_type, edge_type> > newVerticesSection = new ALE::Section<point_type, edge_type>(mesh->comm());
+ std::map<edge_type, std::vector<int> > bdedge2rank;
+
+ for(std::map<edge_type, point_type>::const_iterator e_iter = edge2vertex.begin(); e_iter != edge2vertex.end(); ++e_iter) {
+ const point_type left = e_iter->first.first;
+ const point_type right = e_iter->first.second;
+
+ if (sendOverlap->capContains(left) && sendOverlap->capContains(right)) {
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& leftRanksSeq = sendOverlap->support(left);
+ std::list<int> leftRanks(leftRanksSeq->begin(), leftRanksSeq->end());
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& rightRanks = sendOverlap->support(right);
+ std::list<int> ranks;
+ std::set_intersection(leftRanks.begin(), leftRanks.end(), rightRanks->begin(), rightRanks->end(),
+ std::insert_iterator<std::list<int> >(ranks, ranks.begin()));
+
+ if(ranks.size()) {
+ newVerticesSection->addFiberDimension(std::min(e_iter->first.first, e_iter->first.second)+localOffset, 1);
+ for(std::list<int>::const_iterator r_iter = ranks.begin(); r_iter != ranks.end(); ++r_iter) {
+ bdedge2rank[e_iter->first].push_back(*r_iter);
+ }
+ }
+ }
+ }
+ newVerticesSection->allocatePoint();
+ const ALE::Section<point_type, edge_type>::chart_type& chart = newVerticesSection->getChart();
+
+ for(ALE::Section<point_type, edge_type>::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
+ typedef ALE::Section<point_type, edge_type>::value_type value_type;
+ const point_type p = *c_iter;
+ const int dim = newVerticesSection->getFiberDimension(p);
+ int v = 0;
+ value_type *values = new value_type[dim];
+
+ for(std::map<edge_type, std::vector<int> >::const_iterator e_iter = bdedge2rank.begin(); e_iter != bdedge2rank.end() && v < dim; ++e_iter) {
+ if (std::min(e_iter->first.first, e_iter->first.second)+localOffset == p) {
+ values[v++] = edge_type(std::max(e_iter->first.first, e_iter->first.second)+localOffset, edge2vertex[e_iter->first]);
+ }
+ }
+ newVerticesSection->updatePoint(p, values);
+ delete [] values;
+ }
+ // Copy across overlap
+ typedef ALE::Pair<int, point_type> overlap_point_type;
+ Obj<ALE::Section<overlap_point_type, edge_type> > overlapVertices = new ALE::Section<overlap_point_type, edge_type>(mesh->comm());
+
+ ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newVerticesSection, overlapVertices);
+ // Merge by translating edge to local points, finding edge in edge2vertex, and adding (local new vetex, remote new vertex) to overlap
+ for(std::map<edge_type, std::vector<int> >::const_iterator e_iter = bdedge2rank.begin(); e_iter != bdedge2rank.end(); ++e_iter) {
+ const point_type localPoint = edge2vertex[e_iter->first];
+
+ for(std::vector<int>::const_iterator r_iter = e_iter->second.begin(); r_iter != e_iter->second.end(); ++r_iter) {
+ point_type remoteLeft = -1, remoteRight = -1;
+ const int rank = *r_iter;
+
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& leftRanks = newSendOverlap->support(e_iter->first.first+localOffset);
+ for(mesh_type::send_overlap_type::traits::supportSequence::iterator lr_iter = leftRanks->begin(); lr_iter != leftRanks->end(); ++lr_iter) {
+ if (rank == *lr_iter) {
+ remoteLeft = lr_iter.color();
+ break;
+ }
+ }
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& rightRanks = newSendOverlap->support(e_iter->first.second+localOffset);
+ for(mesh_type::send_overlap_type::traits::supportSequence::iterator rr_iter = rightRanks->begin(); rr_iter != rightRanks->end(); ++rr_iter) {
+ if (rank == *rr_iter) {
+ remoteRight = rr_iter.color();
+ break;
+ }
+ }
+ const point_type remoteMin = std::min(remoteLeft, remoteRight);
+ const point_type remoteMax = std::max(remoteLeft, remoteRight);
+ const int remoteSize = overlapVertices->getFiberDimension(overlap_point_type(rank, remoteMin));
+ const edge_type *remoteVals = overlapVertices->restrictPoint(overlap_point_type(rank, remoteMin));
+ point_type remotePoint = -1;
+
+ for(int d = 0; d < remoteSize; ++d) {
+ if (remoteVals[d].first == remoteMax) {
+ remotePoint = remoteVals[d].second;
+ break;
+ }
+ }
+ newSendOverlap->addArrow(localPoint, rank, remotePoint);
+ newRecvOverlap->addArrow(rank, localPoint, remotePoint);
+ }
+ }
+#endif
+} // _calcNewOverlap
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh 2010-09-24 21:33:18 UTC (rev 17213)
@@ -0,0 +1,114 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/topology/MeshRefiner.hh
+ *
+ * @brief Object for tri3 refinement of cells.
+ */
+
+#if !defined(pylith_topology_meshrefiner_hh)
+#define pylith_topology_meshrefiner_hh
+
+// Include directives ---------------------------------------------------
+#include "topologyfwd.hh" // forward declarations
+
+// RefineTri3 --------------------------------------------------------
+/// Object for tri3 refinement of cells.
+class ALE::MeshRefiner
+{ // MeshRefiner
+ typedef IMesh<> mesh_type;
+ typedef mesh_type::point_type point_type;
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Constructor
+ *
+ * @param mesh Finite-element mesh.
+ */
+ MeshRefiner(void);
+
+ /// Destructor
+ ~MeshRefiner(void);
+
+ /** Refine mesh.
+ *
+ * @param newMesh New mesh.
+ * @param mesh Current mesh.
+ * @param refiner Cell refiner.
+ */
+ void refine(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ CellRefinerTri3& refiner);
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /** Refine mesh.
+ *
+ * @param newMesh New mesh.
+ * @param mesh Current mesh.
+ * @param refiner Cell refiner.
+ */
+ void _refine(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ CellRefinerTri3& refiner);
+
+ /** Refine mesh with a censored depth.
+ *
+ * @param newMesh New mesh.
+ * @param mesh Current mesh.
+ * @param refiner Cell refiner.
+ */
+ void _refineCensored(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ CellRefinerTri3& refiner);
+
+ /** Stratify mesh.
+ *
+ * @param mesh Mesh to stratify.
+ */
+ void _stratify(const Obj<mesh_type>& mesh);
+
+ /** Calculate new overlap.
+ *
+ * @param newMesh New (refined) mesh.
+ * @param mesh Current (unrefined) mesh with overlap.
+ */
+ void _calcNewOverlap(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ MeshOrder* _orderOldMesh; ///< Order of entities in old mesh.
+ MeshOrder* _orderNewMesh; ///< Order of entities in new mesh.
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ MeshRefiner(const MeshRefiner&); ///< Not implemented
+ const MeshRefiner& operator=(const MeshRefiner&); ///< Not implemented
+
+}; // MeshRefiner
+
+#endif // pylith_topology_meshrefiner_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc 2010-09-24 03:35:36 UTC (rev 17212)
+++ short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc 2010-09-24 21:33:18 UTC (rev 17213)
@@ -22,6 +22,9 @@
#include "Mesh.hh" // USES Mesh
+#include "CellRefinerTri3.hh" // USES CellRefinerTri3
+#include "MeshRefiner.hh" // USES MeshRefiner
+
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
#include <cassert> // USES assert()
@@ -30,20 +33,6 @@
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
// ----------------------------------------------------------------------
-template<typename Point>
-class Edge : public std::pair<Point, Point> {
-public:
- Edge() : std::pair<Point, Point>() {};
- Edge(const Point l) : std::pair<Point, Point>(l, l) {};
- Edge(const Point l, const Point r) : std::pair<Point, Point>(l, r) {};
- ~Edge() {};
- friend std::ostream& operator<<(std::ostream& stream, const Edge& edge) {
- stream << "(" << edge.first << ", " << edge.second << ")";
- return stream;
- };
-};
-
-// ----------------------------------------------------------------------
// Constructor
pylith::topology::RefineUniform::RefineUniform(void)
{ // constructor
@@ -65,7 +54,6 @@
assert(0 != newMesh);
typedef SieveMesh::point_type point_type;
- typedef Edge<point_type> edge_type;
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
assert(!sieveMesh.isNull());
@@ -75,11 +63,10 @@
new SieveMesh::sieve_type(mesh.comm(), mesh.debug());
newSieveMesh->setSieve(newSieve);
- ALE::MeshBuilder<SieveMesh>::CellRefiner<SieveMesh,edge_type> refiner(*sieveMesh);
+ ALE::CellRefinerTri3 cellSplitter(*sieveMesh);
+ ALE::MeshRefiner refiner;
+ refiner.refine(newSieveMesh, sieveMesh, cellSplitter);
- ALE::MeshBuilder<SieveMesh>::refineGeneral< SieveMesh,
- ALE::MeshBuilder<SieveMesh>::CellRefiner<SieveMesh,edge_type> >(*sieveMesh, *newSieveMesh, refiner);
-
// Set material ids
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->heightStratum(0);
@@ -106,7 +93,7 @@
cNew_iter = newCellsBegin;
c_iter != cellsEnd;
++c_iter) {
- const int numNewCellsPerCell = refiner.numNewCells(*c_iter);
+ const int numNewCellsPerCell = cellSplitter.numNewCells(*c_iter);
const int material = sieveMesh->getValue(materials, *c_iter);
for(int i=0; i < numNewCellsPerCell; ++i, ++cNew_iter)
@@ -119,6 +106,7 @@
const ALE::Obj<std::set<std::string> >& sectionNames =
sieveMesh->getIntSections();
+#if 0
ALE::MeshBuilder<SieveMesh>::CellRefiner<SieveMesh,edge_type>::edge_map_type& edge2vertex =
refiner.getEdgeToVertex();
@@ -171,6 +159,7 @@
newGroup->updatePoint(e_iter->second, group->restrictPoint(vertexA));
} // for
} // for
+#endif
} // refine
Modified: short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh 2010-09-24 03:35:36 UTC (rev 17212)
+++ short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh 2010-09-24 21:33:18 UTC (rev 17213)
@@ -33,6 +33,10 @@
template<typename point_type,
typename value_type,
typename allocator> class IUniformSectionDS;
+
+ class MeshRefiner;
+ class CellRefinerTri3;
+ class MeshOrder;
} // ALE
namespace pylith {
@@ -53,7 +57,6 @@
class Distributor;
- class MeshRefiner;
class RefineUniform;
class ReverseCuthillMcKee;
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc 2010-09-24 03:35:36 UTC (rev 17212)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc 2010-09-24 21:33:18 UTC (rev 17213)
@@ -130,6 +130,8 @@
faultB.adjustTopology(mesh, &firstFaultVertex,
&firstLagrangeVertex, &firstFaultCell);
} // if
+
+ mesh->view("UNREFINED MESH");
} // _setupMesh
// ----------------------------------------------------------------------
@@ -194,7 +196,7 @@
// Normal cells
ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> pV(sieve->getMaxConeSize());
- const int offset = numCells;
+ const int offset = numCells - data.numCellsCohesive;
SieveMesh::label_sequence::iterator c_iter = cells->begin();
for(int iCell=0, i=0; iCell < data.numCells; ++iCell, ++c_iter) {
pV.clear();
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.hh 2010-09-24 03:35:36 UTC (rev 17212)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.hh 2010-09-24 21:33:18 UTC (rev 17213)
@@ -52,9 +52,11 @@
CPPUNIT_TEST( testConstructor );
CPPUNIT_TEST( testRefineTri3Level2 );
+#if 0
CPPUNIT_TEST( testRefineTri3Level2Fault1 );
CPPUNIT_TEST( testRefineTet4Level2 );
CPPUNIT_TEST( testRefineTet4Level2Fault1 );
+#endif
CPPUNIT_TEST_SUITE_END();
More information about the CIG-COMMITS
mailing list