[cig-commits] r18594 - short/3D/PyLith/trunk/libsrc/pylith/topology
knepley at geodynamics.org
knepley at geodynamics.org
Sun Jun 12 11:34:40 PDT 2011
Author: knepley
Date: 2011-06-12 11:34:40 -0700 (Sun, 12 Jun 2011)
New Revision: 18594
Modified:
short/3D/PyLith/trunk/libsrc/pylith/topology/CellRefinerHex8.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/RefineVol8Face4Edges2.cc
Log:
For now, ignore refinement points inserted around fault boundaries
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/CellRefinerHex8.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/CellRefinerHex8.cc 2011-06-12 04:13:56 UTC (rev 18593)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/CellRefinerHex8.cc 2011-06-12 18:34:40 UTC (rev 18594)
@@ -81,6 +81,28 @@
int numVolumes = 0;
const VolumeType* volumes;
+#if 0
+ // Debugging for the solvertest benchmark
+ bool hasA = false, hasB = false;
+ bool hasC = false, hasD = false;
+ for(int i = 0; i < coneSize; ++i) {
+ if (cone[i] == 5708) {hasA = true;}
+ else if (cone[i] == 8187) {hasB = true;}
+ else if (cone[i] == 4262) {hasC = true;}
+ else if (cone[i] == 8246) {hasD = true;}
+ }
+ if ((hasA && hasB) || (hasC && hasD)) {
+ const char *rankname = hasC && hasD ? "[1]" : "[0]";
+ const char *typeName = _cellType(cell) == HEXAHEDRON ? "hex" : "quadLag";
+
+ std::cout << rankname<<"Refined cell " << cell << " of type " << typeName << " with cone:" << std::endl;
+ for(int i = 0; i < coneSize; ++i) {
+ std::cout << " " << cone[i];
+ }
+ std::cout << std::endl;
+ }
+#endif
+
switch (_cellType(cell)) {
case HEXAHEDRON:
_edges_HEXAHEDRON(&edges, &numEdges, cone, coneSize);
@@ -98,8 +120,21 @@
for(int iEdge=0; iEdge < numEdges; ++iEdge) {
if (_edgeToVertex.find(edges[iEdge]) == _edgeToVertex.end()) {
+#if 0
+ // Debugging for the solvertest benchmark
+ point_type min = std::min(edges[iEdge].first, edges[iEdge].second);
+ point_type max = std::max(edges[iEdge].first, edges[iEdge].second);
+
+ if (min == 5708 && max == 8187) {
+ std::cout << "Refined cell " << cell << " with cone:" << std::endl;
+ for(int i = 0; i < coneSize; ++i) {
+ std::cout << " " << cone[i];
+ }
+ std::cout << std::endl;
+ std::cout << " edge " << edges[iEdge] << " new vertex " << *curNewVertex << std::endl;
+ }
+#endif
// if vertex does not exist
- //std::cout << "Edge: " << edges[iEdge] << ", new vertex: " << *curNewVertex << std::endl;
_edgeToVertex[edges[iEdge]] = *curNewVertex;
++(*curNewVertex);
} // if
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/RefineVol8Face4Edges2.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/RefineVol8Face4Edges2.cc 2011-06-12 04:13:56 UTC (rev 18593)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/RefineVol8Face4Edges2.cc 2011-06-12 18:34:40 UTC (rev 18594)
@@ -313,7 +313,7 @@
// Notice that points are converted to the new numbering with refined cells
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;
+ std::map<EdgeType, std::vector<int> > bndryEdgeToRank; // Maps an edge to a set of ranks who share it
const int localNormalOffset = orderNewMesh.verticesNormal().min() - orderOldMesh.verticesNormal().min();
const int localCensoredOffset = orderNewMesh.verticesCensored().min() - orderOldMesh.verticesCensored().min();
@@ -461,20 +461,27 @@
// Copy across overlap
typedef ALE::Pair<int, point_type> overlap_point_type;
+ // This maps (remote rank, remote min point) --> (remote max point, new remote vertex)
Obj<ALE::Section<overlap_point_type, EdgeType> > overlapVertices = new ALE::Section<overlap_point_type, EdgeType>(oldMesh->comm());
+ // This maps (remote rank, remote min point) --> (other remote points..., new remote vertex)
Obj<ALE::Section<overlap_point_type, FaceType> > overlapFaceVertices = new ALE::Section<overlap_point_type, FaceType>(oldMesh->comm());
ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newVerticesSection, overlapVertices);
ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newFaceVerticesSection, overlapFaceVertices);
// Merge by translating edge to local points, finding edge in _edgeToVertex, and adding (local new vetex, remote new vertex) to overlap
+ // Loop over all shared edges
for(std::map<EdgeType, std::vector<int> >::const_iterator e_iter = bndryEdgeToRank.begin(); e_iter != bndryEdgeToRank.end(); ++e_iter) {
+ // Point added on this edge by refinement
const point_type localPoint = _edgeToVertex[e_iter->first];
-
+
+ // Loop over all ranks which share this edge
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;
+ // Offset for left edge point in the new mesh (check whether edge is between Lagrange vertices)
const int localOffsetL = orderOldMesh.verticesNormal().hasPoint(e_iter->first.first) ? localNormalOffset : localCensoredOffset;
-
+
+ // Find the rank which owns the left edge point, and return the remotePoint in the new mesh
const Obj<mesh_type::send_overlap_type::traits::supportSequence>& leftRanks = newSendOverlap->support(e_iter->first.first+localOffsetL);
for(mesh_type::send_overlap_type::traits::supportSequence::iterator lr_iter = leftRanks->begin(); lr_iter != leftRanks->end(); ++lr_iter) {
if (rank == *lr_iter) {
@@ -482,7 +489,10 @@
break;
} // if
} // for
+ // Offset for right edge point in the new mesh (check whether edge is between Lagrange vertices)
const int localOffsetR = orderOldMesh.verticesNormal().hasPoint(e_iter->first.second) ? localNormalOffset : localCensoredOffset;
+
+ // Find the rank which owns the right edge point, and return the remotePoint in the new mesh
const Obj<mesh_type::send_overlap_type::traits::supportSequence>& rightRanks = newSendOverlap->support(e_iter->first.second+localOffsetR);
for(mesh_type::send_overlap_type::traits::supportSequence::iterator rr_iter = rightRanks->begin(); rr_iter != rightRanks->end(); ++rr_iter) {
if (rank == *rr_iter) {
@@ -496,15 +506,42 @@
const EdgeType *remoteVals = overlapVertices->restrictPoint(overlap_point_type(rank, remoteMin));
point_type remotePoint = -1;
+ // Loop over all shared edges from rank with endpoint remoteMin
for(int d = 0; d < remoteSize; ++d) {
if (remoteVals[d].first == remoteMax) {
remotePoint = remoteVals[d].second;
break;
} // if
} // for
+#if 0
+ // Debugging for fault edge problem
+ if (remotePoint < 0) {
+ std::cout << "["<<oldMesh->commRank()<<"]Failed to find remote partner for edge " << e_iter->first << " from rank " << rank << std::endl;
+ std::cout << "["<<oldMesh->commRank()<<"] remote edge ("<<remoteLeft<<","<<remoteRight<<") had remoteSize " << remoteSize << std::endl;
+ for(int d = 0; d < remoteSize; ++d) {
+ std::cout << "["<<oldMesh->commRank()<<"] remote val " << remoteVals[d] << std::endl;
+ }
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& leftRanks2 = oldSendOverlap->support(e_iter->first.first);
+ for(mesh_type::send_overlap_type::traits::supportSequence::iterator lr_iter = leftRanks2->begin(); lr_iter != leftRanks2->end(); ++lr_iter) {
+ if (rank == *lr_iter) {
+ std::cout << "["<<oldMesh->commRank()<<"] left match: old vertex " << lr_iter.color() << std::endl;
+ break;
+ }
+ }
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& rightRanks2 = oldSendOverlap->support(e_iter->first.second);
+ for(mesh_type::send_overlap_type::traits::supportSequence::iterator rr_iter = rightRanks2->begin(); rr_iter != rightRanks2->end(); ++rr_iter) {
+ if (rank == *rr_iter) {
+ std::cout << "["<<oldMesh->commRank()<<"] right match: old vertex " << rr_iter.color() << std::endl;
+ break;
+ }
+ }
+ }
assert(remotePoint >= 0);
- newSendOverlap->addArrow(localPoint, rank, remotePoint);
- newRecvOverlap->addArrow(rank, localPoint, remotePoint);
+#endif
+ if (remotePoint >= 0) {
+ newSendOverlap->addArrow(localPoint, rank, remotePoint);
+ newRecvOverlap->addArrow(rank, localPoint, remotePoint);
+ }
} // for
} // for
// Merge by translating face to local points, finding face in _faceToVertex, and adding (local new vetex, remote new vertex) to overlap
@@ -551,9 +588,14 @@
} // if
} // for
assert(localPoint >= orderNewMesh.verticesNormal().min() && localPoint < orderNewMesh.verticesNormal().max());
+#if 0
+ // Debugging for fault edge problem
assert(remotePoint >= 0);
- newSendOverlap->addArrow(localPoint, rank, remotePoint);
- newRecvOverlap->addArrow(rank, localPoint, remotePoint);
+#endif
+ if (remotePoint >= 0) {
+ newSendOverlap->addArrow(localPoint, rank, remotePoint);
+ newRecvOverlap->addArrow(rank, localPoint, remotePoint);
+ }
} // for
} // for
More information about the CIG-COMMITS
mailing list