From 5eb8344be15500e3df8c0cbab99195d501d8a09b Mon Sep 17 00:00:00 2001 From: Thomas Witkowski Date: Fri, 11 Feb 2011 13:49:34 +0000 Subject: [PATCH] Fixed 3D refinement problem in parallel computations with mesh repartitioning. --- AMDiS/src/CoarseningManager3d.cc | 168 ++++++++++++++++------ AMDiS/src/Debug.cc | 13 ++ AMDiS/src/Debug.h | 5 + AMDiS/src/ElInfo3d.cc | 44 +----- AMDiS/src/FiniteElemSpace.cc | 18 ++- AMDiS/src/FiniteElemSpace.h | 8 ++ AMDiS/src/Mesh.cc | 28 ---- AMDiS/src/Mesh.h | 7 - AMDiS/src/RCNeighbourList.cc | 4 +- AMDiS/src/RefinementManager3d.cc | 195 +++++++++++++++++++++----- AMDiS/src/RefinementManager3d.h | 22 ++- AMDiS/src/Tetrahedron.cc | 6 +- AMDiS/src/Tetrahedron.h | 3 + AMDiS/src/parallel/MeshDistributor.cc | 69 ++++++--- AMDiS/src/parallel/MeshDistributor.h | 7 +- AMDiS/src/parallel/ParallelDebug.cc | 13 +- 16 files changed, 423 insertions(+), 187 deletions(-) diff --git a/AMDiS/src/CoarseningManager3d.cc b/AMDiS/src/CoarseningManager3d.cc index bbf40633..d4c2bded 100644 --- a/AMDiS/src/CoarseningManager3d.cc +++ b/AMDiS/src/CoarseningManager3d.cc @@ -19,36 +19,40 @@ #include "RCNeighbourList.h" #include "FixVec.h" #include "DOFIndexed.h" +#include "Debug.h" namespace AMDiS { - void CoarseningManager3d::coarsenFunction(ElInfo *el_info) + void CoarseningManager3d::coarsenFunction(ElInfo *elInfo) { + FUNCNAME("CoarseningManager3d::coarsenFunction()"); + Tetrahedron *el = - dynamic_cast(const_cast(el_info->getElement())); + dynamic_cast(const_cast(elInfo->getElement())); + // If element must not be coarsend, return. if (el->getMark() >= 0) - return; // el must not be coarsend, return :-( + return; + // Single leaves don't get coarsened. if (el->isLeaf()) - return; // single leaves don't get coarsened + return; if (el->getChild(0)->getMark() >= 0 || el->getChild(1)->getMark() >= 0) { - // one of the children must not be coarsend; return :-( + // One of the children must not be coarsend, so return. el->setMark(0); return; } if (!(el->getChild(0)->isLeaf()) || !(el->getChild(1)->isLeaf())) { - // one of the children is not a leaf element; try again later on + // One of the children is not a leaf element. Try again later on. doMore = true; return; } DegreeOfFreedom *edge[2]; int n_neigh, bound = 0; - ElInfo *elinfo = el_info; /****************************************************************************/ /* get a list for storing all elements at the coarsening edge and fill it */ @@ -60,24 +64,95 @@ namespace AMDiS { /* give the refinement edge the right orientation */ /****************************************************************************/ - if (el->getDof(0,0) < el->getDof(1,0)) { - edge[0] = const_cast(el->getDof(0)); - edge[1] = const_cast(el->getDof(1)); + if (el->getDof(0, 0) < el->getDof(1, 0)) { + edge[0] = const_cast(el->getDof(0)); + edge[1] = const_cast(el->getDof(1)); } else { - edge[1] = const_cast(el->getDof(0)); - edge[0] = const_cast(el->getDof(1)); + edge[1] = const_cast(el->getDof(0)); + edge[0] = const_cast(el->getDof(1)); } coarsenList.setElement(0, el, true); n_neigh = 1; coarsenList.setOppVertex(0, 0, 0); - coarsenList.setElType(0, el_info->getType()); + coarsenList.setElType(0, elInfo->getType()); bound = false; - if (getCoarsenPatch(elinfo, edge, 0, coarsenList, &n_neigh)) { - getCoarsenPatch(elinfo, edge, 1, coarsenList, &n_neigh); + if (getCoarsenPatch(elInfo, edge, 0, coarsenList, &n_neigh)) { + getCoarsenPatch(elInfo, edge, 1, coarsenList, &n_neigh); bound = true; } + + +#if HAVE_PARALLEL_DOMAIN_AMDIS + Element *otherEl = NULL; + int otherEdge = -1; + FixRefinementPatch::getOtherEl(stack, &otherEl, otherEdge); + + // === If the refinement edge must be fixed, add also the other part of this === + // === edge to the refinement patch. === + + if (otherEl) { + // TODO: Remove these two lines and make something more meaningful!! + el->setMark(0); + return; + + TraverseStack stack2; + ElInfo *elInfo2 = + stack2.traverseFirstOneMacro(mesh, otherEl->getIndex(), -1, + Mesh::CALL_EVERY_EL_PREORDER | + Mesh::FILL_NEIGH | + Mesh::FILL_BOUND); + + bool foundEdge = false; + + while (elInfo2) { + Element *el2 = elInfo2->getElement(); + for (int i = 0; i < 6; i++) { + DofEdge edge2 = elInfo2->getElement()->getEdge(i); + if (edge2.first == *(edge[0]) && edge2.second == *(edge[1]) && !el2->isLeaf()) { + + if (!el2->getChild(0)->isLeaf() || !el2->getChild(1)->isLeaf()) { + int edgeNo0 = el->getEdgeOfChild(0, i, elInfo2->getType()); + int edgeNo1 = el->getEdgeOfChild(1, i, elInfo2->getType()); + + bool refineChildFirst = + !(i > 0 && + (edgeNo0 >= 0 && !el2->getChild(0)->isLeaf()) || + (edgeNo1 >= 0 && !el2->getChild(1)->isLeaf())); + + if (refineChildFirst) { + // TODO: WAS SOLL ICH NUR MACHEN??? + el->setMark(0); + return; + } + } else { + coarsenList.setElType(n_neigh, elInfo2->getType()); + coarsenList.setElement(n_neigh, elInfo2->getElement(), true); + n_neigh++; + + foundEdge = true; + + TraverseStack *tmpStack = stack; + stack = &stack2; + + if (getCoarsenPatch(elInfo2, edge, 0, coarsenList, &n_neigh)) { + getCoarsenPatch(elInfo2, edge, 1, coarsenList, &n_neigh); + bound = true; + } + stack = tmpStack; + break; + } + } + } + + elInfo2 = stack2.traverseNext(elInfo2); + } + + TEST_EXIT_DBG(foundEdge)("Should not happen!\n"); + } +#endif + coarsenList.fillNeighbourRelations(n_neigh, bound); /****************************************************************************/ @@ -95,8 +170,8 @@ namespace AMDiS { while (edge[0] != NULL) { coarsenList.periodicSplit(edge, next_edge, - &n_neigh, &n_neigh_periodic, - periodicList); + &n_neigh, &n_neigh_periodic, + periodicList); coarsenPatch(periodicList, n_neigh_periodic, bound); @@ -119,22 +194,19 @@ namespace AMDiS { Tetrahedron *el = dynamic_cast(const_cast(coarsenList.getElement(index))); Tetrahedron *child[2]; - Tetrahedron *neigh; - int dir, el_type, i, node, opp_v; child[0] = dynamic_cast(const_cast(el->getChild(0))); child[1] = dynamic_cast(const_cast(el->getChild(1))); - el_type = coarsenList.getType(index); + int el_type = coarsenList.getType(index); /****************************************************************************/ /* Information about patch neighbours is still valid! But edge and face */ /* dof's in a common face of patch neighbours have to be removed */ /****************************************************************************/ - for (dir = 0; dir < 2; dir++) { - neigh = + for (int dir = 0; dir < 2; dir++) { + Tetrahedron *neigh = dynamic_cast(const_cast(coarsenList.getNeighbourElement(index, dir))); - opp_v = coarsenList.getOppVertex(index, dir); if (!neigh || neigh->isLeaf()) { /****************************************************************************/ @@ -142,11 +214,11 @@ namespace AMDiS { /****************************************************************************/ if (mesh->getNumberOfDofs(EDGE)) { - node = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][0][dir]; + int node = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][0][dir]; mesh->freeDof(const_cast(child[0]->getDof(node)), EDGE); } if (mesh->getNumberOfDofs(FACE)) { - node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir]; + int node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir]; mesh->freeDof(const_cast(child[0]->getDof(node)), FACE); node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir]; mesh->freeDof(const_cast(child[1]->getDof(node)), FACE); @@ -160,14 +232,14 @@ namespace AMDiS { /****************************************************************************/ if (mesh->getNumberOfDofs(FACE)) { - node = mesh->getNode(FACE); + int node = mesh->getNode(FACE); mesh->freeDof(const_cast(child[0]->getDof(node)), FACE); } if (mesh->getNumberOfDofs(CENTER)) { - node = mesh->getNode(CENTER); - for (i = 0; i < 2; i++) + int node = mesh->getNode(CENTER); + for (int i = 0; i < 2; i++) mesh->freeDof(const_cast(child[i]->getDof(node)), CENTER); } @@ -191,7 +263,7 @@ namespace AMDiS { /****************************************************************************/ /* get_coarse_patch: gets the patch for coarsening starting on element */ - /* el_info->el in direction of neighbour [3-dir]; returns 1 if a boundary */ + /* elInfo->el in direction of neighbour [3-dir]; returns 1 if a boundary */ /* reached and 0 if we come back to the starting element. */ /* */ /* if NEIGH_IN_EL we only can find the complete coarsening patch if the */ @@ -205,16 +277,13 @@ namespace AMDiS { /* the starting element before we return */ /****************************************************************************/ - bool CoarseningManager3d::getCoarsenPatch(ElInfo *el_info, + bool CoarseningManager3d::getCoarsenPatch(ElInfo *elInfo, DegreeOfFreedom *edge[2], int dir, RCNeighbourList &coarsenList, int *n_neigh) { - FUNCNAME("CoarseningManager3d::getCoarsenPatch"); - ElInfo *neigh_info; - Tetrahedron *el, *neigh; - int i, j, k, opp_v, edge_no; + FUNCNAME("CoarseningManager3d::getCoarsenPatch()"); static unsigned char next_el[6][2] = {{3,2}, {1,3}, @@ -223,15 +292,16 @@ namespace AMDiS { {0,2}, {0,1}}; - el = dynamic_cast(const_cast(el_info->getElement())); - neigh = - dynamic_cast(const_cast(el_info->getNeighbour(3 - dir))); + Tetrahedron *el = + dynamic_cast(const_cast(elInfo->getElement())); + Tetrahedron *neigh = + dynamic_cast(const_cast(elInfo->getNeighbour(3 - dir))); if (neigh == NULL) return true; - opp_v = el_info->getOppVertex(3 - dir); + int opp_v = elInfo->getOppVertex(3 - dir); + ElInfo *neigh_info = stack->traverseNeighbour3d(elInfo, 3 - dir); - neigh_info = stack->traverseNeighbour3d(el_info, 3 - dir); TEST_EXIT_DBG(neigh == neigh_info->getElement()) ("neigh %d and neigh_info->el %d are not identical\n", neigh->getIndex(), neigh_info->getElement()->getIndex()); @@ -244,6 +314,7 @@ namespace AMDiS { coarsenList.setElType(*n_neigh, neigh_info->getType()); int n_vertices = mesh->getGeo(VERTEX); + int i, j, k, edge_no; while (neigh != el) { for (j = 0; j < n_vertices; j++) @@ -261,10 +332,15 @@ namespace AMDiS { if (mesh->associated(neigh->getDof(k, 0), edge[1][0])) break; - TEST_EXIT_DBG(j < n_vertices && k < n_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 < n_vertices) + ("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(k < n_vertices) + ("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)); } edge_no = Tetrahedron::edgeOfDofs[j][k]; coarsenList.setCoarsePatch(*n_neigh, edge_no == 0); @@ -366,6 +442,12 @@ namespace AMDiS { // === And now start to coarsen the patch: remove dof's of the coarsening edge. === + WorldVector c; + FiniteElemSpace *feSpace = FiniteElemSpace::provideFeSpace(mesh); + mesh->getDofIndexCoords(el->getChild(0)->getDof(3), feSpace, c); + + // MSG("Free dof %p %f %f %f\n", el->getChild(0)->getDof(3), c[0], c[1], c[2]); + mesh->freeDof(const_cast(el->getChild(0)->getDof(3)), VERTEX); mesh->incrementNumberOfVertices(-1); diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc index 47eb7e39..23d3cc3d 100644 --- a/AMDiS/src/Debug.cc +++ b/AMDiS/src/Debug.cc @@ -74,6 +74,19 @@ namespace AMDiS { } + void colorEdgeInMesh(FiniteElemSpace *feSpace, + Element *el, + int localEdgeNo, + std::string filename) + { + DOFVector tmp(feSpace, "tmp"); + tmp.set(0.0); + tmp[el->getEdge(localEdgeNo).first] = 1.0; + tmp[el->getEdge(localEdgeNo).second] = 1.0; + VtkWriter::writeFile(tmp, filename); + } + + void writeElementIndexMesh(Mesh *mesh, std::string filename, int level) { FUNCNAME("debug::writeElementIndexMesh()"); diff --git a/AMDiS/src/Debug.h b/AMDiS/src/Debug.h index 9d0d141b..55e912f0 100644 --- a/AMDiS/src/Debug.h +++ b/AMDiS/src/Debug.h @@ -66,6 +66,11 @@ namespace AMDiS { */ void writeDofIndexMesh(FiniteElemSpace *feSpace); + void colorEdgeInMesh(FiniteElemSpace *feSpace, + Element *el, + int localEdgeNo, + std::string filename); + /** \brief * Creates a vtu file where all elements in the mesh are colored by the global * element indices. diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc index 62deae09..7a7b55a2 100644 --- a/AMDiS/src/ElInfo3d.cc +++ b/AMDiS/src/ElInfo3d.cc @@ -483,7 +483,7 @@ namespace AMDiS { normal[1] /= det; normal[2] /= det; } else { - MSG("not implemented for DIM_OF_WORLD = %d in 3d\n", dimOfWorld); + MSG("Not implemented for DIM_OF_WORLD = %d in 3D!\n", dimOfWorld); } return det; @@ -593,37 +593,11 @@ namespace AMDiS { if ((nb = const_cast((*neigh_old)[cv[i]]))) { TEST_EXIT_DBG(nb->getChild(0))("nonconforming triangulation\n"); - if (elOld->getIndex() == 15376 && ichild == 0) { - MSG("NEIGH %d OF EL is %d\n", i, nb->getIndex()); - MSG("EL HAS DOFS: %d %d %d %d\n", - elOld->getDof(0, 0), - elOld->getDof(1, 0), - elOld->getDof(2, 0), - elOld->getDof(3, 0)); - MSG("NEIGH-EL HAS DOFS: %d %d %d %d\n", - nb->getDof(0, 0), - nb->getDof(1, 0), - nb->getDof(2, 0), - nb->getDof(3, 0)); - } - - int k; for (k = 0; k < 2; k++) { /* look at both childs of old neighbour */ nbk = const_cast(nb->getChild(k)); - if (elOld->getIndex() == 15376 && ichild == 0) { - MSG("TEST k = %d, nbk = %d, el = %d\n", k, - nbk->getDof(0,0), elOld->getDof(ichild,0)); - MSG("TEST PTR nbk = %p, el = %p\n", - nbk->getDof(0), elOld->getDof(ichild)); - } - if (nbk->getDof(0) == elOld->getDof(ichild)) { - if (elOld->getIndex() == 15376 && ichild == 0) { - MSG("DRIN!\n"); - } - /* opp. vertex */ dof = const_cast(nb->getDof(elInfoOld->oppVertex[cv[i]])); @@ -641,17 +615,13 @@ namespace AMDiS { } (*neigh_local)[i] = nbk->getChild(0); oppVertex[i] = 3; - if (elOld->getIndex() == 15376 && ichild == 0) { - MSG("BREAK 1\n"); - } + break; } } else { if (dof != nbk->getDof(2)) { ov = -1; - if (elOld->getIndex() == 15376 && ichild == 0) { - MSG("BREAK 2\n"); - } + break; } ov = 2; @@ -663,18 +633,12 @@ namespace AMDiS { (*neigh_local)[i] = nbk; oppVertex[i] = ov; - if (elOld->getIndex() == 15376 && ichild == 0) { - MSG("BREAK 3\n"); - } + break; } } /* end for k */ - if (elOld->getIndex() == 15376 && ichild == 0) { - MSG("k = %d ov = %d\n", k, ov); - } - // === Test if periodic. === diff --git a/AMDiS/src/FiniteElemSpace.cc b/AMDiS/src/FiniteElemSpace.cc index 3df125c4..46c76d0f 100644 --- a/AMDiS/src/FiniteElemSpace.cc +++ b/AMDiS/src/FiniteElemSpace.cc @@ -89,7 +89,7 @@ namespace AMDiS { } - FiniteElemSpace *FiniteElemSpace::provideFeSpace(DOFAdmin *admin, + FiniteElemSpace* FiniteElemSpace::provideFeSpace(DOFAdmin *admin, const BasisFunction *basFcts, Mesh *mesh, std::string name_) @@ -104,6 +104,22 @@ namespace AMDiS { } +#if DEBUG + FiniteElemSpace* FiniteElemSpace::provideFeSpace(Mesh *mesh) + { + FUNCNAME("FiniteElemSpace::provideFeSpace()"); + + for (unsigned int i = 0; i < feSpaces.size(); i++) + if (feSpaces[i]->mesh == mesh) + return feSpaces[i]; + + ERROR_EXIT("FE space not found!\n"); + + return NULL; + } +#endif + + int FiniteElemSpace::calcMemoryUsage() { int result = sizeof(FiniteElemSpace); diff --git a/AMDiS/src/FiniteElemSpace.h b/AMDiS/src/FiniteElemSpace.h index c0d2a442..a899e75a 100644 --- a/AMDiS/src/FiniteElemSpace.h +++ b/AMDiS/src/FiniteElemSpace.h @@ -49,6 +49,14 @@ namespace AMDiS { Mesh *mesh, std::string name = ""); +#if DEBUG + /// For debugging it may be useful to get some FE space for a given mesh at a + /// position in code where it is not possible to access the FE space directly. The + /// function assumes that there is only one FE space defined for the mesh. + + static FiniteElemSpace *provideFeSpace(Mesh *mesh); +#endif + /// Destructor. ~FiniteElemSpace(); diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc index 5befd72e..0cb73339 100644 --- a/AMDiS/src/Mesh.cc +++ b/AMDiS/src/Mesh.cc @@ -1426,34 +1426,6 @@ namespace AMDiS { } - void Mesh::computeMatrixFillin(const FiniteElemSpace *feSpace, - std::vector &nnzInRow, int &overall, int &average) - { - std::map dofCounter; - - TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL); - ElementDofIterator elDofIter(feSpace); - while (elInfo) { - elDofIter.reset(elInfo->getElement()); - do { - DegreeOfFreedom dof = elDofIter.getDof(); - if (dofCounter.count(dof) == 0) { - dofCounter[dof] = 1; - } else { - dofCounter[dof]++; - } - } while (elDofIter.next()); - elInfo = stack.traverseNext(elInfo); - } - - overall = 0; - for (std::map::iterator it = dofCounter.begin(); - it != dofCounter.end(); ++it) - overall += it->second * 15; - } - - int Mesh::calcMemoryUsage() { int result = sizeof(Mesh); diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h index 9edd9306..0a86767c 100644 --- a/AMDiS/src/Mesh.h +++ b/AMDiS/src/Mesh.h @@ -595,13 +595,6 @@ namespace AMDiS { /// void clearMacroFileInfo(); - /** \brief - * Traverse this mesh to compute the number of non zero elements the assembled - * matrix will have in each row. - */ - void computeMatrixFillin(const FiniteElemSpace *feSpace, - std::vector &nnzInRow, int &overall, int &average); - /// int calcMemoryUsage(); diff --git a/AMDiS/src/RCNeighbourList.cc b/AMDiS/src/RCNeighbourList.cc index e91f8904..aac8d0ee 100644 --- a/AMDiS/src/RCNeighbourList.cc +++ b/AMDiS/src/RCNeighbourList.cc @@ -320,6 +320,8 @@ namespace AMDiS { int *n_neigh_periodic, RCNeighbourList &periodicList) { + FUNCNAME("RCNeighbourList::periodicSplit()"); + *n_neigh_periodic = 0; int count = 0; int n_neigh_old = *n_neigh; @@ -356,7 +358,7 @@ namespace AMDiS { if (edge[0]) { TEST_EXIT_DBG((dof0 == edge[0] && dof1 == edge[1]) || (dof1 == edge[0] && dof0 == edge[1])) - ("invalid macro file?\n"); + ("Invalid macro file?\n"); } } diff --git a/AMDiS/src/RefinementManager3d.cc b/AMDiS/src/RefinementManager3d.cc index df8d3441..d3cd7035 100644 --- a/AMDiS/src/RefinementManager3d.cc +++ b/AMDiS/src/RefinementManager3d.cc @@ -29,11 +29,15 @@ namespace AMDiS { + FixRefinementPatch::ConnectedEdges FixRefinementPatch::connectedEdges(0); + void RefinementManager3d::bisectTetrahedron(RCNeighbourList& refineList, int index, DegreeOfFreedom* dof[3], DegreeOfFreedom *edge[2]) { + FUNCNAME("RefinementManager3d::bisectTetrahedron()"); + Tetrahedron *el = dynamic_cast(const_cast(refineList.getElement(index))); Tetrahedron *child[2]; @@ -169,7 +173,6 @@ namespace AMDiS { } if (!neigh || neigh->isLeaf()) { - /****************************************************************************/ /* get new dof's in the midedge of the face of el and for the two midpoints*/ /* of the sub-faces. If face is an interior face those pointers have to be */ @@ -190,7 +193,7 @@ namespace AMDiS { node1 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir]; (const_cast(el->getSecondChild()))->setDof(node1, mesh->getDof(FACE)); } - } else { /* if (!neigh || !neigh->child[0]) */ + } else { /****************************************************************************/ /* interior face and neighbour has been refined, look for position at the */ /* refinement edge */ @@ -244,24 +247,24 @@ namespace AMDiS { } - void RefinementManager3d::newCoordsFct(ElInfo *el_info, RCNeighbourList &refineList) + void RefinementManager3d::newCoordsFct(ElInfo *elInfo, RCNeighbourList &refineList) { FUNCNAME("RefinementManager3d::newCoordsFct()"); - Element *el = el_info->getElement(); + Element *el = elInfo->getElement(); DegreeOfFreedom *edge[2]; - ElInfo *elinfo = el_info; + ElInfo *elinfo = elInfo; int dow = Global::getGeo(WORLD); - Projection *projector = el_info->getProjection(0); + Projection *projector = elInfo->getProjection(0); if (!projector || projector->getType() != VOLUME_PROJECTION) - projector = el_info->getProjection(4); + projector = elInfo->getProjection(4); if (el->getFirstChild() && projector && (!el->isNewCoordSet())) { WorldVector *new_coord = new WorldVector; for (int j = 0; j < dow; j++) - (*new_coord)[j] = (el_info->getCoord(0)[j] + el_info->getCoord(1)[j]) * 0.5; + (*new_coord)[j] = (elInfo->getCoord(0)[j] + elInfo->getCoord(1)[j]) * 0.5; projector->project(*new_coord); @@ -271,11 +274,11 @@ namespace AMDiS { /* get the refinement patch */ /****************************************************************************/ refineList.setElement(0, el); - refineList.setElType(0, el_info->getType()); + refineList.setElType(0, elInfo->getType()); int n_neigh = 1; for (int i = 0; i < 2; i++) - edge[i] = const_cast(el_info->getElement()->getDof(i)); + edge[i] = const_cast(elInfo->getElement()->getDof(i)); if (getRefinePatch(&elinfo, edge, 0, refineList, &n_neigh)) { // Domain's boundary was reached while looping around the refinement edge. @@ -322,6 +325,8 @@ namespace AMDiS { RCNeighbourList &refineList, int n_neigh, bool bound) { + FUNCNAME("RefinementManager3d::refinePatch()"); + Tetrahedron *el = dynamic_cast(const_cast(refineList.getElement(0))); /* first element in the list */ @@ -393,24 +398,24 @@ namespace AMDiS { } - int RefinementManager3d::getRefinePatch(ElInfo **el_info, - DegreeOfFreedom *edge[2], - int direction, - RCNeighbourList &refineList, - int *n_neigh) + bool RefinementManager3d::getRefinePatch(ElInfo **elInfo, + DegreeOfFreedom *edge[2], + int direction, + RCNeighbourList &refineList, + int *n_neigh) { FUNCNAME("RefinementManager3d::getRefinePatch()"); int localNeighbour = 3 - direction; Tetrahedron *el = - dynamic_cast(const_cast((*el_info)->getElement())); + dynamic_cast(const_cast((*elInfo)->getElement())); - if ((*el_info)->getNeighbour(localNeighbour) == NULL) - return 1; + if ((*elInfo)->getNeighbour(localNeighbour) == NULL) + return true; - int oppVertex = (*el_info)->getOppVertex(localNeighbour); - int testIndex = (*el_info)->getNeighbour(localNeighbour)->getIndex(); - ElInfo *neighInfo = stack->traverseNeighbour3d((*el_info), localNeighbour); + int oppVertex = (*elInfo)->getOppVertex(localNeighbour); + int testIndex = (*elInfo)->getNeighbour(localNeighbour)->getIndex(); + ElInfo *neighInfo = stack->traverseNeighbour3d((*elInfo), localNeighbour); int neighElType = neighInfo->getType(); TEST_EXIT_DBG(neighInfo->getElement()->getIndex() == testIndex) @@ -447,11 +452,15 @@ namespace AMDiS { if (mesh->indirectlyAssociated(neigh->getDof(edgeDof1, 0), edge[1][0])) break; - 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), + TEST_EXIT_DBG(edgeDof0 < vertices) + ("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", + edge[1][0], neigh->getIndex(), neigh->getDof(0, 0), + neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0)); } } @@ -567,9 +576,9 @@ namespace AMDiS { } if (neigh == el) { - (*el_info) = neighInfo; + (*elInfo) = neighInfo; - return 0; + return false; } @@ -597,20 +606,22 @@ namespace AMDiS { ("Should not happen!\n"); } while (neighInfo->getElement() != el); - (*el_info) = neighInfo; + (*elInfo) = neighInfo; - return 1; + return true; } - ElInfo* RefinementManager3d::refineFunction(ElInfo* el_info) + ElInfo* RefinementManager3d::refineFunction(ElInfo* elInfo) { FUNCNAME("RefinementManager3d::refineFunction()"); - Element *el = el_info->getElement(); + Element *el = elInfo->getElement(); if (el->getMark() <= 0) - return el_info; + return elInfo; + + // MSG("Refine: %d\n", elInfo->getElement()->getIndex()); int bound = false; DegreeOfFreedom *edge[2]; @@ -618,12 +629,12 @@ namespace AMDiS { // === Get memory for a list of all elements at the refinement edge. === RCNeighbourList refineList(mesh->getMaxEdgeNeigh()); - refineList.setElType(0, el_info->getType()); + refineList.setElType(0, elInfo->getType()); refineList.setElement(0, el); int n_neigh = 1; - if (el_info->getProjection(0) && - el_info->getProjection(0)->getType() == VOLUME_PROJECTION) + if (elInfo->getProjection(0) && + elInfo->getProjection(0)->getType() == VOLUME_PROJECTION) newCoords = true; @@ -637,15 +648,65 @@ namespace AMDiS { edge[0] = const_cast(el->getDof(1)); } +#if HAVE_PARALLEL_DOMAIN_AMDIS + Element *otherEl = NULL; + int otherEdge = -1; + FixRefinementPatch::getOtherEl(stack, &otherEl, otherEdge); +#endif // === Traverse and refine the refinement patch. ==== - if (getRefinePatch(&el_info, edge, 0, refineList, &n_neigh)) { + if (getRefinePatch(&elInfo, edge, 0, refineList, &n_neigh)) { // Domain's boundary was reached while looping around the refinement edge - getRefinePatch(&el_info, edge, 1, refineList, &n_neigh); + getRefinePatch(&elInfo, edge, 1, refineList, &n_neigh); bound = true; } +#if HAVE_PARALLEL_DOMAIN_AMDIS + // === If the refinement edge must be fixed, add also the other part of this === + // === edge to the refinement patch. === + + if (otherEl) { + TraverseStack stack2; + ElInfo *elInfo2 = + stack2.traverseFirstOneMacro(mesh, otherEl->getIndex(), -1, + Mesh::CALL_LEAF_EL | + Mesh::FILL_NEIGH | + Mesh::FILL_BOUND); + + bool foundEdge = false; + + while (elInfo2) { + for (int i = 0; i < 6; i++) { + 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; + } + + stack = tmpStack; + break; + } + } + + elInfo2 = stack2.traverseNext(elInfo2); + } + + TEST_EXIT_DBG(foundEdge)("Should not happen!\n"); + } +#endif + + // fill neighbour information inside the patch in the refinement list refineList.fillNeighbourRelations(n_neigh, bound); @@ -712,6 +773,66 @@ namespace AMDiS { stack->update(); - return el_info; + return elInfo; + } + + + void FixRefinementPatch::getOtherEl(TraverseStack *stack, + Element **otherEl, + int &otherEdge) + { + FUNCNAME("FixRefinementPatch::getOtherEl()"); + + if (!FixRefinementPatch::connectedEdges.empty()) { + // === Get stack of current traverse. === + std::vector elInfos; + std::vector infos; + int stackUsed = stack->getStackData(elInfos, infos); + int checkIndex = stackUsed; + int localEdgeNo = 0; + ElInfo *elInfo = elInfos[stackUsed]; + + // === Calculate the refinement edge number on the macro element level. === + + for (int i = 0; i < elInfo->getLevel(); i++) { + TEST_EXIT_DBG(checkIndex >= 1)("Should not happen!\n"); + + int parentType = elInfos[checkIndex - 1]->getType(); + int ithParentChild = 0; + if (elInfos[checkIndex - 1]->getElement()->getChild(1) == + elInfos[checkIndex]->getElement()) + ithParentChild = 1; + localEdgeNo = Tetrahedron::childEdge[parentType][ithParentChild][localEdgeNo]; + + if (localEdgeNo >= 4) { + localEdgeNo = -1; + break; + } + + checkIndex--; + } + + // If the refinement edge is part of an edge on the macro level, we must check + // if the refinement edge is part of an edge that must be fixed. + if (localEdgeNo >= 0) { + TEST_EXIT_DBG(elInfos[checkIndex]->getLevel() == 0)("Should not happen!\n"); + TEST_EXIT_DBG(elInfos[checkIndex]->getElement()->getIndex() == + elInfo->getMacroElement()->getIndex())("Should not happen!\n"); + TEST_EXIT_DBG(localEdgeNo <= 5)("Should not happen!\n"); + + for (unsigned int i = 0; i < FixRefinementPatch::connectedEdges.size(); i++) { + if (FixRefinementPatch::connectedEdges[i].first.first->getIndex() == + elInfos[checkIndex]->getElement()->getIndex() && + FixRefinementPatch::connectedEdges[i].first.second == + localEdgeNo) { + // Okay, we have found that this edge must be fixed. + *otherEl = FixRefinementPatch::connectedEdges[i].second.first; + otherEdge = FixRefinementPatch::connectedEdges[i].second.second; + + break; + } + } + } + } } } diff --git a/AMDiS/src/RefinementManager3d.h b/AMDiS/src/RefinementManager3d.h index 57d4a968..38fd7dba 100644 --- a/AMDiS/src/RefinementManager3d.h +++ b/AMDiS/src/RefinementManager3d.h @@ -56,11 +56,11 @@ namespace AMDiS { * \param[in] direction Determines the direction of the first neighbour * */ - int getRefinePatch(ElInfo **el_info, - DegreeOfFreedom *edge[2], - int direction, - RCNeighbourList &refineList, - int *n_neigh); + bool 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, @@ -77,6 +77,18 @@ namespace AMDiS { void fillPatchConnectivity(RCNeighbourList &refineList, int index); }; + class FixRefinementPatch + { + public: + typedef std::pair EdgeInEl; + typedef std::pair EdgesMap; + typedef std::vector ConnectedEdges; + + static ConnectedEdges connectedEdges; + + static void getOtherEl(TraverseStack *stack, Element **otherEl, int &otherEdge); + }; + } #endif // AMDIS_REFINEMENT_MANAGER_3D_H diff --git a/AMDiS/src/Tetrahedron.cc b/AMDiS/src/Tetrahedron.cc index 5372118a..f2b40d1a 100644 --- a/AMDiS/src/Tetrahedron.cc +++ b/AMDiS/src/Tetrahedron.cc @@ -32,9 +32,9 @@ namespace AMDiS { {{0,2,3,4},{1,2,3,4}}, {{0,2,3,4},{1,2,3,4}}}; - const int Tetrahedron::childEdge[3][2][6] = {{{1,2,0,5,5,4},{4,3,0,5,5,4}}, - {{1,2,0,5,4,5},{3,4,0,5,4,5}}, - {{1,2,0,5,4,5},{3,4,0,5,4,5}}}; + const int Tetrahedron::childEdge[3][2][6] = {{{1,2,0,5,5,4}, {4,3,0,5,5,4}}, + {{1,2,0,5,4,5}, {3,4,0,5,4,5}}, + {{1,2,0,5,4,5}, {3,4,0,5,4,5}}}; const unsigned char Tetrahedron::adjacentChild[2][2] = {{0,1}, {1,0}}; diff --git a/AMDiS/src/Tetrahedron.h b/AMDiS/src/Tetrahedron.h index 0f77b8bd..9adb1c65 100644 --- a/AMDiS/src/Tetrahedron.h +++ b/AMDiS/src/Tetrahedron.h @@ -182,6 +182,9 @@ namespace AMDiS { * new edge 2 is half of old edge 0, * new edges 4,5 are really new edges, and value is different: * childEdge[][][4,5] = index of same edge in other child + * + * Note: el_type, i.e., the first index of the array, defines the element type + * of the parent elements. */ static const int childEdge[3][2][6]; diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index 37790f16..b8b0783c 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -46,6 +46,7 @@ #include "MeshStructure.h" #include "ProblemVec.h" #include "ProblemInstat.h" +#include "RefinementManager3d.h" #include "Debug.h" namespace AMDiS { @@ -201,8 +202,6 @@ namespace AMDiS { ParallelDebug::printBoundaryInfo(*this); #endif - MSG("HERE 1\n"); - for (deque::iterator it = mesh->firstMacroElement(); it != mesh->endOfMacroElements(); ++it) { for (int i = 0; i < mesh->getGeo(NEIGH); i++) { @@ -529,7 +528,8 @@ namespace AMDiS { elInfo = stack.traverseNext(elInfo); } - bool meshIsValid = true; + + FixRefinementPatch::connectedEdges.clear(); for (std::set::iterator it = allEdges.begin(); it != allEdges.end(); ++it) { vector& edgeEls = elObjects.getElements(*it); @@ -538,12 +538,18 @@ namespace AMDiS { ("No edge %d/%d in elObjDB!\n", it->first, it->second); std::set el0, el1; - for (unsigned int i = 0; i < edgeEls.size(); i++) - if (rankMacroEls.count(edgeEls[i].elIndex)) + std::map edgeNoInEl; + + for (unsigned int i = 0; i < edgeEls.size(); i++) { + if (rankMacroEls.count(edgeEls[i].elIndex)) { if (el0.empty()) el0.insert(edgeEls[i].elIndex); else el1.insert(edgeEls[i].elIndex); + + edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject; + } + } TEST_EXIT_DBG(!el0.empty())("Should not happen!\n"); @@ -569,18 +575,42 @@ namespace AMDiS { } while (found); if (!el1.empty()) { - MSG("Unvalid 3D mesh with %d elements on this edge!\n", edgeEls.size()); - meshIsValid = false; - } - } + MSG("Unvalid 3D mesh with %d (%d - %d) elements on this edge!\n", + edgeEls.size(), el0.size(), el1.size()); - if (!meshIsValid) { - MSG("Error: mesh is not a valid 3D mesh on this rank!\n"); - } + for (std::set::iterator elIt0 = el0.begin(); elIt0 != el0.end(); ++elIt0) { + for (std::set::iterator elIt1 = el1.begin(); elIt1 != el1.end(); ++elIt1) { + + TEST_EXIT_DBG(macroElIndexMap.count(*elIt0))("Should not happen!\n"); + TEST_EXIT_DBG(macroElIndexMap.count(*elIt1))("Should not happen!\n"); + TEST_EXIT_DBG(edgeNoInEl.count(*elIt0))("Should not happen!\n"); + TEST_EXIT_DBG(edgeNoInEl.count(*elIt1))("Should not happen!\n"); - int foundError = !meshIsValid; - mpi::globalAdd(foundError); - TEST_EXIT(foundError == 0)("Error found on at least on rank!\n"); + pair edge0 = + make_pair(macroElIndexMap[*elIt0], edgeNoInEl[*elIt0]); + pair edge1 = + make_pair(macroElIndexMap[*elIt1], edgeNoInEl[*elIt1]); + + DofEdge dofEdge0 = edge0.first->getEdge(edge0.second); + DofEdge dofEdge1 = edge1.first->getEdge(edge1.second); + + WorldVector c0, c1; + 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]); + + TEST_EXIT_DBG(dofEdge0 == dofEdge1)("Should noth happen!\n"); + + FixRefinementPatch::connectedEdges.push_back(make_pair(edge0, edge1)); + FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0)); + } + } + } + } } @@ -651,7 +681,7 @@ namespace AMDiS { } - void MeshDistributor::checkMeshChange() + void MeshDistributor::checkMeshChange(bool tryRepartition) { FUNCNAME("MeshDistributor::checkMeshChange()"); @@ -695,7 +725,7 @@ namespace AMDiS { elCode.init(it->rankObj); MeshManipulation mm(feSpace); - meshChanged |= !(mm.fitElementToMeshCode(elCode, it->neighObj)); + meshChanged |= mm.fitElementToMeshCode(elCode, it->neighObj); } } else { if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || @@ -748,7 +778,8 @@ namespace AMDiS { INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n", MPI::Wtime() - first); - if (repartitioningAllowed && + if (tryRepartition && + repartitioningAllowed && nMeshChangesAfterLastRepartitioning >= repartitionIthChange) { repartitionMesh(); nMeshChangesAfterLastRepartitioning = 0; @@ -1977,7 +2008,7 @@ namespace AMDiS { { FUNCNAME("MeshDistributor::serialize()"); - checkMeshChange(); + checkMeshChange(false); partitioner->serialize(out); diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h index 2b14fced..59fbf050 100644 --- a/AMDiS/src/parallel/MeshDistributor.h +++ b/AMDiS/src/parallel/MeshDistributor.h @@ -108,8 +108,13 @@ namespace AMDiS { * all ranks. Furthermore the function \ref updateLocalGlobalNumbering() is called * to update the dof numberings and mappings on all rank due to the new mesh * structure. + * + * \param[in] tryRepartition If this parameter is true, repartitioning may be + * done. This depends on several other parameters. If + * the parameter is false, the mesh is only checked + * and adapted but never repartitioned. */ - void checkMeshChange(); + void checkMeshChange(bool tryRepartition = true); /** \brief * Test, if the mesh consists of macro elements only. The mesh partitioning of diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc index 7219cf1e..a4bc7ea0 100644 --- a/AMDiS/src/parallel/ParallelDebug.cc +++ b/AMDiS/src/parallel/ParallelDebug.cc @@ -144,6 +144,8 @@ namespace AMDiS { { FUNCNAME("ParallelDebug::testPeriodicBoundary()"); + return; + // === 1. check: All periodic DOFs should have at least a correct number === // === of periodic associations. === @@ -179,6 +181,7 @@ namespace AMDiS { // === The boundary DOFs are checked only on the zero rank. === if (pdb.mpiRank == 0) { + // Stores to each rank the periodic DOF mappings of this rank. std::map rankToMaps; PeriodicDofMap dofMap = pdb.periodicDof; rankToMaps[0] = dofMap; @@ -598,8 +601,9 @@ namespace AMDiS { cMap[c] = elInfo->getElement()->getDof(i, 0); } else { if (cMap[c] != elInfo->getElement()->getDof(i, 0)) { - MSG("[DBG] Found two DOFs %d and %d with the same coords %f %f!\n", - cMap[c], elInfo->getElement()->getDof(i, 0), c[0], c[1]); + MSG("[DBG] Found two DOFs %d and %d with the same coords %f %f %f!\n", + cMap[c], elInfo->getElement()->getDof(i, 0), + c[0], c[1], mesh->getDim() == 3 ? c[2] : 0.0); foundError = 1; } } @@ -727,6 +731,11 @@ namespace AMDiS { { FUNCNAME("ParallelDebug::printBoundaryInfo()"); + int tmp = 1; + GET_PARAMETER(0, "parallel->debug->print boundary info", "%d", &tmp); + if (tmp <= 0) + return; + for (InteriorBoundary::iterator it(pdb.myIntBoundary); !it.end(); ++it) { MSG("Rank owned boundary with rank %d: \n", it.getRank()); MSG(" ranks obj-ind: %d sub-obj: %d ith-obj: %d\n", -- GitLab