Commit 9bb9b920 authored by Thomas Witkowski's avatar Thomas Witkowski

Fixed some problems in 3D refinement with reversed mesh traversel stratagey.

parent acfe9a60
......@@ -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);
/****************************************************************************/
......
......@@ -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()");
......
......@@ -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).
......
......@@ -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) {
......
......@@ -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);
......
......@@ -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()
......
......@@ -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);
......
......@@ -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 */