[cig-commits] r17243 - in short/3D/PyLith/trunk/libsrc: . topology
brad at geodynamics.org
brad at geodynamics.org
Fri Oct 8 09:41:39 PDT 2010
Author: brad
Date: 2010-10-08 09:41:39 -0700 (Fri, 08 Oct 2010)
New Revision: 17243
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.cc
short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.hh
short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc
short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh
short/3D/PyLith/trunk/libsrc/topology/Makefile.am
short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh
Log:
Refactored tri3 and tet4 refinement. Use common base class RefineEdges2 for edge refinement data structure and common operations.
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2010-10-07 14:59:43 UTC (rev 17242)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2010-10-08 16:41:39 UTC (rev 17243)
@@ -143,7 +143,8 @@
libpylith_la_SOURCES += \
topology/CellRefinerTri3.cc \
topology/CellRefinerTet4.cc \
- topology/MeshOrder.cc
+ topology/MeshOrder.cc \
+ topology/RefineEdges2.cc
libpylith_la_LDFLAGS = $(AM_LDFLAGS) $(PYTHON_LA_LDFLAGS)
Modified: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.cc 2010-10-07 14:59:43 UTC (rev 17242)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.cc 2010-10-08 16:41:39 UTC (rev 17243)
@@ -28,7 +28,7 @@
// ----------------------------------------------------------------------
// Constructor
ALE::CellRefinerTet4::CellRefinerTet4(const mesh_type& mesh) :
- _mesh(mesh)
+ RefineEdges2(mesh)
{ // constructor
assert(3 == mesh.getDimension());
} // constructor
@@ -157,241 +157,6 @@
} // 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
-
-// ----------------------------------------------------------------------
-// Calculate new overlap.
-void
-ALE::CellRefinerTet4::overlapAddNewVertices(const Obj<mesh_type>& newMesh,
- const MeshOrder& orderNewMesh,
- const Obj<mesh_type>& oldMesh,
- const MeshOrder& orderOldMesh)
-{ // overlapAddNewVertices
- assert(!newMesh.isNull());
- assert(!oldMesh.isNull());
-
- Obj<mesh_type::send_overlap_type> newSendOverlap = newMesh->getSendOverlap();
- assert(!newSendOverlap.isNull());
- Obj<mesh_type::recv_overlap_type> newRecvOverlap = newMesh->getRecvOverlap();
- assert(!newRecvOverlap.isNull());
- const Obj<mesh_type::send_overlap_type>& oldSendOverlap = oldMesh->getSendOverlap();
- assert(!oldSendOverlap.isNull());
- const Obj<mesh_type::recv_overlap_type>& oldRecvOverlap = oldMesh->getRecvOverlap();
- assert(!oldRecvOverlap.isNull());
-
-
- // Check edges in edgeToVertex 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, EdgeType> > newVerticesSection = new ALE::Section<point_type, EdgeType>(oldMesh->comm());
- assert(!newVerticesSection.isNull());
- std::map<EdgeType, std::vector<int> > bndryEdgeToRank;
-
- const int localOffset = orderNewMesh.verticesNormal().min() - orderOldMesh.verticesNormal().min();
-
- for(std::map<EdgeType, point_type>::const_iterator e_iter = _edgeToVertex.begin(); e_iter != _edgeToVertex.end(); ++e_iter) {
- const point_type left = e_iter->first.first;
- const point_type right = e_iter->first.second;
-
- if (oldSendOverlap->capContains(left) && oldSendOverlap->capContains(right)) {
- const Obj<mesh_type::send_overlap_type::traits::supportSequence>& leftRanksSeq = oldSendOverlap->support(left);
- assert(!leftRanksSeq.isNull());
- std::list<int> leftRanks(leftRanksSeq->begin(), leftRanksSeq->end());
- const Obj<mesh_type::send_overlap_type::traits::supportSequence>& rightRanks = oldSendOverlap->support(right);
- assert(!rightRanks.isNull());
- 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) {
- bndryEdgeToRank[e_iter->first].push_back(*r_iter);
- } // for
- } // if
- } // if
- } // for
- newVerticesSection->allocatePoint();
- const ALE::Section<point_type, EdgeType>::chart_type& chart = newVerticesSection->getChart();
-
- for(ALE::Section<point_type, EdgeType>::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
- typedef ALE::Section<point_type, EdgeType>::value_type value_type;
- const point_type p = *c_iter;
- const int dim = newVerticesSection->getFiberDimension(p);
- int v = 0;
- value_type* values = (dim > 0) ? new value_type[dim] : 0;
-
- for(std::map<EdgeType, std::vector<int> >::const_iterator e_iter = bndryEdgeToRank.begin(); e_iter != bndryEdgeToRank.end() && v < dim; ++e_iter) {
- if (std::min(e_iter->first.first, e_iter->first.second)+localOffset == p) {
- values[v++] = EdgeType(std::max(e_iter->first.first, e_iter->first.second)+localOffset, _edgeToVertex[e_iter->first]);
- } // if
- } // for
- newVerticesSection->updatePoint(p, values);
- delete [] values;
- } // for
- // Copy across overlap
- typedef ALE::Pair<int, point_type> overlap_point_type;
- Obj<ALE::Section<overlap_point_type, EdgeType> > overlapVertices = new ALE::Section<overlap_point_type, EdgeType>(oldMesh->comm());
-
- ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newVerticesSection, overlapVertices);
- // Merge by translating edge to local points, finding edge in _edgeToVertex, and adding (local new vetex, remote new vertex) to overlap
- for(std::map<EdgeType, std::vector<int> >::const_iterator e_iter = bndryEdgeToRank.begin(); e_iter != bndryEdgeToRank.end(); ++e_iter) {
- const point_type localPoint = _edgeToVertex[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;
- } // if
- } // for
- 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;
- } // if
- } // for
- 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 EdgeType *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;
- } // if
- } // for
- newSendOverlap->addArrow(localPoint, rank, remotePoint);
- newRecvOverlap->addArrow(rank, localPoint, remotePoint);
- } // for
- } // for
-
- oldSendOverlap->view("OLD SEND OVERLAP");
- oldRecvOverlap->view("OLD RECV OVERLAP");
- newSendOverlap->view("NEW SEND OVERLAP");
- newRecvOverlap->view("NEW RECV OVERLAP");
-} // overlapAddNewVertices
-
-// ----------------------------------------------------------------------
// Get cell type.
ALE::CellRefinerTet4::CellEnum
ALE::CellRefinerTet4::_cellType(const point_type cell)
Modified: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.hh 2010-10-07 14:59:43 UTC (rev 17242)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.hh 2010-10-08 16:41:39 UTC (rev 17243)
@@ -26,17 +26,12 @@
#define pylith_topology_cellrefinertet4_hh
// Include directives ---------------------------------------------------
-#include "topologyfwd.hh" // forward declarations
+#include "RefineEdges2.hh" // ISA RefineEdges2
-#include <list> // USES std::pair
-
// CellRefinerTet4 ------------------------------------------------------
/// Object for tet4 refinement of cells.
-class ALE::CellRefinerTet4
+class ALE::CellRefinerTet4 : public RefineEdges2
{ // CellRefinerTet4
- typedef IMesh<> mesh_type;
- typedef mesh_type::point_type point_type;
-
// PUBLIC MEMBERS ///////////////////////////////////////////////////////
public :
@@ -100,71 +95,6 @@
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);
-
- /** Calculate new overlap.
- *
- * @param newMesh New (refined) mesh.
- * @param orderNewMesh Order in new mesh.
- * @param oldMesh Current (unrefined) mesh with overlap.
- * @param orderOldMesh Order in old mesh.
- */
- void overlapAddNewVertices(const Obj<mesh_type>& newMesh,
- const MeshOrder& orderNewMesh,
- const Obj<mesh_type>& oldMesh,
- const MeshOrder& orderOldMesh);
-
-// 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 :
@@ -240,12 +170,6 @@
const int coneVertexOffsetNormal,
const int coneVertexOffsetCensored);
-// PRIVATE MEMBERS //////////////////////////////////////////////////////
-private :
-
- const mesh_type& _mesh;
- edge_map_type _edgeToVertex;
-
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc 2010-10-07 14:59:43 UTC (rev 17242)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc 2010-10-08 16:41:39 UTC (rev 17243)
@@ -28,7 +28,7 @@
// ----------------------------------------------------------------------
// Constructor
ALE::CellRefinerTri3::CellRefinerTri3(const mesh_type& mesh) :
- _mesh(mesh)
+ RefineEdges2(mesh)
{ // constructor
assert(2 == mesh.getDimension());
} // constructor
@@ -157,241 +157,6 @@
} // 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
-
-// ----------------------------------------------------------------------
-// Add space for new vertices in group.
-void
-ALE::CellRefinerTri3::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::CellRefinerTri3::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::CellRefinerTri3::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
-
-// ----------------------------------------------------------------------
-// Calculate new overlap.
-void
-ALE::CellRefinerTri3::overlapAddNewVertices(const Obj<mesh_type>& newMesh,
- const MeshOrder& orderNewMesh,
- const Obj<mesh_type>& oldMesh,
- const MeshOrder& orderOldMesh)
-{ // overlapAddNewVertices
- assert(!newMesh.isNull());
- assert(!oldMesh.isNull());
-
- Obj<mesh_type::send_overlap_type> newSendOverlap = newMesh->getSendOverlap();
- assert(!newSendOverlap.isNull());
- Obj<mesh_type::recv_overlap_type> newRecvOverlap = newMesh->getRecvOverlap();
- assert(!newRecvOverlap.isNull());
- const Obj<mesh_type::send_overlap_type>& oldSendOverlap = oldMesh->getSendOverlap();
- assert(!oldSendOverlap.isNull());
- const Obj<mesh_type::recv_overlap_type>& oldRecvOverlap = oldMesh->getRecvOverlap();
- assert(!oldRecvOverlap.isNull());
-
-
- // Check edges in edgeToVertex 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, EdgeType> > newVerticesSection = new ALE::Section<point_type, EdgeType>(oldMesh->comm());
- assert(!newVerticesSection.isNull());
- std::map<EdgeType, std::vector<int> > bndryEdgeToRank;
-
- const int localOffset = orderNewMesh.verticesNormal().min() - orderOldMesh.verticesNormal().min();
-
- for(std::map<EdgeType, point_type>::const_iterator e_iter = _edgeToVertex.begin(); e_iter != _edgeToVertex.end(); ++e_iter) {
- const point_type left = e_iter->first.first;
- const point_type right = e_iter->first.second;
-
- if (oldSendOverlap->capContains(left) && oldSendOverlap->capContains(right)) {
- const Obj<mesh_type::send_overlap_type::traits::supportSequence>& leftRanksSeq = oldSendOverlap->support(left);
- assert(!leftRanksSeq.isNull());
- std::list<int> leftRanks(leftRanksSeq->begin(), leftRanksSeq->end());
- const Obj<mesh_type::send_overlap_type::traits::supportSequence>& rightRanks = oldSendOverlap->support(right);
- assert(!rightRanks.isNull());
- 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) {
- bndryEdgeToRank[e_iter->first].push_back(*r_iter);
- } // for
- } // if
- } // if
- } // for
- newVerticesSection->allocatePoint();
- const ALE::Section<point_type, EdgeType>::chart_type& chart = newVerticesSection->getChart();
-
- for(ALE::Section<point_type, EdgeType>::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
- typedef ALE::Section<point_type, EdgeType>::value_type value_type;
- const point_type p = *c_iter;
- const int dim = newVerticesSection->getFiberDimension(p);
- int v = 0;
- value_type* values = (dim > 0) ? new value_type[dim] : 0;
-
- for(std::map<EdgeType, std::vector<int> >::const_iterator e_iter = bndryEdgeToRank.begin(); e_iter != bndryEdgeToRank.end() && v < dim; ++e_iter) {
- if (std::min(e_iter->first.first, e_iter->first.second)+localOffset == p) {
- values[v++] = EdgeType(std::max(e_iter->first.first, e_iter->first.second)+localOffset, _edgeToVertex[e_iter->first]);
- } // if
- } // for
- newVerticesSection->updatePoint(p, values);
- delete [] values;
- } // for
- // Copy across overlap
- typedef ALE::Pair<int, point_type> overlap_point_type;
- Obj<ALE::Section<overlap_point_type, EdgeType> > overlapVertices = new ALE::Section<overlap_point_type, EdgeType>(oldMesh->comm());
-
- ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newVerticesSection, overlapVertices);
- // Merge by translating edge to local points, finding edge in _edgeToVertex, and adding (local new vetex, remote new vertex) to overlap
- for(std::map<EdgeType, std::vector<int> >::const_iterator e_iter = bndryEdgeToRank.begin(); e_iter != bndryEdgeToRank.end(); ++e_iter) {
- const point_type localPoint = _edgeToVertex[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;
- } // if
- } // for
- 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;
- } // if
- } // for
- 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 EdgeType *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;
- } // if
- } // for
- newSendOverlap->addArrow(localPoint, rank, remotePoint);
- newRecvOverlap->addArrow(rank, localPoint, remotePoint);
- } // for
- } // for
-
- oldSendOverlap->view("OLD SEND OVERLAP");
- oldRecvOverlap->view("OLD RECV OVERLAP");
- newSendOverlap->view("NEW SEND OVERLAP");
- newRecvOverlap->view("NEW RECV OVERLAP");
-} // overlapAddNewVertces
-
-// ----------------------------------------------------------------------
// Get cell type.
ALE::CellRefinerTri3::CellEnum
ALE::CellRefinerTri3::_cellType(const point_type cell)
Modified: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh 2010-10-07 14:59:43 UTC (rev 17242)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh 2010-10-08 16:41:39 UTC (rev 17243)
@@ -26,17 +26,12 @@
#define pylith_topology_cellrefinertri3_hh
// Include directives ---------------------------------------------------
-#include "topologyfwd.hh" // forward declarations
+#include "RefineEdges2.hh" // ISA RefineEdges2
-#include <list> // USES std::pair
-
// CellRefinerTri3 ------------------------------------------------------
/// Object for tri3 refinement of cells.
-class ALE::CellRefinerTri3
+class ALE::CellRefinerTri3 : public RefineEdges2
{ // CellRefinerTri3
- typedef IMesh<> mesh_type;
- typedef mesh_type::point_type point_type;
-
// PUBLIC MEMBERS ///////////////////////////////////////////////////////
public :
@@ -100,71 +95,6 @@
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);
-
- /** Calculate new overlap.
- *
- * @param newMesh New (refined) mesh.
- * @param orderNewMesh Order in new mesh.
- * @param oldMesh Current (unrefined) mesh with overlap.
- * @param orderOldMesh Order in old mesh.
- */
- void overlapAddNewVertices(const Obj<mesh_type>& newMesh,
- const MeshOrder& orderNewMesh,
- const Obj<mesh_type>& oldMesh,
- const MeshOrder& orderOldMesh);
-
-// 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 :
@@ -240,12 +170,6 @@
const int coneVertexOffsetNormal,
const int coneVertexOffsetCensored);
-// PRIVATE MEMBERS //////////////////////////////////////////////////////
-private :
-
- const mesh_type& _mesh;
- edge_map_type _edgeToVertex;
-
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/topology/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2010-10-07 14:59:43 UTC (rev 17242)
+++ short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2010-10-08 16:41:39 UTC (rev 17243)
@@ -52,7 +52,8 @@
MeshRefiner.cc \
CellRefinerTri3.hh \
CellRefinerTet4.hh \
- MeshOrder.hh
+ MeshOrder.hh \
+ RefineEdges2.hh
# export
Modified: short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh 2010-10-07 14:59:43 UTC (rev 17242)
+++ short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh 2010-10-08 16:41:39 UTC (rev 17243)
@@ -35,6 +35,7 @@
typename allocator> class IUniformSectionDS;
template<typename cellrefiner_type> class MeshRefiner;
+ class RefineEdges2;
class CellRefinerTri3;
class CellRefinerTet4;
class MeshOrder;
More information about the CIG-COMMITS
mailing list