From 9bb9b920ccfbb16b8eab683c25e8dd8416793974 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 30 Jul 2010 07:02:43 +0000
Subject: [PATCH] Fixed some problems in 3D refinement with reversed mesh
traversel stratagey.
---
AMDiS/src/CoarseningManager3d.cc | 2 +-
AMDiS/src/Debug.cc | 46 +++-
AMDiS/src/Debug.h | 2 +
AMDiS/src/MacroReader.cc | 6 +-
AMDiS/src/MeshStructure.cc | 27 +++
AMDiS/src/MeshStructure.h | 24 +--
AMDiS/src/RefinementManager.cc | 2 +-
AMDiS/src/RefinementManager3d.cc | 290 ++++++++++++++------------
AMDiS/src/RefinementManager3d.h | 20 +-
AMDiS/src/Tetrahedron.cc | 8 +-
AMDiS/src/Tetrahedron.h | 2 +-
AMDiS/src/Traverse.cc | 60 ++++--
AMDiS/src/Traverse.h | 7 +-
AMDiS/src/parallel/MeshDistributor.cc | 94 ++++++++-
AMDiS/src/parallel/ParallelDebug.cc | 21 +-
15 files changed, 394 insertions(+), 217 deletions(-)
diff --git a/AMDiS/src/CoarseningManager3d.cc b/AMDiS/src/CoarseningManager3d.cc
index 749a6a92..e05f1a7b 100644
--- a/AMDiS/src/CoarseningManager3d.cc
+++ b/AMDiS/src/CoarseningManager3d.cc
@@ -257,7 +257,7 @@ namespace AMDiS {
edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDOF(0,0),
neigh->getDOF(1,0), neigh->getDOF(2,0), neigh->getDOF(3,0));
}
- edge_no = Tetrahedron::edgeOfDOFs[j][k];
+ edge_no = Tetrahedron::edgeOfDofs[j][k];
coarsenList->setCoarsePatch(*n_neigh, edge_no == 0);
/****************************************************************************/
diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc
index 0765ccf3..2d647673 100644
--- a/AMDiS/src/Debug.cc
+++ b/AMDiS/src/Debug.cc
@@ -231,21 +231,22 @@ namespace AMDiS {
void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof)
{
+ FUNCNAME("debug::printInfoByDof()");
+
Element *el = getDofIndexElement(feSpace, dof);
Element *parEl = getLevel0ParentElement(feSpace->getMesh(), el);
- std::cout << "[DBG] DOF-INFO: dof = " << dof
- << " elidx = " << el->getIndex()
- << " pelidx = " << parEl->getIndex() << std::endl;
+ MSG("[DBG] DOF-INFO: dof = %d elidx = %d pelidx = %d\n",
+ dof, el->getIndex(), parEl->getIndex());
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
Mesh::CALL_EVERY_EL_PREORDER);
while (elInfo) {
- if (elInfo->getElement()->getIndex() == parEl->getIndex())
- std::cout << "[DBG] EL INFO TO " << parEl->getIndex() << ": "
- << " type = " << elInfo->getType() << std::endl;
-
+ if (elInfo->getElement()->getIndex() == parEl->getIndex()) {
+ MSG("[DBG] EL INFO TO %d: type = %d\n", parEl->getIndex(), elInfo->getType());
+ }
+
elInfo = stack.traverseNext(elInfo);
}
}
@@ -470,6 +471,37 @@ namespace AMDiS {
}
+ int getLocalNeighbourIndex(Mesh *mesh, int elIndex, int neighIndex)
+ {
+ FUNCNAME("debug::getLocalNeighbourIndex()");
+
+ TraverseStack stack;
+ ElInfo *elInfo =
+ stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH);
+ while (elInfo) {
+ if (elInfo->getElement()->getIndex() == elIndex) {
+ int returnValue = -1;
+ for (int i = 0; i <= mesh->getDim(); i++) {
+ if (elInfo->getNeighbour(i)) {
+ MSG("NEIGH %d of El-Idx %d: %d\n", i, elIndex, elInfo->getNeighbour(i)->getIndex());
+ } else {
+ MSG("NEIGH %d of El-Idx %d: %d\n", i, elIndex, -1);
+ }
+
+ if (elInfo->getNeighbour(i) &&
+ elInfo->getNeighbour(i)->getIndex() == neighIndex)
+ returnValue = i;
+ }
+
+ return returnValue;
+ }
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ return -2;
+ }
+
+
void createSortedDofs(Mesh *mesh, ElementIdxToDofs &elMap)
{
FUNCNAME("debug::dbgCreateElementMap()");
diff --git a/AMDiS/src/Debug.h b/AMDiS/src/Debug.h
index e65ae42a..8432077f 100644
--- a/AMDiS/src/Debug.h
+++ b/AMDiS/src/Debug.h
@@ -135,6 +135,8 @@ namespace AMDiS {
void printElementHierarchie(Mesh *mesh, int elIndex);
+ int getLocalNeighbourIndex(Mesh *mesh, int elIndex, int neighIndex);
+
/** \brief
* Traverse a mesh and store for each element all its vertex DOFs in local sorted
* order (by values).
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index b7cc591e..4f5fe6c4 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -1346,7 +1346,7 @@ namespace AMDiS {
if (!mesh->getNumberOfDOFs(EDGE) && nei->getIndex() < mel_index)
return false;
- edge_no = Tetrahedron::edgeOfDOFs[j][k];
+ edge_no = Tetrahedron::edgeOfDofs[j][k];
TEST_EXIT(*n_neigh < max_no_list_el)
("too many neigbours for local list\n");
@@ -1428,7 +1428,7 @@ namespace AMDiS {
if (nei->getIndex() < mel_index)
return false;
- edge_no = Tetrahedron::edgeOfDOFs[j][k];
+ edge_no = Tetrahedron::edgeOfDofs[j][k];
TEST_EXIT(*n_neigh < max_no_list_el)("too many neigbours for local list\n");
@@ -1440,7 +1440,7 @@ namespace AMDiS {
("dof %d on element %d is already set, but not on element %d\n",
node + edge_no, nei->getIndex(), mel_index);
- nei->element->setDOF(node+edge_no,edge_dof);
+ nei->element->setDOF(node+edge_no, edge_dof);
}
if (next_el[edge_no][0] != opp_v) {
diff --git a/AMDiS/src/MeshStructure.cc b/AMDiS/src/MeshStructure.cc
index 3d91553e..859fb2d4 100644
--- a/AMDiS/src/MeshStructure.cc
+++ b/AMDiS/src/MeshStructure.cc
@@ -308,6 +308,33 @@ namespace AMDiS {
}
+ void MeshStructure::print(bool resetCode)
+ {
+ FUNCNAME("MeshStructure::print()");
+
+ std::stringstream oss;
+
+ if (empty()) {
+ oss << "-" << std::endl;
+ } else {
+ if (resetCode)
+ reset();
+
+ bool cont = true;
+ while (cont) {
+ if (isLeafElement())
+ oss << "0";
+ else
+ oss << "1";
+
+ cont = nextElement();
+ }
+ }
+
+ MSG("Mesh structure code: %s\n", oss.str().c_str());
+ }
+
+
bool MeshStructure::compare(MeshStructure &other)
{
return (other.getCode() == code);
diff --git a/AMDiS/src/MeshStructure.h b/AMDiS/src/MeshStructure.h
index 9348c7bc..b7ee9812 100644
--- a/AMDiS/src/MeshStructure.h
+++ b/AMDiS/src/MeshStructure.h
@@ -114,29 +114,7 @@ namespace AMDiS {
int macroElIndex = -1);
/// Prints the mesh structure code.
- void print()
- {
- FUNCNAME("MeshStructure::print()");
-
- std::stringstream oss;
-
- if (empty()) {
- oss << "-" << std::endl;
- } else {
- reset();
- bool cont = true;
- while (cont) {
- if (isLeafElement())
- oss << "0";
- else
- oss << "1";
-
- cont = nextElement();
- }
- }
-
- MSG("Mesh structure code: %s\n", oss.str().c_str());
- }
+ void print(bool resetCode = true);
/// Returns the mesh structure code.
inline const std::vector<unsigned long int>& getCode()
diff --git a/AMDiS/src/RefinementManager.cc b/AMDiS/src/RefinementManager.cc
index 9386be79..133151ad 100644
--- a/AMDiS/src/RefinementManager.cc
+++ b/AMDiS/src/RefinementManager.cc
@@ -53,7 +53,7 @@ namespace AMDiS {
stack->traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
while (elInfo) {
- if (elInfo->getElement()->getMark() > 0) {
+ if (elInfo->getElement()->getMark() > 0) {
doMoreRecursiveRefine =
doMoreRecursiveRefine || (elInfo->getElement()->getMark() > 1);
elInfo = refineFunction(elInfo);
diff --git a/AMDiS/src/RefinementManager3d.cc b/AMDiS/src/RefinementManager3d.cc
index 40069bd6..40c9e0cd 100644
--- a/AMDiS/src/RefinementManager3d.cc
+++ b/AMDiS/src/RefinementManager3d.cc
@@ -146,13 +146,13 @@ namespace AMDiS {
int el_type = ref_list->getType(index);
int n_type = 0;
int dir, adjc, i, j, i_neigh, j_neigh;
- int node0, node1, opp_v = 0;
+ int node0, node1, oppVertex = 0;
for (dir = 0; dir < 2; dir++) {
neigh = ref_list->getNeighbourElement(index, dir);
if (neigh) {
n_type = ref_list->getType(ref_list->getNeighbourNr(index, dir));
- opp_v = ref_list->getOppVertex(index, dir);
+ oppVertex = ref_list->getOppVertex(index, dir);
}
if (!neigh || neigh->isLeaf()) {
@@ -199,7 +199,7 @@ namespace AMDiS {
j = Tetrahedron::adjacentChild[adjc][i];
i_neigh = Tetrahedron::nChildFace[el_type][i][dir];
- j_neigh = Tetrahedron::nChildFace[n_type][j][opp_v - 2];
+ j_neigh = Tetrahedron::nChildFace[n_type][j][oppVertex - 2];
/****************************************************************************/
/* adjust dof pointer in the edge in the common face of el and neigh and */
@@ -208,7 +208,7 @@ namespace AMDiS {
if (mesh->getNumberOfDOFs(EDGE)) {
node0 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][i][dir];
- node1 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[n_type][j][opp_v - 2];
+ node1 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[n_type][j][oppVertex - 2];
TEST_EXIT_DBG(neigh->getChild(j)->getDOF(node1))
("no dof on neighbour %d at node %d\n",
@@ -267,9 +267,8 @@ namespace AMDiS {
edge[i] = const_cast<int*>(el_info->getElement()->getDOF(i));
if (getRefinePatch(&elinfo, edge, 0, refList, &n_neigh)) {
- /****************************************************************************/
- /* domain's boundary was reached while looping around the refinement edge */
- /****************************************************************************/
+ // Domain's boundary was reached while looping around the refinement edge.
+
getRefinePatch(&elinfo, edge, 1, refList, &n_neigh);
}
@@ -383,173 +382,209 @@ namespace AMDiS {
int RefinementManager3d::getRefinePatch(ElInfo **el_info,
DegreeOfFreedom *edge[2],
- int dir,
+ int direction,
RCNeighbourList* refineList,
int *n_neigh)
{
FUNCNAME("RefinementManager3d::getRefinePatch()");
- int i, j, k;
-
+ int localNeighbour = 3 - direction;
Tetrahedron *el =
dynamic_cast<Tetrahedron*>(const_cast<Element*>((*el_info)->getElement()));
- if ((*el_info)->getNeighbour(3 - dir) == NULL)
+ if ((*el_info)->getNeighbour(localNeighbour) == NULL)
return 1;
- int opp_v = (*el_info)->getOppVertex(3-dir);
- ElInfo *neigh_info = stack->traverseNeighbour3d((*el_info), 3 - dir);
- int neigh_el_type = neigh_info->getType();
+ int oppVertex = (*el_info)->getOppVertex(localNeighbour);
+ int testIndex = (*el_info)->getNeighbour(localNeighbour)->getIndex();
+ ElInfo *neighInfo = stack->traverseNeighbour3d((*el_info), localNeighbour);
+ int neighElType = neighInfo->getType();
+
+ TEST_EXIT_DBG(neighInfo->getElement()->getIndex() == testIndex)
+ ("Should not happen!\n");
Tetrahedron *neigh =
- dynamic_cast<Tetrahedron*>(const_cast<Element*>(neigh_info->getElement()));
-
+ dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
int vertices = mesh->getGeo(VERTEX);
while (neigh != el) {
- for (j = 0; j < vertices; j++)
- if (neigh->getDOF(j) == edge[0])
+ // === Determine the common edge of the element and its neighbour. ===
+
+ int edgeDof0, edgeDof1;
+ for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
+ if (neigh->getDOF(edgeDof0) == edge[0])
break;
- for (k = 0; k < vertices; k++)
- if (neigh->getDOF(k) == edge[1])
+ for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
+ if (neigh->getDOF(edgeDof1) == edge[1])
break;
-
- if (j > 3 || k > 3) {
- for (j = 0; j < vertices; j++)
- if (mesh->associated(neigh->getDOF(j, 0), edge[0][0]))
+ if (edgeDof0 > 3 || edgeDof1 > 3) {
+ for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
+ if (mesh->associated(neigh->getDOF(edgeDof0, 0), edge[0][0]))
break;
- for (k = 0; k < vertices; k++)
- if (mesh->associated(neigh->getDOF(k, 0), edge[1][0]))
+ for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
+ if (mesh->associated(neigh->getDOF(edgeDof1, 0), edge[1][0]))
break;
- if (j > 3 || k > 3) {
- for (j = 0; j < vertices; j++)
- if (mesh->indirectlyAssociated(neigh->getDOF(j, 0), edge[0][0]))
+ if (edgeDof0 > 3 || edgeDof1 > 3) {
+ for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
+ if (mesh->indirectlyAssociated(neigh->getDOF(edgeDof0, 0), edge[0][0]))
break;
- for (k = 0; k < vertices; k++)
- if (mesh->indirectlyAssociated(neigh->getDOF(k, 0), edge[1][0]))
+ for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
+ if (mesh->indirectlyAssociated(neigh->getDOF(edgeDof1, 0), edge[1][0]))
break;
- TEST_EXIT_DBG(j < vertices && k < vertices)
- ("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));
+ TEST_EXIT_DBG(edgeDof0 < vertices && edgeDof1 < vertices)
+ ("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));
}
}
- TEST_EXIT_DBG(j < mesh->getGeo(VERTEX) &&
- k < mesh->getGeo(VERTEX))
+ TEST_EXIT_DBG(edgeDof0 < vertices && edgeDof1 < vertices)
("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));
+ edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDOF(0, 0),
+ neigh->getDOF(1, 0), neigh->getDOF(2, 0), neigh->getDOF(3, 0));
+
+ int edgeNo = Tetrahedron::edgeOfDofs[edgeDof0][edgeDof1];
+
+ if (edgeNo) {
+ // Only 0 can be a compatible commen refinement edge. Thus, neigh has not
+ // a compatible refinement edge. Refine it first.
- int edge_no = Tetrahedron::edgeOfDOFs[j][k];
- if (edge_no) {
- /****************************************************************************/
- /* neigh has not a compatible refinement edge */
- /****************************************************************************/
neigh->setMark(max(neigh->getMark(), 1));
- neigh_info = refineFunction(neigh_info);
-
- /****************************************************************************/
- /* now, go to a child at the edge and determine the opposite vertex for */
- /* this child; continue the looping around the edge with this element */
- /****************************************************************************/
+
+ neighInfo = refineFunction(neighInfo);
+
+ // === Now, go to a child at the edge and determine the opposite vertex ===
+ // === for this child; continue the looping around the edge with this ===
+ // === element. ===
- neigh_info = stack->traverseNext(neigh_info);
- neigh_el_type = neigh_info->getType();
+ neighInfo = stack->traverseNext(neighInfo);
+ neighElType = neighInfo->getType();
+ bool reverseMode = stack->getTraverseFlag().isSet(Mesh::CALL_REVERSE_MODE);
- switch (edge_no) {
+ switch (edgeNo) {
case 1:
- opp_v = opp_v == 1 ? 3 : 2;
+ if (reverseMode) {
+ neighInfo = stack->traverseNext(neighInfo);
+ neighElType = neighInfo->getType();
+ }
+
+ oppVertex = (oppVertex == 1 ? 3 : 2);
break;
case 2:
- opp_v = opp_v == 2 ? 1 : 3;
+ if (reverseMode) {
+ neighInfo = stack->traverseNext(neighInfo);
+ neighElType = neighInfo->getType();
+ }
+
+ oppVertex = (oppVertex == 2 ? 1 : 3);
break;
case 3:
- neigh_info = stack->traverseNext(neigh_info);
- neigh_el_type = neigh_info->getType();
- if (neigh_el_type != 1)
- opp_v = opp_v == 0 ? 3 : 2;
+ if (!reverseMode) {
+ neighInfo = stack->traverseNext(neighInfo);
+ neighElType = neighInfo->getType();
+ }
+
+ if (neighElType != 1)
+ oppVertex = (oppVertex == 0 ? 3 : 2);
else
- opp_v = opp_v == 0 ? 3 : 1;
+ oppVertex = (oppVertex == 0 ? 3 : 1);
break;
case 4:
- neigh_info = stack->traverseNext(neigh_info);
- neigh_el_type = neigh_info->getType();
- if (neigh_el_type != 1)
- opp_v = opp_v == 0 ? 3 : 1;
+ if (!reverseMode) {
+ neighInfo = stack->traverseNext(neighInfo);
+ neighElType = neighInfo->getType();
+ }
+
+ if (neighElType != 1)
+ oppVertex = (oppVertex == 0 ? 3 : 1);
else
- opp_v = opp_v == 0 ? 3 : 2;
+ oppVertex = (oppVertex == 0 ? 3 : 2);
break;
case 5:
- if (neigh_el_type != 1) {
- neigh_info = stack->traverseNext(neigh_info);
- neigh_el_type = neigh_info->getType();
+ if (neighElType != 1) {
+ if (!reverseMode) {
+ neighInfo = stack->traverseNext(neighInfo);
+ neighElType = neighInfo->getType();
+ }
}
- opp_v = 3;
+ oppVertex = 3;
break;
}
+
neigh =
- dynamic_cast<Tetrahedron*>(const_cast<Element*>(neigh_info->getElement()));
+ dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
} else {
- /****************************************************************************/
- /* neigh is compatible devisible; put neigh to the list of patch elements; */
- /* go to next neighbour */
- /****************************************************************************/
+ // Neigh is compatible devisible. Put neigh to the list of patch elements
+ // and go to next neighbour.
+
TEST_EXIT_DBG(*n_neigh < mesh->getMaxEdgeNeigh())
("too many neighbours %d in refine patch\n", *n_neigh);
-
-
+
refineList->setElement(*n_neigh, neigh);
- refineList->setElType(*n_neigh, neigh_el_type);
-
- /****************************************************************************/
- /* we have to go back to the starting element via opp_v values */
- /* correct information is produce by get_neigh_on_patch() */
- /****************************************************************************/
- refineList->setOppVertex(*n_neigh, 0, opp_v);
+ refineList->setElType(*n_neigh, neighElType);
- ++*n_neigh;
+ // We have to go back to the starting element via oppVertex values.
- if (opp_v != 3)
- i = 3;
- else
- i = 2;
+ refineList->setOppVertex(*n_neigh, 0, oppVertex);
- if (neigh_info->getNeighbour(i)) {
- opp_v = neigh_info->getOppVertex(i);
- neigh_info = stack->traverseNeighbour3d(neigh_info, i);
- neigh_el_type = neigh_info->getType();
+ (*n_neigh)++;
+
+ int i = (oppVertex != 3 ? 3 : 2);
+
+ if (neighInfo->getNeighbour(i)) {
+ oppVertex = neighInfo->getOppVertex(i);
+ int testIndex = neighInfo->getNeighbour(i)->getIndex();
+
+ neighInfo = stack->traverseNeighbour3d(neighInfo, i);
+
+ TEST_EXIT_DBG(neighInfo->getElement()->getIndex() == testIndex)
+ ("Should not happen!\n");
+
+ neighElType = neighInfo->getType();
neigh =
- dynamic_cast<Tetrahedron*>(const_cast<Element*>(neigh_info->getElement()));
- } else
+ dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighInfo->getElement()));
+ } else {
break;
+ }
}
}
if (neigh == el) {
- (*el_info) = neigh_info;
+ (*el_info) = neighInfo;
+
return 0;
}
- /****************************************************************************/
- /* the domain's boundary is reached; loop back to the starting el */
- /****************************************************************************/
+ // The domain's boundary is reached. Loop back to the starting el.
- i = *n_neigh - 1;
- opp_v = refineList->getOppVertex(i, 0);
+ int i = *n_neigh - 1;
+ oppVertex = refineList->getOppVertex(i, 0);
do {
- TEST_EXIT_DBG(neigh_info->getNeighbour(opp_v) && i > 0)
+ TEST_EXIT_DBG(neighInfo->getNeighbour(oppVertex) && i > 0)
("While looping back domains boundary was reached or i == 0\n");
- opp_v = refineList->getOppVertex(i--, 0);
- neigh_info = stack->traverseNeighbour3d(neigh_info, opp_v);
- } while (neigh_info->getElement() != el);
+ oppVertex = refineList->getOppVertex(i--, 0);
+
+ int testIndex = neighInfo->getNeighbour(oppVertex)->getIndex();
+ neighInfo = stack->traverseNeighbour3d(neighInfo, oppVertex);
+
+ int edgeDof0, edgeDof1;
+ for (edgeDof0 = 0; edgeDof0 < vertices; edgeDof0++)
+ if (neigh->getDOF(edgeDof0) == edge[0])
+ break;
+ for (edgeDof1 = 0; edgeDof1 < vertices; edgeDof1++)
+ if (neigh->getDOF(edgeDof1) == edge[1])
+ break;
+
+ TEST_EXIT_DBG(neighInfo->getElement()->getIndex() == testIndex)
+ ("Should not happen!\n");
+ } while (neighInfo->getElement() != el);
- (*el_info) = neigh_info;
+ (*el_info) = neighInfo;
return 1;
}
@@ -565,9 +600,8 @@ namespace AMDiS {
if (el_info->getElement()->getMark() <= 0)
return el_info; /* element may not be refined */
- /****************************************************************************/
- /* get memory for a list of all elements at the refinement edge */
- /****************************************************************************/
+ // === Get memory for a list of all elements at the refinement edge. ===
+
RCNeighbourList *ref_list = new RCNeighbourList(mesh->getMaxEdgeNeigh());
ref_list->setElType(0, el_info->getType());
ref_list->setElement(0, el_info->getElement());
@@ -577,37 +611,31 @@ namespace AMDiS {
el_info->getProjection(0)->getType() == VOLUME_PROJECTION)
newCoords = true;
- /****************************************************************************/
- /* give the refinement edge the right orientation */
- /****************************************************************************/
- if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1,0)) {
- edge[0] = const_cast<int*>(el_info->getElement()->getDOF(0));
- edge[1] = const_cast<int*>(el_info->getElement()->getDOF(1));
+ // === Give the refinement edge the right orientation. ===
+
+ if (el_info->getElement()->getDOF(0, 0) < el_info->getElement()->getDOF(1, 0)) {
+ edge[0] = const_cast<DegreeOfFreedom*>(el_info->getElement()->getDOF(0));
+ edge[1] = const_cast<DegreeOfFreedom*>(el_info->getElement()->getDOF(1));
} else {
- edge[1] = const_cast<int*>(el_info->getElement()->getDOF(0));
- edge[0] = const_cast<int*>(el_info->getElement()->getDOF(1));
+ edge[1] = const_cast<DegreeOfFreedom*>(el_info->getElement()->getDOF(0));
+ edge[0] = const_cast<DegreeOfFreedom*>(el_info->getElement()->getDOF(1));
}
- /****************************************************************************/
- /* get the refinement patch */
- /****************************************************************************/
+
+ // === Traverse and refine the refinement patch. ====
+
if (getRefinePatch(&el_info, edge, 0, ref_list, &n_neigh)) {
- /****************************************************************************/
- /* domain's boundary was reached while looping around the refinement edge */
- /****************************************************************************/
+ // Domain's boundary was reached while looping around the refinement edge
getRefinePatch(&el_info, edge, 1, ref_list, &n_neigh);
bound = true;
}
- /****************************************************************************/
- /* fill neighbour information inside the patch in the refinement list */
- /****************************************************************************/
+ // fill neighbour information inside the patch in the refinement list
ref_list->getNeighOnPatch(n_neigh, bound);
- // ==========================================================================
- // === check for periodic boundary ==========================================
- // ==========================================================================
+
+ // ============ Check for periodic boundary ============
DegreeOfFreedom *next_edge[2];
DegreeOfFreedom *first_edge[2] = {edge[0], edge[1]};
@@ -671,12 +699,6 @@ namespace AMDiS {
}
}
-
- /****************************************************************************/
- /* and now refine the patch */
- /****************************************************************************/
- //refinePatch(edge, ref_list, n_neigh, bound);
-
stack->update();
delete ref_list;
diff --git a/AMDiS/src/RefinementManager3d.h b/AMDiS/src/RefinementManager3d.h
index 401bc311..875f4cef 100644
--- a/AMDiS/src/RefinementManager3d.h
+++ b/AMDiS/src/RefinementManager3d.h
@@ -47,15 +47,19 @@ namespace AMDiS {
void setNewCoords();
/** \brief
- * gets the elements around the refinement edge with vertices node[0] and
- * node[1] ; refines those elements at this edge are not compatible
- * devisible;
- * dir determines the direction of the first neighbour
- * get_refine_patch returns 1 if the domain's boundary is reached while
- * looping around the refinement edge, otherwise 0
+ * Gets the elements around the refinement edge with vertices node[0] and
+ * node[1]. Refines those elements at this edge that are not compatible
+ * devisible. The function returns 1 if the domain's boundary is reached
+ * while looping around the refinement edge, otherwise 0.
+ *
+ * \param[in] direction Determines the direction of the first neighbour
+ *
*/
- int getRefinePatch(ElInfo **el_info, DegreeOfFreedom *edge[2], int dir,
- RCNeighbourList* refineList, int *n_neigh);
+ int getRefinePatch(ElInfo **el_info,
+ DegreeOfFreedom *edge[2],
+ int direction,
+ RCNeighbourList* refineList,
+ int *n_neigh);
/// Refines all elements in the patch.
DegreeOfFreedom refinePatch(DegreeOfFreedom *edge[2], RCNeighbourList* refineList,
diff --git a/AMDiS/src/Tetrahedron.cc b/AMDiS/src/Tetrahedron.cc
index 25ce025b..3609f7e2 100644
--- a/AMDiS/src/Tetrahedron.cc
+++ b/AMDiS/src/Tetrahedron.cc
@@ -30,10 +30,10 @@ namespace AMDiS {
{1,-1},
{1,-1}};
- const unsigned char Tetrahedron::edgeOfDOFs[4][4] = {{255,0,1,2},
- {0,255,3,4},
- {1,3,255,5},
- {2,4,5,255}};
+ const unsigned char Tetrahedron::edgeOfDofs[4][4] = {{255, 0, 1, 2},
+ {0, 255, 3, 4},
+ {1, 3, 255, 5},
+ {2, 4, 5, 255}};
const int Tetrahedron::vertexOfEdge[6][2] = {{0,1},
{0,2},
diff --git a/AMDiS/src/Tetrahedron.h b/AMDiS/src/Tetrahedron.h
index 67547218..1ac2cf38 100644
--- a/AMDiS/src/Tetrahedron.h
+++ b/AMDiS/src/Tetrahedron.h
@@ -200,7 +200,7 @@ namespace AMDiS {
static const signed char childOrientation[3][2];
/// edgeOfDOFs[i][j]: gives the local index of edge with vertices i and j
- static const unsigned char edgeOfDOFs[4][4];
+ static const unsigned char edgeOfDofs[4][4];
///
static const int edgeOfFace[4][3];
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index e6bda8db..b8085a50 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -10,6 +10,7 @@
#include "FixVec.h"
#include "Parametric.h"
#include "OpenMP.h"
+#include "Debug.h"
namespace AMDiS {
@@ -73,9 +74,9 @@ namespace AMDiS {
else if (traverse_fill_flag.isSet(Mesh::CALL_MG_LEVEL))
elinfo = traverseMultiGridLevel();
else
- if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_PREORDER))
+ if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_PREORDER)) {
elinfo = traverseEveryElementPreorder();
- else if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_INORDER))
+ } else if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_INORDER))
elinfo = traverseEveryElementInorder();
else if (traverse_fill_flag.isSet(Mesh::CALL_EVERY_EL_POSTORDER))
elinfo = traverseEveryElementPostorder();
@@ -347,10 +348,13 @@ namespace AMDiS {
enlargeTraverseStack();
int fillIthChild = info_stack[stack_used];
+
info_stack[stack_used]++;
if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
fillIthChild = 1 - fillIthChild;
+
elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
+
stack_used++;
TEST_EXIT_DBG(stack_used < stack_size)
@@ -453,18 +457,17 @@ namespace AMDiS {
ElInfo* TraverseStack::traverseNeighbour3d(ElInfo* elinfo_old, int neighbour)
{
- FUNCNAME("TraverseStackLLtraverseNeighbour3d()");
+ FUNCNAME("TraverseStack::traverseNeighbour3d()");
Element *el2;
ElInfo *elinfo2;
- const DegreeOfFreedom *dof;
int stack2_used;
int sav_neighbour = neighbour;
- // father.neigh[coarse_nb[i][j]] == child[i-1].neigh[j]
- static int coarse_nb[3][3][4] = {{{-2,-2,-2,-2}, {-1,2,3,1}, {-1,3,2,0}},
- {{-2,-2,-2,-2}, {-1,2,3,1}, {-1,2,3,0}},
- {{-2,-2,-2,-2}, {-1,2,3,1}, {-1,2,3,0}}};
+ // father.neigh[coarse_nb[i][j]] == child[i - 1].neigh[j]
+ static int coarse_nb[3][3][4] = {{{-2, -2, -2, -2}, {-1, 2, 3, 1}, {-1, 3, 2, 0}},
+ {{-2, -2, -2, -2}, {-1, 2, 3, 1}, {-1, 2, 3, 0}},
+ {{-2, -2, -2, -2}, {-1, 2, 3, 1}, {-1, 2, 3, 0}}};
TEST_EXIT_DBG(stack_used > 0)("no current element\n");
@@ -479,7 +482,7 @@ namespace AMDiS {
Element *el = elinfo_stack[stack_used]->getElement();
int sav_index = el->getIndex();
- /* first, goto to leaf level, if necessary... */
+ // First, goto to leaf level, if necessary ...
if ((traverse_fill_flag & Mesh::CALL_LEAF_EL).isAnySet()) {
if (el->getChild(0) && neighbour < 2) {
if (stack_used >= stack_size - 1)
@@ -517,12 +520,12 @@ namespace AMDiS {
elinfo_stack[stack_used + 1]->getElement())
("Should not happen!\n");
-
nb = coarse_nb[typ][elIsIthChild][nb];
+
if (nb == -1)
break;
- TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
+ TEST_EXIT_DBG(nb >= 0)("Invalid coarse_nb %d!\n", nb);
}
for (int i = stack_used; i <= save_stack_used; i++) {
@@ -537,6 +540,7 @@ namespace AMDiS {
// Go to macro element neighbour.
int i = traverse_mel->getOppVertex(nb);
+
traverse_mel = traverse_mel->getNeighbour(nb);
if (traverse_mel == NULL)
return NULL;
@@ -562,7 +566,7 @@ namespace AMDiS {
elinfo2 = save_elinfo_stack[stack2_used];
el2 = elinfo2->getElement();
- if (stack_used >= stack_size-1)
+ if (stack_used >= stack_size - 1)
enlargeTraverseStack();
int i = 2 - info_stack[stack_used];
@@ -570,7 +574,7 @@ namespace AMDiS {
int fillIthChild = i;
if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
fillIthChild = 1 - fillIthChild;
- elinfo_stack[stack_used+1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
+ elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
stack_used++;
info_stack[stack_used] = 0;
nb = 0;
@@ -588,11 +592,16 @@ namespace AMDiS {
if (stack_used >= stack_size - 1)
enlargeTraverseStack();
- info_stack[stack_used] = 2 - nb;
+
+
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) {
+ int t = 2 - nb;
+ info_stack[stack_used] = (t == 2 ? 1 : 2);
+ } else {
+ info_stack[stack_used] = 2 - nb;
+ }
int fillIthChild = 1 - nb;
- if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
- fillIthChild = 1 - fillIthChild;
elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
stack_used++;
@@ -618,18 +627,23 @@ namespace AMDiS {
else if (traverse_mesh->associated(el->getDOF(1, 0), el2->getDOF(0, 0)))
i = 2 - save_info_stack[stack2_used];
else {
- ERROR_EXIT("no common refinement edge\n");
+ ERROR_EXIT("No common refinement edge! %d\n", traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE));
}
}
+ int testChild = i;
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
+ testChild = 1 - testChild;
+
if (el->getChild(0) &&
- (el->getChild(i)->getDOF(1) == el->getDOF(nb) ||
- traverse_mesh->associated(el->getChild(i)->getDOF(1, 0), el->getDOF(nb, 0))))
+ (el->getChild(testChild)->getDOF(1) == el->getDOF(nb) ||
+ traverse_mesh->associated(el->getChild(testChild)->getDOF(1, 0), el->getDOF(nb, 0))))
nb = 1;
else
nb = 2;
- info_stack[stack_used] = i + 1;
+ info_stack[stack_used] = i + 1;
+
if (stack_used >= stack_size - 1)
enlargeTraverseStack();
@@ -647,8 +661,9 @@ namespace AMDiS {
stack2_used++;
elinfo2 = save_elinfo_stack[stack2_used];
el2 = elinfo2->getElement();
+
if (save_stack_used > stack2_used) {
- dof = el2->getDOF(1);
+ const DegreeOfFreedom *dof = el2->getDOF(1);
if (dof != el->getDOF(1) && dof != el->getDOF(2) &&
!traverse_mesh->associated(dof[0], el->getDOF(1, 0)) &&
!traverse_mesh->associated(dof[0], el->getDOF(2, 0))) {
@@ -663,6 +678,7 @@ namespace AMDiS {
elinfo = elinfo_stack[stack_used];
el = elinfo->getElement();
+
break;
}
}
@@ -706,7 +722,7 @@ namespace AMDiS {
int sav_neighbour = neighbour;
// father.neigh[coarse_nb[i][j]] == child[i-1].neigh[j]
- static int coarse_nb[3][3] = {{-2,-2,-2}, {2,-1,1}, {-1,2,0}};
+ static int coarse_nb[3][3] = {{-2, -2, -2}, {2, -1, 1}, {-1, 2, 0}};
TEST_EXIT_DBG(stack_used > 0)("no current element");
diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h
index 1fc1a23f..2ab64a45 100644
--- a/AMDiS/src/Traverse.h
+++ b/AMDiS/src/Traverse.h
@@ -62,7 +62,7 @@ namespace AMDiS {
stack_used(0),
save_stack_used(0),
myThreadId(0),
- maxThreads(1)
+ maxThreads(1)
{}
/// Destructor
@@ -150,6 +150,11 @@ namespace AMDiS {
return elinfo_stack[stack_used];
}
+ inline Flag getTraverseFlag()
+ {
+ return traverse_fill_flag;
+ }
+
private:
/// Enlargement of the stack
void enlargeTraverseStack();
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 154b16c1..fff97781 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -551,6 +551,9 @@ namespace AMDiS {
std::map<int, MeshCodeVec> sendCodes;
for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
+
+ // MSG("HAVE %d BOUNDARIES WITH RANK %d!\n", it->second.size(), it->first);
+
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end(); ++boundIt) {
MeshStructure elCode;
@@ -568,6 +571,15 @@ namespace AMDiS {
bool meshFitTogether = true;
+// if (mpiRank == 0) {
+// debug::highlightElementIndexMesh(mesh, 139, "rank0_139.vtu");
+// debug::highlightElementIndexMesh(mesh, 161, "rank0_161.vtu");
+// debug::printElementHierarchie(mesh, 139);
+// } else {
+// debug::highlightElementIndexMesh(mesh, 137, "rank1_137.vtu");
+// debug::highlightElementIndexMesh(mesh, 163, "rank1_163.vtu");
+// }
+
for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
@@ -577,21 +589,56 @@ namespace AMDiS {
boundIt != it->second.end(); ++boundIt) {
MeshStructure elCode;
+// if (i == 13)
+// elCode.setDebugMode(true);
+
elCode.init(boundIt->rankObj);
if (elCode.getCode() != recvCodes[i].getCode()) {
TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
+ /*
+ MSG("----------------------------------------------------------------\n");
+ MSG("CODE %d: MESH SUB OBJ = %d INDEX = %d ITH = %d REVERSE = %d\n",
+ i,
+ boundIt->rankObj.subObj,
+ boundIt->rankObj.el->getIndex(),
+ boundIt->rankObj.ithObj,
+ boundIt->rankObj.reverseMode);
+
+ MSG("RANK CODE: \n");
+ elCode.print();
+ elCode.reset();
+
+ MSG("OTHER CODE: \n");
+ recvCodes[i].print();
+ recvCodes[i].reset();
+ */
bool b = fitElementToMeshCode(recvCodes[i],
boundIt->rankObj.el,
boundIt->rankObj.subObj,
boundIt->rankObj.ithObj,
boundIt->rankObj.elType,
- boundIt->rankObj.reverseMode);
+ boundIt->neighObj.reverseMode);
+
+ // MSG("Mesh changed in bound from %d to %d: %d\n", mpiRank, it->first, b);
if (b)
meshFitTogether = false;
- }
+ } else {
+ /*
+ MSG("----------------------------------------------------------------\n");
+ MSG("CODE %d is equal!\n", i);
+
+ MSG("RANK CODE: \n");
+ elCode.print();
+ elCode.reset();
+
+ MSG("OTHER CODE: \n");
+ recvCodes[i].print();
+ recvCodes[i].reset();
+ */
+ }
i++;
}
@@ -615,6 +662,13 @@ namespace AMDiS {
if (code.empty())
return false;
+// MSG("FIT EL %d!\n", el->getIndex());
+// if (mpiRank == 0) {
+// debug::highlightElementIndexMesh(mesh, 139, "test.vtu");
+// }
+
+// MSG("REVERSE MODE = %d\n", reverseMode);
+
int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
@@ -629,12 +683,19 @@ namespace AMDiS {
traverseFlag |= Mesh::CALL_REVERSE_MODE;
}
+ // MSG("SUB-OBJ (SWAPED): %d %d\n", s1, s2);
+
if (s1 != -1 && s2 != -1) {
+ // MSG("GO TO FIT 1\n");
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
while (elInfo && elInfo->getElement() != el)
elInfo = stack.traverseNext(elInfo);
+ TEST_EXIT_DBG(elInfo->getElement() == el)("This should not happen!\n");
+
+ // MSG("START WITH EL = %d\n", elInfo->getElement()->getIndex());
+
meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
return meshChanged;
}
@@ -643,6 +704,8 @@ namespace AMDiS {
if (code.getNumElements() == 1 && code.isLeafElement())
return false;
+ // MSG("GO TO FIT 2\n");
+
meshChanged = true;
TraverseStack stack;
@@ -666,6 +729,8 @@ namespace AMDiS {
std::swap(child0, child1);
if (s1 != -1) {
+ // MSG("GO TO FIT 3\n");
+
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
@@ -676,6 +741,8 @@ namespace AMDiS {
meshChanged |=
fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
} else {
+ // MSG("GO TO FIT 4\n");
+
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
@@ -702,6 +769,11 @@ namespace AMDiS {
ElInfo *elInfo = stack.getElInfo();
+ // MSG("IN %d\n", elInfo->getElement()->getIndex());
+
+// MeshStructure c = code;
+// c.print(false);
+
bool value = false;
if (!elInfo)
return value;
@@ -709,11 +781,18 @@ namespace AMDiS {
Element *el = elInfo->getElement();
if (code.isLeafElement()) {
+ // MSG("HERE 1\n");
int level = elInfo->getLevel();
do {
elInfo = stack.traverseNext(elInfo);
} while (elInfo && elInfo->getLevel() > level);
+
+// if (elInfo) {
+// MSG("NEXT ELEMENT IS: %d\n", elInfo->getElement()->getIndex());
+// } else {
+// MSG("NO NEXT ELEMENT!\n");
+// }
return value;
}
@@ -722,6 +801,10 @@ namespace AMDiS {
return value;
if (el->isLeaf()) {
+// MSG("HERE 2: %d\n", el->getIndex());
+// if (mpiRank == 0) {
+// debug::highlightElementIndexMesh(mesh, 1592, "test_1592.vtu");
+// }
el->setMark(1);
refineManager->setMesh(el->getMesh());
refineManager->setStack(&stack);
@@ -739,22 +822,27 @@ namespace AMDiS {
}
if (s1 != -1) {
+ // MSG("HERE 3\n");
stack.traverseNext(elInfo);
code.nextElement();
value |= fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
elInfo = stack.getElInfo();
} else {
+ // MSG("HERE 4\n");
do {
elInfo = stack.traverseNext(elInfo);
- } while (elInfo && elInfo->getElement() != child1);
+ } while (elInfo && elInfo->getElement() != child1);
+ // MSG("S1 == -1 -> SWITCHED TO EL: %d\n", elInfo->getElement()->getIndex());
}
TEST_EXIT_DBG(elInfo->getElement() == child1)("This should not happen!\n");
if (s2 != -1) {
+ // MSG("HERE 5\n");
code.nextElement();
value |= fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
} else {
+ // MSG("HERE 6\n");
int level = elInfo->getLevel();
do {
diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc
index d43c4f0f..09442845 100644
--- a/AMDiS/src/parallel/ParallelDebug.cc
+++ b/AMDiS/src/parallel/ParallelDebug.cc
@@ -231,22 +231,25 @@ namespace AMDiS {
// === Print error message if the coordinates are not the same. ===
if (printCoords && !c) {
- std::cout.precision(5);
- std::cout << "[DBG] i = " << i << std::endl;
- std::cout << "[DBG] Rank " << pdb.mpiRank << " from rank " << it->first
+ MSG("[DBG] i = %d\n", i);
+
+ std::stringstream oss;
+ oss.precision(5);
+ oss << "[DBG] Rank " << pdb.mpiRank << " from rank " << it->first
<< " expect coords (";
for (int k = 0; k < dimOfWorld; k++) {
- std::cout << recvCoords[it->first][i][k];
+ oss << recvCoords[it->first][i][k];
if (k + 1 < dimOfWorld)
- std::cout << " / ";
+ oss << " / ";
}
- std::cout << ") received coords (";
+ oss << ") received coords (";
for (int k = 0; k < dimOfWorld; k++) {
- std::cout << (it->second)[i][k];
+ oss << (it->second)[i][k];
if (k + 1 < dimOfWorld)
- std::cout << " / ";
+ oss << " / ";
}
- std::cout << ")" << std::endl;
+ oss << ")";
+ MSG("%s\n", oss.str().c_str());
debug::printInfoByDof(pdb.feSpace, *(pdb.recvDofs[it->first][i]));
}
--
GitLab