[cig-commits] r17226 - in short/3D/PyLith/trunk: libsrc/topology unittests/libtests/topology
brad at geodynamics.org
brad at geodynamics.org
Mon Sep 27 18:09:19 PDT 2010
Author: brad
Date: 2010-09-27 18:09:18 -0700 (Mon, 27 Sep 2010)
New Revision: 17226
Modified:
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/MeshRefiner.cc
short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc
Log:
Fixed bug in parallel global uniform refinement. Processor boundary may be cohesive cell, so we must split cohesive cells without adding censored (Lagrange multiplier) vertices and then split them allowing Lagrange multiplier vertices.
Modified: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.cc 2010-09-27 22:25:15 UTC (rev 17225)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.cc 2010-09-28 01:09:18 UTC (rev 17226)
@@ -90,6 +90,42 @@
} // splitCell
// ----------------------------------------------------------------------
+// Split cell into smaller cells of same type.
+void
+ALE::CellRefinerTet4::splitCellUncensored(const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ point_type* curNewVertex)
+{ // splitCellUncensored
+ assert(curNewVertex);
+
+ int numEdges = 0;
+ const EdgeType* edges;
+
+ const bool uncensored = true;
+
+ switch (_cellType(cell)) {
+ case TETRAHEDRON:
+ // No censored vertices on normal cell.
+ break;
+ case TRIANGLE_COHESIVE_LAGRANGE:
+ _edges_TRIANGLE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize, uncensored);
+ 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
+} // splitCellUncensored
+
+// ----------------------------------------------------------------------
// Get refined cells.
void
ALE::CellRefinerTet4::getNewCells(const point_type** cells,
@@ -237,6 +273,16 @@
} // 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
+} // overlapAddNewVertices
+
+// ----------------------------------------------------------------------
// Get cell type.
ALE::CellRefinerTet4::CellEnum
ALE::CellRefinerTet4::_cellType(const point_type cell)
@@ -290,22 +336,39 @@
ALE::CellRefinerTet4::_edges_TRIANGLE_COHESIVE_LAGRANGE(const EdgeType** edges,
int* numEdges,
const point_type cone[],
- const int coneSize)
+ const int coneSize,
+ const bool uncensored)
{ // _edges_TRIANGLE_COHESIVE_LAGRANGE
- static EdgeType splitEdges[9];
+ if (uncensored) {
+ // Use all vertices
+ 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;
+ 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;
+ } else {
+ // Omit edges with censored (Lagrange multipler) vertices.
+ static EdgeType splitEdges[6];
+
+ 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]));
+ *numEdges = 6;
+ *edges = splitEdges;
+ } // if/else
} // _edges_TRIANGLE_COHESIVE_LAGRANGE
// ----------------------------------------------------------------------
@@ -406,7 +469,7 @@
int numEdges = 0;
const EdgeType *edges;
- _edges_TRIANGLE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize);
+ _edges_TRIANGLE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize, true);
assert(numEdgesTriPrism9 == numEdges);
static point_type newCells[numNewCells*coneSizeTriPrism9];
Modified: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.hh 2010-09-27 22:25:15 UTC (rev 17225)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTet4.hh 2010-09-28 01:09:18 UTC (rev 17226)
@@ -57,7 +57,8 @@
*/
int numNewCells(const point_type cell);
- /** Split cell into smaller cells of same type.
+ /** Split cell into smaller cells of same type. Do not create
+ * censored vertices on censored cells.
*
* @param cell Original cell.
* @param cone Vertices in cell (original mesh).
@@ -69,6 +70,19 @@
const int coneSize,
point_type* curNewVertex);
+ /** Split cell into smaller cells of same type. Create only censored
+ * vertices on censored cells.
+ *
+ * @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 splitCellUncensored(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).
@@ -120,6 +134,18 @@
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 :
@@ -175,11 +201,13 @@
* @param numEdges Number of edges.
* @param cone Vertices in cell (original mesh).
* @param coneSize Number of vertices in cell.
+ * @param uncensored True if including edges with censored vertices.
*/
void _edges_TRIANGLE_COHESIVE_LAGRANGE(const EdgeType** edges,
int* numEdges,
const point_type cone[],
- const int coneSize);
+ const int coneSize,
+ const bool uncensored =false);
/** Get new cells from refinement of a triangular cell.
*
Modified: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc 2010-09-27 22:25:15 UTC (rev 17225)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc 2010-09-28 01:09:18 UTC (rev 17226)
@@ -90,6 +90,42 @@
} // splitCell
// ----------------------------------------------------------------------
+// Split cell into smaller cells of same type.
+void
+ALE::CellRefinerTri3::splitCellUncensored(const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ point_type* curNewVertex)
+{ // splitCellUncensored
+ assert(curNewVertex);
+
+ int numEdges = 0;
+ const EdgeType* edges;
+
+ const bool uncensored = true;
+
+ switch (_cellType(cell)) {
+ case TRIANGLE:
+ // No censored vertices on normal cell.
+ break;
+ case LINE_COHESIVE_LAGRANGE:
+ _edges_LINE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize, uncensored);
+ 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
+} // splitCellUncensored
+
+// ----------------------------------------------------------------------
// Get refined cells.
void
ALE::CellRefinerTri3::getNewCells(const point_type** cells,
@@ -237,6 +273,125 @@
} // 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)
@@ -287,16 +442,29 @@
ALE::CellRefinerTri3::_edges_LINE_COHESIVE_LAGRANGE(const EdgeType** edges,
int* numEdges,
const point_type cone[],
- const int coneSize)
+ const int coneSize,
+ const bool uncensored)
{ // _edges_LINE_COHESIVE_LAGRANGE
- static EdgeType lineEdges[6];
+ if (uncensored) {
+ // Include all edges
+ static EdgeType lineEdges[3];
- 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;
+ 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;
+ } else {
+ // Omit edges with censored (Lagrange multiplier) vertices.
+ static EdgeType lineEdges[2];
+
+ 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]));
+ *numEdges = 2;
+ *edges = lineEdges;
+ } // if/else
} // _edges_LINE_COHESIVE_LAGRANGE
// ----------------------------------------------------------------------
@@ -369,7 +537,7 @@
int numEdges = 0;
const EdgeType *edges;
- _edges_LINE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize);
+ _edges_LINE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize, true);
assert(numEdgesLine6 == numEdges);
static point_type lineCells[numNewCells*coneSizeLine6];
Modified: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh 2010-09-27 22:25:15 UTC (rev 17225)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh 2010-09-28 01:09:18 UTC (rev 17226)
@@ -57,7 +57,8 @@
*/
int numNewCells(const point_type cell);
- /** Split cell into smaller cells of same type.
+ /** Split cell into smaller cells of same type. Do not create
+ * censored vertices on censored cells.
*
* @param cell Original cell.
* @param cone Vertices in cell (original mesh).
@@ -69,6 +70,19 @@
const int coneSize,
point_type* curNewVertex);
+ /** Split cell into smaller cells of same type. Create only censored
+ * vertices on censored cells.
+ *
+ * @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 splitCellUncensored(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).
@@ -120,6 +134,18 @@
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 :
@@ -175,11 +201,13 @@
* @param numEdges Number of edges.
* @param cone Vertices in cell (original mesh).
* @param coneSize Number of vertices in cell.
+ * @param uncensored True if including edges with censored vertices.
*/
void _edges_LINE_COHESIVE_LAGRANGE(const EdgeType** edges,
int* numEdges,
const point_type cone[],
- const int coneSize);
+ const int coneSize,
+ const bool uncensored =false);
/** Get new cells from refinement of a triangular cell.
*
Modified: short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc 2010-09-27 22:25:15 UTC (rev 17225)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc 2010-09-28 01:09:18 UTC (rev 17226)
@@ -231,9 +231,10 @@
// Count number of remaining cells (other cells).
const int numSkip = cellsNormal->size();
- mesh_type::label_sequence::iterator c_iter = cells->begin();
+ mesh_type::label_sequence::const_iterator c_iter = cells->begin();
for (int i=0; i < numSkip; ++i)
++c_iter;
+ const mesh_type::label_sequence::const_iterator cellsCensoredBegin = c_iter;
int newNumCellsCensored = 0;
for (; c_iter != cellsEnd; ++c_iter)
newNumCellsCensored += refiner.numNewCells(*c_iter);
@@ -251,19 +252,21 @@
sieve->cone(*c_iter, cV);
refiner.splitCell(*c_iter, cV.getPoints(), cV.getSize(), &curNewVertex);
} // for
+ for(mesh_type::label_sequence::iterator c_iter = cellsCensoredBegin; c_iter != cellsEnd; ++c_iter) {
+ cV.clear();
+ sieve->cone(*c_iter, cV);
+ refiner.splitCell(*c_iter, cV.getPoints(), cV.getSize(), &curNewVertex);
+ } // for
const int newNumVerticesNormal = curNewVertex - counterBegin + oldNumVerticesNormal;
// Count number of remaining vertices (other vertices).
const int oldNumVerticesCensored = vertices->size() - oldNumVerticesNormal;
counterBegin = newNumCellsNormal + newNumVerticesNormal + oldNumVerticesCensored;
curNewVertex = counterBegin;
- c_iter = cells->begin();
- for (int i=0; i < numSkip; ++i)
- ++c_iter;
- for (; c_iter != cellsEnd; ++c_iter) {
+ for(mesh_type::label_sequence::iterator c_iter = cellsCensoredBegin; c_iter != cellsEnd; ++c_iter) {
cV.clear();
sieve->cone(*c_iter, cV);
- refiner.splitCell(*c_iter, cV.getPoints(), cV.getSize(), &curNewVertex);
+ refiner.splitCellUncensored(*c_iter, cV.getPoints(), cV.getSize(), &curNewVertex);
} // for
const int newNumVerticesCensored = curNewVertex - counterBegin + oldNumVerticesCensored;
@@ -627,19 +630,23 @@
// 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;
+
+ // Get offsets for points across processors and add points in overlap from original mesh to new overlap
+ int* oldVerticesStartP = (mesh->commSize() > 0) ? new int[mesh->commSize()] : 0;
+ int* newVerticesStartP = (newMesh->commSize() > 0) ? new int[newMesh->commSize()] : 0;
+ int ierr = 0;
- 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);
+ const int oldVerticesStart = _orderOldMesh->verticesNormal().min();
+ const int newVerticesStart = _orderNewMesh->verticesNormal().min();
+
+ ierr = MPI_Allgather((void *) &oldVerticesStart, 1, MPI_INT, oldVerticesStartP, 1, MPI_INT, mesh->comm());CHKERRXX(ierr);
+ ierr = MPI_Allgather((void *) &newVerticesStart, 1, MPI_INT, newVerticesStartP, 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()];
+ const mesh_type::send_overlap_type::source_type localOffset = newVerticesStart - oldVerticesStart;
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);
@@ -648,11 +655,11 @@
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];
+ const mesh_type::send_overlap_type::source_type remoteOffset = newVerticesStartP[rank] - oldVerticesStartP[rank];
newSendOverlap->addArrow(localPoint+localOffset, rank, remotePoint+remoteOffset);
- }
- }
+ } // for
+ } // for
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) {
@@ -662,101 +669,16 @@
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];
+ const mesh_type::recv_overlap_type::target_type remoteOffset = newVerticesStartP[rank] - oldVerticesStartP[rank];
newRecvOverlap->addArrow(rank, localPoint+localOffset, remotePoint+remoteOffset);
- }
- }
+ } // for
+ } // for
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
+ delete [] oldVerticesStartP; oldVerticesStartP = PETSC_NULL;
+ delete [] newVerticesStartP; newVerticesStartP = PETSC_NULL;
+
+ refiner.overlapAddNewVertices(newMesh, *_orderNewMesh, mesh, *_orderOldMesh);
} // _calcNewOverlap
Modified: short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc 2010-09-27 22:25:15 UTC (rev 17225)
+++ short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc 2010-09-28 01:09:18 UTC (rev 17226)
@@ -114,6 +114,8 @@
default :
throw std::logic_error("Unknown dimension.");
} // switch
+
+ newMesh->view("REFINED MESH");
} // refine
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc 2010-09-27 22:25:15 UTC (rev 17225)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc 2010-09-28 01:09:18 UTC (rev 17226)
@@ -130,8 +130,6 @@
faultB.adjustTopology(mesh, &firstFaultVertex,
&firstLagrangeVertex, &firstFaultCell);
} // if
-
- mesh->view("UNREFINED MESH");
} // _setupMesh
// ----------------------------------------------------------------------
@@ -148,8 +146,6 @@
Mesh newMesh(data.cellDim);
refiner.refine(&newMesh, mesh, data.refineLevel);
- newMesh.view("REFINED MESH");
-
// Check mesh dimension
CPPUNIT_ASSERT_EQUAL(data.cellDim, newMesh.dimension());
More information about the CIG-COMMITS
mailing list