From 147d0144a065c892204f6017a0fcc0d11ba57fa8 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 11 Apr 2011 15:13:49 +0000
Subject: [PATCH] Fixed bug in 3d refinement in parallel computations.
---
AMDiS/src/Debug.h | 2 +-
AMDiS/src/Global.h | 6 +++
AMDiS/src/RefinementManager3d.cc | 66 +++++++++++++++++++--------
AMDiS/src/parallel/MeshDistributor.cc | 23 ++++++----
AMDiS/src/parallel/MeshPartitioner.cc | 2 +-
AMDiS/src/parallel/ParallelDebug.cc | 3 ++
6 files changed, 70 insertions(+), 32 deletions(-)
diff --git a/AMDiS/src/Debug.h b/AMDiS/src/Debug.h
index 55e912f0..bdabf133 100644
--- a/AMDiS/src/Debug.h
+++ b/AMDiS/src/Debug.h
@@ -60,7 +60,7 @@ namespace AMDiS {
/** \brief
* Create a vtu file with name 'dofindex.vtu'. All nodes in the mesh are colored
- * by the global dof index.
+ * by the global DOF index.
*
* \param[in] feSpace The FE space to be used.
*/
diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h
index 4435b243..abbdbfde 100644
--- a/AMDiS/src/Global.h
+++ b/AMDiS/src/Global.h
@@ -304,6 +304,12 @@ namespace AMDiS {
/// prints a message
#define MSG Msg::print_funcname(funcName), Msg::print
+#if (DEBUG == 0)
+ #define MSG_DBG
+#else
+ #define MSG_DBG Msg::print_funcname(funcName), Msg::print
+#endif
+
/// prints a message, if min(Msg::msgInfo, info) >= noinfo
#define INFO(info,noinfo) \
if (Msg::getMsgInfo() && (std::min(Msg::getMsgInfo(), (info)) >= (noinfo))) MSG
diff --git a/AMDiS/src/RefinementManager3d.cc b/AMDiS/src/RefinementManager3d.cc
index d3cd7035..80927697 100644
--- a/AMDiS/src/RefinementManager3d.cc
+++ b/AMDiS/src/RefinementManager3d.cc
@@ -453,19 +453,19 @@ namespace AMDiS {
break;
TEST_EXIT_DBG(edgeDof0 < vertices)
- ("dof %d not found on element %d with nodes (%d %d %d %d)\n",
+ ("DOF %d not found on element %d with nodes (%d %d %d %d)\n",
edge[0][0], neigh->getIndex(), neigh->getDof(0, 0),
neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
TEST_EXIT_DBG(edgeDof1 < vertices)
- ("dof %d not found on element %d with nodes (%d %d %d %d)\n",
+ ("DOF %d not found on element %d with nodes (%d %d %d %d)\n",
edge[1][0], neigh->getIndex(), neigh->getDof(0, 0),
neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
}
}
TEST_EXIT_DBG(edgeDof0 < vertices && edgeDof1 < vertices)
- ("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n",
+ ("DOF %d or DOF %d not found on element %d with nodes (%d %d %d %d)\n",
edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDof(0, 0),
neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
@@ -621,8 +621,6 @@ namespace AMDiS {
if (el->getMark() <= 0)
return elInfo;
- // MSG("Refine: %d\n", elInfo->getElement()->getIndex());
-
int bound = false;
DegreeOfFreedom *edge[2];
@@ -681,24 +679,52 @@ namespace AMDiS {
DofEdge edge2 = elInfo2->getElement()->getEdge(i);
if (edge2.first == *(edge[0]) &&
edge2.second == *(edge[1])) {
- refineList.setElType(n_neigh, elInfo2->getType());
- refineList.setElement(n_neigh, elInfo2->getElement());
- n_neigh++;
- foundEdge = true;
- TraverseStack *tmpStack = stack;
- stack = &stack2;
-
- if (getRefinePatch(&elInfo2, edge, 0, refineList, &n_neigh)) {
- getRefinePatch(&elInfo2, edge, 1, refineList, &n_neigh);
- bound = true;
- }
+ // We found the edge on the other element side on leaf level. Two
+ // cases may occure: either it is already a refinement edge, i.e,
+ // it has local edge number 0, or it is not a refinement edge. In
+ // the first case, we are fine and this leaf element may be set
+ // to all elements on the refinement edge. Otherwise, the element
+ // must be refined at least once to get a refinement edge.
+
+ if (i == 0) {
+ // Edge is refinement edge, so add it to refine list.
+
+ refineList.setElType(n_neigh, elInfo2->getType());
+ refineList.setElement(n_neigh, elInfo2->getElement());
+ n_neigh++;
+
+ foundEdge = true;
+ TraverseStack *tmpStack = stack;
+ stack = &stack2;
+
+ if (getRefinePatch(&elInfo2, edge, 0, refineList, &n_neigh)) {
+ getRefinePatch(&elInfo2, edge, 1, refineList, &n_neigh);
+ bound = true;
+ }
+
+ stack = tmpStack;
+ break;
+ } else {
+ // Edge i not refinement edge, so refine the element further.
- stack = tmpStack;
- break;
+ Element *el2 = elInfo->getElement();
+ el2->setMark(std::max(el2->getMark(), 1));
+
+ TraverseStack *tmpStack = stack;
+ stack = &stack2;
+
+ elInfo2 = refineFunction(elInfo2);
+
+ stack = tmpStack;
+ break;
+ }
}
}
+ if (foundEdge)
+ break;
+
elInfo2 = stack2.traverseNext(elInfo2);
}
@@ -778,8 +804,8 @@ namespace AMDiS {
void FixRefinementPatch::getOtherEl(TraverseStack *stack,
- Element **otherEl,
- int &otherEdge)
+ Element **otherEl,
+ int &otherEdge)
{
FUNCNAME("FixRefinementPatch::getOtherEl()");
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index db0137c0..5a38b94d 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -170,7 +170,7 @@ namespace AMDiS {
}
#if (DEBUG != 0)
- ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat");
+ ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat");
#endif
initialized = true;
@@ -558,6 +558,7 @@ namespace AMDiS {
FixRefinementPatch::connectedEdges.clear();
+ bool valid3dMesh = true;
for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
vector<ElementObjectData>& edgeEls = elObjects.getElements(*it);
@@ -603,8 +604,10 @@ namespace AMDiS {
} while (found);
if (!el1.empty()) {
- MSG("Unvalid 3D mesh with %d (%d - %d) elements on this edge!\n",
- edgeEls.size(), el0.size(), el1.size());
+ valid3dMesh = false;
+
+ MSG_DBG("Unvalid 3D mesh with %d (%d - %d) elements on this edge!\n",
+ edgeEls.size(), el0.size(), el1.size());
for (std::set<int>::iterator elIt0 = el0.begin(); elIt0 != el0.end(); ++elIt0) {
for (std::set<int>::iterator elIt1 = el1.begin(); elIt1 != el1.end(); ++elIt1) {
@@ -626,10 +629,10 @@ namespace AMDiS {
mesh->getDofIndexCoords(dofEdge0.first, feSpace, c0);
mesh->getDofIndexCoords(dofEdge0.second, feSpace, c1);
- MSG("FOUND EDGE: %d %d with coords %f %f %f %f %f %f\n",
- dofEdge0.first, dofEdge0.second,
- c0[0], c0[1], c0[2],
- c1[0], c1[1], c1[2]);
+ MSG_DBG("FOUND EDGE: %d %d with coords %f %f %f %f %f %f\n",
+ dofEdge0.first, dofEdge0.second,
+ c0[0], c0[1], c0[2],
+ c1[0], c1[1], c1[2]);
TEST_EXIT_DBG(dofEdge0 == dofEdge1)("Should noth happen!\n");
@@ -639,6 +642,8 @@ namespace AMDiS {
}
}
}
+
+ MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
}
@@ -768,9 +773,7 @@ namespace AMDiS {
// === Check the boundaries and adapt mesh if necessary. ===
-#if (DEBUG != 0)
- MSG("Run checkAndAdaptBoundary ...\n");
-#endif
+ MSG_DBG("Run checkAndAdaptBoundary ...\n");
meshChanged |= checkAndAdaptBoundary(allBound);
diff --git a/AMDiS/src/parallel/MeshPartitioner.cc b/AMDiS/src/parallel/MeshPartitioner.cc
index 5d29a67b..16b81847 100644
--- a/AMDiS/src/parallel/MeshPartitioner.cc
+++ b/AMDiS/src/parallel/MeshPartitioner.cc
@@ -47,7 +47,7 @@ namespace AMDiS {
vertexElements[el->getEdge(0)].insert(elIndex);
} else {
int elInRank = std::min(elIndex / elPerRank, mpiSize - 1);
-
+
elementInRank[elIndex] = (elInRank == mpiRank);
partitionMap[elIndex] = elInRank;
}
diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc
index a98db4c0..90fd2078 100644
--- a/AMDiS/src/parallel/ParallelDebug.cc
+++ b/AMDiS/src/parallel/ParallelDebug.cc
@@ -805,6 +805,9 @@ namespace AMDiS {
// === Write to all elements in ranks mesh the included dofs. ===
+ file << "\n\n";
+ file << "# First line containes number of elements in mesh, second line contain the number of DOFs per element.\n";
+ file << "# Than, each entry contains of two lines. The first is the element index, the second line is a list with the local DOF indices of this element.\n";
file << elDofMap.size() << "\n";
file << basisFcts->getNumber() << "\n";
for (ElDofMap::iterator it = elDofMap.begin(); it != elDofMap.end(); ++it) {
--
GitLab