[cig-commits] r17223 - in short/3D/PyLith/trunk/libsrc: . topology
brad at geodynamics.org
brad at geodynamics.org
Mon Sep 27 14:06:53 PDT 2010
Author: brad
Date: 2010-09-27 14:06:53 -0700 (Mon, 27 Sep 2010)
New Revision: 17223
Added:
short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.cc
short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.hh
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc
short/3D/PyLith/trunk/libsrc/topology/Makefile.am
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/topologyfwd.hh
Log:
Made MeshRefiner templated over cell refiner. Cleaner implementation for creating new labels. Added cell refiner for tet4 meshes.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2010-09-27 19:32:24 UTC (rev 17222)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2010-09-27 21:06:53 UTC (rev 17223)
@@ -141,8 +141,8 @@
# TEMPORARY
libpylith_la_SOURCES += \
- topology/MeshRefiner.cc \
topology/CellRefinerTri3.cc \
+ topology/CellRefinerTet4.cc \
topology/MeshOrder.cc
Added: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.cc 2010-09-27 21:06:53 UTC (rev 17223)
@@ -0,0 +1,466 @@
+// -*- 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 "CellRefinerTet4.hh" // implementation of class methods
+
+#include "MeshOrder.hh" // USES MeshOrder
+
+#include <cassert> // USES assert()
+
+#include <iostream> // TEMPORARY
+// ----------------------------------------------------------------------
+// Constructor
+ALE::CellRefinerTet4::CellRefinerTet4(const mesh_type& mesh) :
+ _mesh(mesh)
+{ // constructor
+ assert(3 == mesh.getDimension());
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+ALE::CellRefinerTet4::~CellRefinerTet4(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Get number of refined cells for each original cell.
+int
+ALE::CellRefinerTet4::numNewCells(const point_type cell)
+{ // numNewCells
+ switch (_cellType(cell)) {
+ case TETRAHEDRON:
+ return 8;
+ case TRIANGLE_COHESIVE_LAGRANGE:
+ return 4;
+ default:
+ assert(0);
+ throw ALE::Exception("Unknown cell type.");
+ } // switch
+} // numNewCells
+
+// ----------------------------------------------------------------------
+// Split cell into smaller cells of same type.
+void
+ALE::CellRefinerTet4::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 TETRAHEDRON:
+ _edges_TETRAHEDRON(&edges, &numEdges, cone, coneSize);
+ break;
+ case TRIANGLE_COHESIVE_LAGRANGE:
+ _edges_TRIANGLE_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::CellRefinerTet4::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 TETRAHEDRON: {
+ const int coneVertexOffset = orderNewMesh.verticesNormal().min() - orderOldMesh.verticesNormal().min();
+ _newCells_TETRAHEDRON(cells, numCells, cone, coneSize, coneVertexOffset);
+ break;
+ } // TETRAHEDRON
+ case TRIANGLE_COHESIVE_LAGRANGE: {
+ const int coneVertexOffsetNormal = orderNewMesh.verticesNormal().min() - orderOldMesh.verticesNormal().min();
+ const int coneVertexOffsetCensored = orderNewMesh.verticesCensored().min() - orderOldMesh.verticesCensored().min();
+ _newCells_TRIANGLE_COHESIVE_LAGRANGE(cells, numCells, cone, coneSize, coneVertexOffsetNormal, coneVertexOffsetCensored);
+ break;
+ } // TRIANGLE_COHESIVE_LAGRANGE
+ default:
+ throw ALE::Exception("Unknown cell type.");
+ } // switch
+} // getNewCells
+
+// ----------------------------------------------------------------------
+// Set coordinates of new vertices.
+void
+ALE::CellRefinerTet4::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
+
+// ----------------------------------------------------------------------
+// Add space for new vertices in group.
+void
+ALE::CellRefinerTet4::groupAddNewVertices(const ALE::Obj<mesh_type::int_section_type>& newGroup,
+ const ALE::Obj<mesh_type::int_section_type>& oldGroup)
+{ // groupAddNewVertices
+ assert(!newGroup.isNull());
+ assert(!oldGroup.isNull());
+
+ 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;
+
+ if (oldGroup->getFiberDimension(edgeVertexA) && oldGroup->getFiberDimension(edgeVertexB)) {
+ if (oldGroup->restrictPoint(edgeVertexA)[0] == oldGroup->restrictPoint(edgeVertexB)[0]) {
+ newGroup->setFiberDimension(newVertex, 1);
+ } // if
+ } // if
+ } // for
+} // groupAddNewVertices
+
+// ----------------------------------------------------------------------
+// Set new vertices in group.
+void
+ALE::CellRefinerTet4::groupSetNewVertices(const ALE::Obj<mesh_type::int_section_type>& newGroup,
+ const ALE::Obj<mesh_type::int_section_type>& oldGroup)
+{ // groupSetNewVertices
+ assert(!newGroup.isNull());
+ assert(!oldGroup.isNull());
+
+ 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;
+
+ if (oldGroup->getFiberDimension(edgeVertexA) && oldGroup->getFiberDimension(edgeVertexB)) {
+ if (oldGroup->restrictPoint(edgeVertexA)[0] == oldGroup->restrictPoint(edgeVertexB)[0]) {
+ newGroup->updatePoint(newVertex, oldGroup->restrictPoint(edgeVertexA));
+ std::cout << "Adding new vertex: " << newVertex << " based on old vertices " << edgeVertexA << " and " << edgeVertexB << std::endl;
+ } // if
+ } // if
+ } // for
+} // groupSetNewVertices
+
+// ----------------------------------------------------------------------
+// Add new vertices to label.
+void
+ALE::CellRefinerTet4::labelAddNewVertices(const ALE::Obj<mesh_type>& newMesh,
+ const ALE::Obj<mesh_type>& oldMesh,
+ const char* labelName)
+{ // labelAddNewVertices
+ assert(!newMesh.isNull());
+ assert(!oldMesh.isNull());
+
+ const Obj<mesh_type::label_sequence>& oldLabelVertices = oldMesh->getLabelStratum(labelName, 0);
+ assert(!oldLabelVertices.isNull());
+
+ const Obj<mesh_type::label_type>& oldLabel = oldMesh->getLabel(labelName);
+ assert(!oldLabel.isNull());
+ const Obj<mesh_type::label_type>& newLabel = newMesh->getLabel(labelName);
+ assert(!newLabel.isNull());
+
+ const int defaultValue = -999;
+
+ 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;
+
+ const int valueA = oldMesh->getValue(oldLabel, edgeVertexA, defaultValue);
+ const int valueB = oldMesh->getValue(oldLabel, edgeVertexB, defaultValue);
+
+ if (valueA != defaultValue && valueA == valueB) {
+ newMesh->setValue(newLabel, newVertex, valueA);
+ } // if
+ } // for
+} // labelAddNewVertices
+
+// ----------------------------------------------------------------------
+// Get cell type.
+ALE::CellRefinerTet4::CellEnum
+ALE::CellRefinerTet4::_cellType(const point_type cell)
+{ // _cellType
+ assert(!_mesh.getSieve().isNull());
+
+ switch (_mesh.getSieve()->getConeSize(cell)) {
+ case 4:
+ return TETRAHEDRON;
+ case 9:
+ return TRIANGLE_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::CellRefinerTet4::_edges_TETRAHEDRON(const EdgeType** edges,
+ int* numEdges,
+ const point_type cone[],
+ const int coneSize)
+{ // _edges_TETRAHEDRON
+ static EdgeType splitEdges[6];
+
+ assert(coneSize == 4);
+ splitEdges[0] = EdgeType(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
+ splitEdges[1] = EdgeType(std::min(cone[1], cone[2]), std::max(cone[1], cone[2]));
+ splitEdges[2] = EdgeType(std::min(cone[2], cone[0]), std::max(cone[2], cone[0]));
+ splitEdges[3] = EdgeType(std::min(cone[0], cone[3]), std::max(cone[0], cone[3]));
+ splitEdges[4] = EdgeType(std::min(cone[1], cone[3]), std::max(cone[1], cone[3]));
+ splitEdges[5] = EdgeType(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]));
+ *numEdges = 6;
+ *edges = splitEdges;
+} // _edges_TETRAHEDRON
+
+// ----------------------------------------------------------------------
+// Get edges of line cohesive cell with Lagrange multipler vertices.
+void
+ALE::CellRefinerTet4::_edges_TRIANGLE_COHESIVE_LAGRANGE(const EdgeType** edges,
+ int* numEdges,
+ const point_type cone[],
+ const int coneSize)
+{ // _edges_TRIANGLE_COHESIVE_LAGRANGE
+ static EdgeType splitEdges[9];
+
+ assert(coneSize == 9);
+ splitEdges[0] = EdgeType(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
+ splitEdges[1] = EdgeType(std::min(cone[1], cone[2]), std::max(cone[1], cone[2]));
+ splitEdges[2] = EdgeType(std::min(cone[2], cone[0]), std::max(cone[2], cone[0]));
+ splitEdges[3] = EdgeType(std::min(cone[3], cone[4]), std::max(cone[3], cone[4]));
+ splitEdges[4] = EdgeType(std::min(cone[4], cone[5]), std::max(cone[4], cone[5]));
+ splitEdges[5] = EdgeType(std::min(cone[5], cone[3]), std::max(cone[5], cone[3]));
+ splitEdges[6] = EdgeType(std::min(cone[6], cone[7]), std::max(cone[6], cone[7]));
+ splitEdges[7] = EdgeType(std::min(cone[7], cone[8]), std::max(cone[7], cone[8]));
+ splitEdges[8] = EdgeType(std::min(cone[8], cone[6]), std::max(cone[8], cone[6]));
+ *numEdges = 9;
+ *edges = splitEdges;
+} // _edges_TRIANGLE_COHESIVE_LAGRANGE
+
+// ----------------------------------------------------------------------
+// Get new cells from refinement of a triangular cell.
+void
+ALE::CellRefinerTet4::_newCells_TETRAHEDRON(const point_type** cells,
+ int *numCells,
+ const point_type cone[],
+ const int coneSize,
+ const int coneVertexOffset)
+{ // _newCells_TETRAHEDRON
+ const int coneSizeTet4 = 4;
+ const int numEdgesTet4 = 6;
+ const int numNewCells = 8;
+ const int numNewVertices = 6;
+
+ int numEdges = 0;
+ const EdgeType *edges;
+ _edges_TETRAHEDRON(&edges, &numEdges, cone, coneSize);
+ assert(numEdgesTet4 == numEdges);
+
+ static point_type newCells[numNewCells*coneSizeTet4];
+ point_type newVertices[numNewVertices];
+ for(int iEdge=0, iNewVertex=0; iEdge < numEdgesTet4; ++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
+ newCells[0*4+0] = cone[0]+coneVertexOffset;
+ newCells[0*4+1] = newVertices[3];
+ newCells[0*4+2] = newVertices[0];
+ newCells[0*4+3] = newVertices[2];
+
+ // new cell 1
+ newCells[1*4+0] = newVertices[0];
+ newCells[1*4+1] = newVertices[1];
+ newCells[1*4+2] = newVertices[2];
+ newCells[1*4+3] = newVertices[3];
+
+ // new cell 2
+ newCells[2*4+0] = newVertices[0];
+ newCells[2*4+1] = newVertices[3];
+ newCells[2*4+2] = newVertices[4];
+ newCells[2*4+3] = newVertices[1];
+
+ // new cell 3
+ newCells[3*4+0] = cone[1]+coneVertexOffset;
+ newCells[3*4+1] = newVertices[4];
+ newCells[3*4+2] = newVertices[1];
+ newCells[3*4+3] = newVertices[0];
+
+ // new cell 4
+ newCells[4*4+0] = newVertices[2];
+ newCells[4*4+1] = newVertices[5];
+ newCells[4*4+2] = newVertices[3];
+ newCells[4*4+3] = newVertices[1];
+
+ // new cell 5
+ newCells[5*4+0] = cone[2]+coneVertexOffset;
+ newCells[5*4+1] = newVertices[5];
+ newCells[5*4+2] = newVertices[2];
+ newCells[5*4+3] = newVertices[1];
+
+ // new cell 6
+ newCells[6*4+0] = newVertices[1];
+ newCells[6*4+1] = newVertices[4];
+ newCells[6*4+2] = newVertices[5];
+ newCells[6*4+3] = newVertices[3];
+
+ // new cell 7
+ newCells[7*4+0] = cone[3]+coneVertexOffset;
+ newCells[7*4+1] = newVertices[3];
+ newCells[7*4+2] = newVertices[5];
+ newCells[7*4+3] = newVertices[4];
+
+ *numCells = numNewCells;
+ *cells = newCells;
+} // _newCells_TETRAHEDRON
+
+// ----------------------------------------------------------------------
+// Get new cells from refinement of a line cohseive cell with Lagrange
+// multiplier vertices.
+void
+ALE::CellRefinerTet4::_newCells_TRIANGLE_COHESIVE_LAGRANGE(const point_type** cells,
+ int *numCells,
+ const point_type cone[],
+ const int coneSize,
+ const int coneVertexOffsetNormal,
+ const int coneVertexOffsetCensored)
+{ // _newCells_TRIANGLE_COHESIVE_LAGRANGE
+ const int coneSizeTriPrism9 = 9;
+ const int numEdgesTriPrism9 = 9;
+ const int numNewCells = 4;
+ const int numNewVertices = 9;
+
+ int numEdges = 0;
+ const EdgeType *edges;
+ _edges_TRIANGLE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize);
+ assert(numEdgesTriPrism9 == numEdges);
+
+ static point_type newCells[numNewCells*coneSizeTriPrism9];
+ point_type newVertices[numNewVertices];
+ for(int iEdge=0, iNewVertex=0; iEdge < numEdgesTriPrism9; ++iEdge) {
+ if (_edgeToVertex.find(edges[iEdge]) == _edgeToVertex.end()) {
+ throw ALE::Exception("Missing edge in refined mesh");
+ } // if
+ newVertices[iNewVertex++] = _edgeToVertex[edges[iEdge]];
+ } // for
+
+ newCells[0*9+0] = cone[0]+coneVertexOffsetNormal; // New cell 0
+ newCells[0*9+1] = newVertices[0];
+ newCells[0*9+2] = newVertices[2];
+ newCells[0*9+3] = cone[3]+coneVertexOffsetNormal;
+ newCells[0*9+4] = newVertices[3];
+ newCells[0*9+5] = newVertices[5];
+ newCells[0*9+6] = cone[6]+coneVertexOffsetCensored;
+ newCells[0*9+7] = newVertices[6];
+ newCells[0*9+8] = newVertices[8];
+
+ newCells[1*9+0] = newVertices[0]; // New cell 1
+ newCells[1*9+1] = newVertices[1];
+ newCells[1*9+2] = newVertices[2];
+ newCells[1*9+3] = newVertices[3];
+ newCells[1*9+4] = newVertices[4];
+ newCells[1*9+5] = newVertices[5];
+ newCells[1*9+6] = newVertices[6];
+ newCells[1*9+7] = newVertices[7];
+ newCells[1*9+8] = newVertices[8];
+
+ newCells[2*9+0] = cone[1]+coneVertexOffsetNormal; // New cell 2
+ newCells[2*9+1] = newVertices[1];
+ newCells[2*9+2] = newVertices[0];
+ newCells[2*9+3] = cone[4]+coneVertexOffsetNormal;
+ newCells[2*9+4] = newVertices[4];
+ newCells[2*9+5] = newVertices[3];
+ newCells[2*9+6] = cone[7]+coneVertexOffsetCensored;
+ newCells[2*9+7] = newVertices[7];
+ newCells[2*9+8] = newVertices[6];
+
+ newCells[3*9+0] = cone[2]+coneVertexOffsetNormal; // New cell 3
+ newCells[3*9+1] = newVertices[2];
+ newCells[3*9+2] = newVertices[1];
+ newCells[3*9+3] = cone[5]+coneVertexOffsetNormal;
+ newCells[3*9+4] = newVertices[5];
+ newCells[3*9+5] = newVertices[4];
+ newCells[3*9+6] = cone[8]+coneVertexOffsetCensored;
+ newCells[3*9+7] = newVertices[8];
+ newCells[3*9+8] = newVertices[7];
+
+ *numCells = numNewCells;
+ *cells = newCells;
+} // _newCells_TRIANGLE_COHESIVE_LAGRANGE
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.hh 2010-09-27 21:06:53 UTC (rev 17223)
@@ -0,0 +1,232 @@
+// -*- 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/CellRefinerTet4.hh
+ *
+ * @brief Object for tet4 refinement of cells.
+ */
+
+#if !defined(pylith_topology_cellrefinertet4_hh)
+#define pylith_topology_cellrefinertet4_hh
+
+// Include directives ---------------------------------------------------
+#include "topologyfwd.hh" // forward declarations
+
+#include <list> // USES std::pair
+
+// CellRefinerTet4 ------------------------------------------------------
+/// Object for tet4 refinement of cells.
+class ALE::CellRefinerTet4
+{ // CellRefinerTet4
+ typedef IMesh<> mesh_type;
+ typedef mesh_type::point_type point_type;
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Constructor
+ *
+ * @param mesh Finite-element mesh.
+ */
+ CellRefinerTet4(const mesh_type& mesh);
+
+ /// Destructor
+ ~CellRefinerTet4(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);
+
+ /** Add space for new vertices in group.
+ *
+ * @param newGroup Group in refine mesh.
+ * @param oldGroup Group in original mesh.
+ */
+ void groupAddNewVertices(const ALE::Obj<mesh_type::int_section_type>& newGroup,
+ const ALE::Obj<mesh_type::int_section_type>& oldGroup);
+
+ /** Set new vertices in group.
+ *
+ * @param newGroup Group in refine mesh.
+ * @param oldGroup Group in original mesh.
+ */
+ void groupSetNewVertices(const ALE::Obj<mesh_type::int_section_type>& newGroup,
+ const ALE::Obj<mesh_type::int_section_type>& oldGroup);
+
+ /** Add new vertices to label.
+ *
+ * @param newMesh Mesh with refined cells.
+ * @param oldMesh Original mesh.
+ * @param labelName Name of label.
+ */
+ void labelAddNewVertices(const ALE::Obj<mesh_type>& newMesh,
+ const ALE::Obj<mesh_type>& oldMesh,
+ const char* labelName);
+
+// 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 {
+ TETRAHEDRON, // Normal tetrahedral cell
+ TRIANGLE_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_TETRAHEDRON(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_TRIANGLE_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_TETRAHEDRON(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 coneVertexOffsetNormal Offset for normal cone vertices.
+ * @param coneVertexOffset Offset for censored cone vertices.
+ */
+ void _newCells_TRIANGLE_COHESIVE_LAGRANGE(const point_type** cells,
+ int *numCells,
+ const point_type cone[],
+ const int coneSize,
+ const int coneVertexOffsetNormal,
+ const int coneVertexOffsetCensored);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ const mesh_type& _mesh;
+ edge_map_type _edgeToVertex;
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ CellRefinerTet4(const CellRefinerTet4&); ///< Not implemented
+ const CellRefinerTet4& operator=(const CellRefinerTet4&); ///< Not implemented
+
+}; // CellRefinerTet4
+
+#endif // pylith_topology_cellrefinertet4_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc 2010-09-27 19:32:24 UTC (rev 17222)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc 2010-09-27 21:06:53 UTC (rev 17223)
@@ -219,19 +219,20 @@
const Obj<mesh_type::label_type>& newLabel = newMesh->getLabel(labelName);
assert(!newLabel.isNull());
+ const int defaultValue = -999;
+
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;
-#if 0
- // NEED TO TEST IF OLD LABEL HAS ENDPOINTS
- if (oldLabel->supportContains<>(edgeVertexA) && oldLabel->supportContains<>(edgeVertexB)) {
- assert(oldMesh->getValue(oldLabel, edgeVertexA) == oldMesh->getValue(oldLabel, edgeVertexB));
- newMesh->setValue(newLabel, newVertex, oldMesh->getValue(oldLabel, edgeVertexA));
+ const int valueA = oldMesh->getValue(oldLabel, edgeVertexA, defaultValue);
+ const int valueB = oldMesh->getValue(oldLabel, edgeVertexB, defaultValue);
+
+ if (valueA != defaultValue && valueA == valueB) {
+ newMesh->setValue(newLabel, newVertex, valueA);
} // if
-#endif
} // for
} // labelAddNewVertices
Modified: short/3D/PyLith/trunk/libsrc/topology/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2010-09-27 19:32:24 UTC (rev 17222)
+++ short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2010-09-27 21:06:53 UTC (rev 17223)
@@ -49,7 +49,9 @@
# TEMPORARY
noinst_HEADERS += \
MeshRefiner.hh \
+ MeshRefiner.cc \
CellRefinerTri3.hh \
+ CellRefinerTet4.hh \
MeshOrder.hh
Modified: short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc 2010-09-27 19:32:24 UTC (rev 17222)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc 2010-09-27 21:06:53 UTC (rev 17223)
@@ -20,14 +20,14 @@
#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) :
+template<typename cellrefiner_type>
+ALE::MeshRefiner<cellrefiner_type>::MeshRefiner(void) :
_orderOldMesh(new MeshOrder),
_orderNewMesh(new MeshOrder)
{ // constructor
@@ -35,7 +35,8 @@
// ----------------------------------------------------------------------
// Destructor
-ALE::MeshRefiner::~MeshRefiner(void)
+template<typename cellrefiner_type>
+ALE::MeshRefiner<cellrefiner_type>::~MeshRefiner(void)
{ // destructor
delete _orderOldMesh; _orderOldMesh = PETSC_NULL;
delete _orderNewMesh; _orderNewMesh = PETSC_NULL;
@@ -43,10 +44,11 @@
// ----------------------------------------------------------------------
// Refine mesh.
+template<typename cellrefiner_type>
void
-ALE::MeshRefiner::refine(const Obj<mesh_type>& newMesh,
- const Obj<mesh_type>& mesh,
- CellRefinerTri3& refiner)
+ALE::MeshRefiner<cellrefiner_type>::refine(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ cellrefiner_type& refiner)
{ // refine
assert(!mesh.isNull());
if (mesh->hasLabel("censored depth")) {
@@ -58,10 +60,11 @@
// ----------------------------------------------------------------------
// Refine mesh.
+template<typename cellrefiner_type>
void
-ALE::MeshRefiner::_refine(const Obj<mesh_type>& newMesh,
- const Obj<mesh_type>& mesh,
- CellRefinerTri3& refiner)
+ALE::MeshRefiner<cellrefiner_type>::_refine(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ cellrefiner_type& refiner)
{ // _refine
typedef Interval<point_type>::const_iterator interval_type;
@@ -182,17 +185,18 @@
refiner.setCoordsNewVertices(newCoordinates, coordinates);
_stratify(newMesh);
- _calcNewOverlap(newMesh, mesh);
+ _calcNewOverlap(newMesh, mesh, refiner);
_createIntSections(newMesh, mesh, refiner);
_createLabels(newMesh, mesh, refiner);
} // _refine
// ----------------------------------------------------------------------
// Refine mesh with a censored depth.
- void
- ALE::MeshRefiner::_refineCensored(const Obj<mesh_type>& newMesh,
- const Obj<mesh_type>& mesh,
- CellRefinerTri3& refiner)
+template<typename cellrefiner_type>
+void
+ALE::MeshRefiner<cellrefiner_type>::_refineCensored(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ cellrefiner_type& refiner)
{ // _refineCensored
typedef Interval<point_type>::const_iterator interval_type;
@@ -385,7 +389,7 @@
refiner.setCoordsNewVertices(newCoordinates, coordinates);
_stratify(newMesh);
- _calcNewOverlap(newMesh, mesh);
+ _calcNewOverlap(newMesh, mesh, refiner);
_createIntSections(newMesh, mesh, refiner);
_createLabels(newMesh, mesh, refiner);
@@ -407,8 +411,9 @@
// ----------------------------------------------------------------------
// Stratify mesh.
+template<typename cellrefiner_type>
void
-ALE::MeshRefiner::_stratify(const Obj<mesh_type>& mesh)
+ALE::MeshRefiner<cellrefiner_type>::_stratify(const Obj<mesh_type>& mesh)
{ // _stratify
typedef Interval<point_type>::const_iterator interval_type;
@@ -449,10 +454,11 @@
// ----------------------------------------------------------------------
// Create integer sections in new mesh.
+template<typename cellrefiner_type>
void
-ALE::MeshRefiner::_createIntSections(const Obj<mesh_type>& newMesh,
- const Obj<mesh_type>& mesh,
- CellRefinerTri3& refiner)
+ALE::MeshRefiner<cellrefiner_type>::_createIntSections(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ cellrefiner_type& refiner)
{ // _createIntSections
assert(!newMesh.isNull());
assert(!mesh.isNull());
@@ -521,16 +527,19 @@
// ----------------------------------------------------------------------
// Create labels in new mesh.
+template<typename cellrefiner_type>
void
-ALE::MeshRefiner::_createLabels(const Obj<mesh_type>& newMesh,
- const Obj<mesh_type>& mesh,
- CellRefinerTri3& refiner)
+ALE::MeshRefiner<cellrefiner_type>::_createLabels(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ cellrefiner_type& refiner)
{ // _createLabels
assert(!newMesh.isNull());
assert(!mesh.isNull());
assert(_orderOldMesh);
assert(_orderNewMesh);
+ typedef ALE::Interval<point_type> interval_type;
+
const mesh_type::labels_type labels = mesh->getLabels();
const mesh_type::labels_type::const_iterator labelsEnd = labels.end();
for (mesh_type::labels_type::const_iterator l_iter=labels.begin(); l_iter != labelsEnd; ++l_iter) {
@@ -544,47 +553,68 @@
const Obj<mesh_type::label_type>& newLabel = newMesh->createLabel(l_iter->first);
assert(!newLabel.isNull());
- typedef mesh_type::label_type::traits::arrow_container_type::set_type arrows_set_type;
+ const int defaultValue = -999;
- arrows_set_type::const_iterator oldArrowsEnd = oldLabel->_arrows.set.end();
- for (arrows_set_type::const_iterator a_iter=oldLabel->_arrows.set.begin(); a_iter != oldArrowsEnd; ++a_iter) {
+ // Update cells
+ // Normal cells
+ interval_type::const_iterator pointsEnd = _orderOldMesh->cellsNormal().end();
+ for (interval_type::const_iterator p_iter = _orderOldMesh->cellsNormal().begin(); p_iter != pointsEnd; ++p_iter) {
+ const int value = mesh->getValue(oldLabel, *p_iter, defaultValue);
+ if (defaultValue == value)
+ continue;
- const mesh_type::point_type pOld = a_iter->target;
- const int value = mesh->getValue(oldLabel, pOld);
+ const int numNewCellsPerCell = refiner.numNewCells(*p_iter);
+ mesh_type::point_type pNew = _orderNewMesh->cellsNormal().min() + (*p_iter - _orderOldMesh->cellsNormal().min())*numNewCellsPerCell;
+ for(int i=0; i < numNewCellsPerCell; ++i, ++pNew)
+ newMesh->setValue(newLabel, pNew, value);
+ } // for
+
+ // Censored cells
+ pointsEnd = _orderOldMesh->cellsCensored().end();
+ for (interval_type::const_iterator p_iter = _orderOldMesh->cellsCensored().begin(); p_iter != pointsEnd; ++p_iter) {
+ const int value = mesh->getValue(oldLabel, *p_iter, defaultValue);
+ if (defaultValue == value)
+ continue;
- if (_orderOldMesh->cellsNormal().hasPoint(pOld)) {
- const int numNewCellsPerCell = refiner.numNewCells(pOld);
- mesh_type::point_type pNew = _orderNewMesh->cellsNormal().min() + (pOld - _orderOldMesh->cellsNormal().min())*numNewCellsPerCell;
- for(int i=0; i < numNewCellsPerCell; ++i, ++pNew)
- newMesh->setValue(newLabel, pNew, value);
-
- } else if (_orderOldMesh->verticesNormal().hasPoint(pOld)) {
- const mesh_type::point_type pNew = _orderNewMesh->verticesNormal().min() + (pOld - _orderOldMesh->verticesNormal().min());
+ const int numNewCellsPerCell = refiner.numNewCells(*p_iter);
+ mesh_type::point_type pNew = _orderNewMesh->cellsCensored().min() + (*p_iter - _orderOldMesh->cellsCensored().min())*numNewCellsPerCell;
+ for(int i=0; i < numNewCellsPerCell; ++i, ++pNew)
newMesh->setValue(newLabel, pNew, value);
-
- } else if (_orderOldMesh->verticesCensored().hasPoint(pOld)) {
- const mesh_type::point_type pNew = _orderNewMesh->verticesCensored().min() + (pOld - _orderOldMesh->verticesCensored().min());
- newMesh->setValue(newLabel, pNew, value);
-
- } else if (_orderOldMesh->cellsCensored().hasPoint(pOld)) {
- const int numNewCellsPerCell = refiner.numNewCells(pOld);
- mesh_type::point_type pNew = _orderNewMesh->cellsCensored().min() + (pOld - _orderOldMesh->cellsCensored().min())*numNewCellsPerCell;
- for(int i=0; i < numNewCellsPerCell; ++i, ++pNew)
- newMesh->setValue(newLabel, pNew, value);
- } else {
- throw ALE::Exception("Unexpected cell encountered when creating labels.");
- } // if/else
} // for
+
+ // Normal vertices
+ pointsEnd = _orderOldMesh->verticesNormal().end();
+ for (interval_type::const_iterator p_iter = _orderOldMesh->verticesNormal().begin(); p_iter != pointsEnd; ++p_iter) {
+ const int value = mesh->getValue(oldLabel, *p_iter, defaultValue);
+ if (defaultValue == value)
+ continue;
+
+ const mesh_type::point_type pNew = _orderNewMesh->verticesNormal().min() + (*p_iter - _orderOldMesh->verticesNormal().min());
+ newMesh->setValue(newLabel, pNew, value);
+ } // for
+ // Censored vertices
+ pointsEnd = _orderOldMesh->verticesCensored().end();
+ for (interval_type::const_iterator p_iter = _orderOldMesh->verticesCensored().begin(); p_iter != pointsEnd; ++p_iter) {
+ const int value = mesh->getValue(oldLabel, *p_iter, defaultValue);
+ if (defaultValue == value)
+ continue;
+
+ const mesh_type::point_type pNew = _orderNewMesh->verticesCensored().min() + (*p_iter - _orderOldMesh->verticesCensored().min());
+ newMesh->setValue(newLabel, pNew, value);
+ } // for
+
refiner.labelAddNewVertices(newMesh, mesh, l_iter->first.c_str());
} // for
} // _createLabels
// ----------------------------------------------------------------------
// Calculate new overlap.
+template<typename cellrefiner_type>
void
-ALE::MeshRefiner::_calcNewOverlap(const Obj<mesh_type>& newMesh,
- const Obj<mesh_type>& mesh)
+ALE::MeshRefiner<cellrefiner_type>::_calcNewOverlap(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ cellrefiner_type& refiner)
{ // _calcNewOverlap
assert(!newMesh.isNull());
assert(!mesh.isNull());
Modified: short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh 2010-09-27 19:32:24 UTC (rev 17222)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh 2010-09-27 21:06:53 UTC (rev 17223)
@@ -19,7 +19,7 @@
/**
* @file libsrc/topology/MeshRefiner.hh
*
- * @brief Object for tri3 refinement of cells.
+ * @brief Object for refinement of cells.
*/
#if !defined(pylith_topology_meshrefiner_hh)
@@ -29,7 +29,8 @@
#include "topologyfwd.hh" // forward declarations
// RefineTri3 --------------------------------------------------------
-/// Object for tri3 refinement of cells.
+/// Object for refinement of cells.
+template<typename cellrefiner_type>
class ALE::MeshRefiner
{ // MeshRefiner
typedef IMesh<> mesh_type;
@@ -55,7 +56,7 @@
*/
void refine(const Obj<mesh_type>& newMesh,
const Obj<mesh_type>& mesh,
- CellRefinerTri3& refiner);
+ cellrefiner_type& refiner);
// PRIVATE METHODS //////////////////////////////////////////////////////
private :
@@ -68,7 +69,7 @@
*/
void _refine(const Obj<mesh_type>& newMesh,
const Obj<mesh_type>& mesh,
- CellRefinerTri3& refiner);
+ cellrefiner_type& refiner);
/** Refine mesh with a censored depth.
*
@@ -78,7 +79,7 @@
*/
void _refineCensored(const Obj<mesh_type>& newMesh,
const Obj<mesh_type>& mesh,
- CellRefinerTri3& refiner);
+ cellrefiner_type& refiner);
/** Stratify mesh.
*
@@ -92,10 +93,11 @@
*
* @param newMesh New (refined) mesh.
* @param mesh Current (unrefined) mesh with integer sections.
+ * @param refiner Cell refiner.
*/
void _createIntSections(const Obj<mesh_type>& newMesh,
const Obj<mesh_type>& mesh,
- CellRefinerTri3& refiner);
+ cellrefiner_type& refiner);
/** Create labels in new mesh.
*
@@ -103,18 +105,21 @@
*
* @param newMesh New (refined) mesh.
* @param mesh Current (unrefined) mesh with integer sections.
+ * @param refiner Cell refiner.
*/
void _createLabels(const Obj<mesh_type>& newMesh,
const Obj<mesh_type>& mesh,
- CellRefinerTri3& refiner);
+ cellrefiner_type& refiner);
/** Calculate new overlap.
*
* @param newMesh New (refined) mesh.
* @param mesh Current (unrefined) mesh with overlap.
+ * @param refiner Cell refiner.
*/
void _calcNewOverlap(const Obj<mesh_type>& newMesh,
- const Obj<mesh_type>& mesh);
+ const Obj<mesh_type>& mesh,
+ cellrefiner_type& refiner);
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
@@ -130,6 +135,8 @@
}; // MeshRefiner
+#include "MeshRefiner.cc"
+
#endif // pylith_topology_meshrefiner_hh
Modified: short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc 2010-09-27 19:32:24 UTC (rev 17222)
+++ short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc 2010-09-27 21:06:53 UTC (rev 17223)
@@ -23,6 +23,7 @@
#include "Mesh.hh" // USES Mesh
#include "CellRefinerTri3.hh" // USES CellRefinerTri3
+#include "CellRefinerTet4.hh" // USES CellRefinerTet4
#include "MeshRefiner.hh" // USES MeshRefiner
#include <stdexcept> // USES std::runtime_error
@@ -63,10 +64,56 @@
new SieveMesh::sieve_type(mesh.comm(), mesh.debug());
newSieveMesh->setSieve(newSieve);
- ALE::CellRefinerTri3 cellSplitter(*sieveMesh);
- ALE::MeshRefiner refiner;
- refiner.refine(newSieveMesh, sieveMesh, cellSplitter);
+ const Obj<SieveMesh::label_sequence>& cells = sieveMesh->heightStratum(0);
+ assert(!cells.isNull());
+ assert(cells->size() > 0);
+ const int numCorners = sieveMesh->getNumCellCorners();
+ const int dim = mesh.dimension();
+
+ switch (dim) {
+ case 0:
+ case 1:
+ throw std::runtime_error("Uniform refinement not implemented.");
+ break;
+
+ case 2:
+ switch (numCorners) {
+ case 3: {
+ ALE::CellRefinerTri3 cellSplitter(*sieveMesh);
+ ALE::MeshRefiner<ALE::CellRefinerTri3> refinement;
+ refinement.refine(newSieveMesh, sieveMesh, cellSplitter);
+ break;
+ } // case 3
+ case 4: {
+ throw std::logic_error("Not implemented.");
+ break;
+ } // case 4
+ default :
+ throw std::runtime_error("Unknown number of corners for cells.");
+ } // switch
+ break;
+
+ case 3:
+ switch (numCorners) {
+ case 4: {
+ ALE::CellRefinerTet4 cellSplitter(*sieveMesh);
+ ALE::MeshRefiner<ALE::CellRefinerTet4> refinement;
+ refinement.refine(newSieveMesh, sieveMesh, cellSplitter);
+ break;
+ } // case 4
+ case 8: {
+ throw std::logic_error("Not implemented.");
+ break;
+ } // case 4
+ default :
+ throw std::runtime_error("Unknown number of corners for cells.");
+ } // switch
+ break;
+
+ default :
+ throw std::logic_error("Unknown dimension.");
+ } // switch
} // refine
Modified: short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh 2010-09-27 19:32:24 UTC (rev 17222)
+++ short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh 2010-09-27 21:06:53 UTC (rev 17223)
@@ -34,8 +34,9 @@
typename value_type,
typename allocator> class IUniformSectionDS;
- class MeshRefiner;
+ template<typename cellrefiner_type> class MeshRefiner;
class CellRefinerTri3;
+ class CellRefinerTet4;
class MeshOrder;
} // ALE
More information about the CIG-COMMITS
mailing list