//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.
#include "CoarseningManager3d.h"
#include "Mesh.h"
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "Traverse.h"
#include "MacroElement.h"
#include "RCNeighbourList.h"
#include "FixVec.h"
#include "DOFIndexed.h"
#include "Debug.h"
namespace AMDiS {
void CoarseningManager3d::coarsenFunction(ElInfo *elInfo)
{
FUNCNAME("CoarseningManager3d::coarsenFunction()");
Tetrahedron *el =
dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getElement()));
// If element must not be coarsend, return.
if (el->getMark() >= 0)
return;
// Single leaves don't get coarsened.
if (el->isLeaf())
return;
if (el->getChild(0)->getMark() >= 0 || el->getChild(1)->getMark() >= 0) {
// 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.
doMore = true;
return;
}
DegreeOfFreedom *edge[2];
int n_neigh, bound = 0;
/****************************************************************************/
/* get a list for storing all elements at the coarsening edge and fill it */
/****************************************************************************/
RCNeighbourList coarsenList(mesh->getMaxEdgeNeigh());
coarsenList.setCoarseningManager(this);
/****************************************************************************/
/* give the refinement edge the right orientation */
/****************************************************************************/
if (el->getDof(0, 0) < el->getDof(1, 0)) {
edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(0));
edge[1] = const_cast<DegreeOfFreedom*>(el->getDof(1));
} else {
edge[1] = const_cast<DegreeOfFreedom*>(el->getDof(0));
edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(1));
}
coarsenList.setElement(0, el, true);
n_neigh = 1;
coarsenList.setOppVertex(0, 0, 0);
coarsenList.setElType(0, elInfo->getType());
bound = false;
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);
/****************************************************************************/
/* check wether we can coarsen the patch or not */
/****************************************************************************/
// ==========================================================================
// === check for periodic boundary ==========================================
// ==========================================================================
if (coarsenList.doCoarsePatch(n_neigh)) {
int n_neigh_periodic;
DegreeOfFreedom *next_edge[2];
RCNeighbourList periodicList;
while (edge[0] != NULL) {
coarsenList.periodicSplit(edge, next_edge,
&n_neigh, &n_neigh_periodic,
periodicList);
coarsenPatch(periodicList, n_neigh_periodic, bound);
edge[0] = next_edge[0];
edge[1] = next_edge[1];
}
}
}
/*****************************************************************************/
/* coarsenTetrahedron: coarses a single element of the coarsening patch; */
/* dofs in the interior of the element are removed and dofs in the faces of */
/* the element are removed if neighbour has been coarsend or if the face */
/* is part of the domains boundary */
/*****************************************************************************/
void CoarseningManager3d::coarsenTetrahedron(RCNeighbourList &coarsenList,
int index)
{
Tetrahedron *el =
dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getElement(index)));
Tetrahedron *child[2];
child[0] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(el->getChild(0)));
child[1] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(el->getChild(1)));
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 (int dir = 0; dir < 2; dir++) {
Tetrahedron *neigh =
dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getNeighbourElement(index, dir)));
if (!neigh || neigh->isLeaf()) {
/****************************************************************************/
/* boundary face or neigh has been coarsend: free the dof's in the face */
/****************************************************************************/
if (mesh->getNumberOfDofs(EDGE)) {
int node = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][0][dir];
mesh->freeDof(const_cast<int*>(child[0]->getDof(node)), EDGE);
}
if (mesh->getNumberOfDofs(FACE)) {
int node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir];
mesh->freeDof(const_cast<int*>(child[0]->getDof(node)), FACE);
node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir];
mesh->freeDof(const_cast<int*>(child[1]->getDof(node)), FACE);
}
}
}
/****************************************************************************/
/* finally remove the interior dof's: in the common face of child[0] and */
/* child[1] and at the two barycenter */
/****************************************************************************/
if (mesh->getNumberOfDofs(FACE)) {
int node = mesh->getNode(FACE);
mesh->freeDof(const_cast<int*>(child[0]->getDof(node)), FACE);
}
if (mesh->getNumberOfDofs(CENTER)) {
int node = mesh->getNode(CENTER);
for (int i = 0; i < 2; i++)
mesh->freeDof(const_cast<int*>(child[i]->getDof(node)), CENTER);
}
/****************************************************************************/
/* get new data on parent and transfer data from children to parent */
/****************************************************************************/
el->coarsenElementData(child[0], child[1], el_type);
el->setFirstChild(NULL);
el->setSecondChild(NULL);
mesh->freeElement(child[0]);
mesh->freeElement(child[1]);
el->incrementMark();
mesh->incrementNumberOfLeaves(-1);
mesh->incrementNumberOfElements(-2);
}
/****************************************************************************/
/* get_coarse_patch: gets the patch for coarsening starting on element */
/* 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 */
/* can be coarsend; otherwise neighbour information is not valid for */
/* parents; in such situation we stop looping around the edge and return 0 */
/* */
/* if !NEIGH_IN_EL we complete the loop also in the case of a incompatible */
/* coarsening patch since then all marks of patch elements are reset by */
/* do_coarse_patch() and this minimizes calls of traverse_neighbour(); */
/* if we reach a boundary while looping around the edge we loop back to */
/* the starting element before we return */
/****************************************************************************/
bool CoarseningManager3d::getCoarsenPatch(ElInfo *elInfo,
DegreeOfFreedom *edge[2],
int dir,
RCNeighbourList &coarsenList,
int *n_neigh)
{
FUNCNAME("CoarseningManager3d::getCoarsenPatch()");
static unsigned char next_el[6][2] = {{3,2},
{1,3},
{1,2},
{0,3},
{0,2},
{0,1}};
Tetrahedron *el =
dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getElement()));
Tetrahedron *neigh =
dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getNeighbour(3 - dir)));
if (neigh == NULL)
return true;
int opp_v = elInfo->getOppVertex(3 - dir);
ElInfo *neigh_info = stack->traverseNeighbour3d(elInfo, 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());
/****************************************************************************/
/* we have to go back to the starting element via opp_v values */
/* correct information is produce by get_neigh_on_patch() */
/****************************************************************************/
coarsenList.setOppVertex(*n_neigh, 0, opp_v);
coarsenList.setElement(*n_neigh, neigh);
coarsenList.setElType(*n_neigh, neigh_info->getType());
int nVertices = mesh->getGeo(VERTEX);
while (neigh != el) {
// === Find the coarsening edge on the current neighbour element. ===
int i, j, k;
// === First, try to identify the edge DOFs directly. ===
for (j = 0; j < nVertices; j++)
if (neigh->getDof(j) == edge[0])
break;
for (k = 0; k < nVertices; k++)
if (neigh->getDof(k) == edge[1])
break;
// === If one of the edge DOFs was not found, try to make use of periodic ===
// === DOF associations. First, use the direct associations between DOFs. ===
// === If this will not work, continue with testing on indirect ===
// === associations. ===
if (j >= nVertices) {
for (j = 0; j < nVertices; j++)
if (mesh->associated(neigh->getDof(j, 0), edge[0][0]))
break;
if (j >= nVertices)
for (j = 0; j < nVertices; j++)
if (mesh->indirectlyAssociated(neigh->getDof(j, 0), edge[0][0]))
break;
TEST_EXIT_DBG(j < nVertices)
("Process element %d: DOF %d not found on element %d with nodes (%d %d %d %d)\n",
el->getIndex(), edge[0][0], neigh->getIndex(), neigh->getDof(0, 0),
neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
}
if (k >= nVertices) {
for (k = 0; k < nVertices; k++)
if (mesh->associated(neigh->getDof(k, 0), edge[1][0]))
break;
if (k >= nVertices)
for (k = 0; k < nVertices; k++)
if (mesh->indirectlyAssociated(neigh->getDof(k, 0), edge[1][0]))
break;
TEST_EXIT_DBG(k < nVertices)
("Process element %d: DOF %d not found on element %d with nodes (%d %d %d %d)\n",
el->getIndex(), 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[j][k];
coarsenList.setCoarsePatch(*n_neigh, edgeNo == 0);
// === Get the direction of the next neighbour. ===
if (next_el[edgeNo][0] != opp_v)
i = next_el[edgeNo][0];
else
i = next_el[edgeNo][1];
++*n_neigh;
opp_v = neigh_info->getOppVertex(i);
neigh =
dynamic_cast<Tetrahedron*>(const_cast<Element*>(neigh_info->getNeighbour(i)));
if (neigh) {
neigh_info = stack->traverseNeighbour3d(neigh_info, i);
TEST_EXIT_DBG(neigh == neigh_info->getElement())
("neigh %d and neigh_info->el %d are not identical\n",
neigh->getIndex(), neigh_info->getElement()->getIndex());
/****************************************************************************/
/* we have to go back to the starting element via opp_v values */
/* correct information is produce by get_neigh_on_patch() */
/****************************************************************************/
coarsenList.setOppVertex(*n_neigh, 0, opp_v);
coarsenList.setElement(*n_neigh, neigh);
coarsenList.setElType(*n_neigh, neigh_info->getType());
} else {
break;
}
}
if (neigh == el)
return false;
/****************************************************************************/
/* the domain's boundary is reached; loop back to the starting el */
/****************************************************************************/
int i = *n_neigh - 1;
opp_v = coarsenList.getOppVertex(i, 0);
do {
TEST_EXIT_DBG(neigh_info->getNeighbour(opp_v) && i > 0)
("while looping back domains boundary was reached or i == 0\n");
opp_v = coarsenList.getOppVertex(i--, 0);
neigh_info = stack->traverseNeighbour3d(neigh_info, opp_v);
} while (neigh_info->getElement() != el);
return true;
}
/****************************************************************************/
/* coarse_patch: first rebuild the dofs on the parents then do restriction */
/* of data (if possible) and finally coarsen the patch elements */
/****************************************************************************/
void CoarseningManager3d::coarsenPatch(RCNeighbourList &coarsenList,
int n_neigh,
int bound)
{
FUNCNAME("CoarseningManager3d::coarsenPatch()");
Tetrahedron *el =
dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getElement(0)));
DegreeOfFreedom *dof = NULL;
TEST_EXIT_DBG(el)("No element!\n");
TEST_EXIT_DBG(el->getChild(0))("No child in element!\n");
if (mesh->getNumberOfDofs(EDGE)) {
// === Get dof for coarsening edge. ===
int node = mesh->getNode(EDGE);
if (!(dof = const_cast<int*>(el->getDof(node))))
dof = mesh->getDof(EDGE);
}
if (mesh->getNumberOfDofs(EDGE) ||
mesh->getNumberOfDofs(FACE) ||
mesh->getNumberOfDofs(CENTER)) {
for (int i = 0; i < n_neigh; i++)
coarsenList.addDOFParent(i, dof);
}
// === Restrict dof vectors to the parents on the patch. ===
int nrAdmin = mesh->getNumberOfDOFAdmin();
for (int iadmin = 0; iadmin < nrAdmin; iadmin++) {
std::list<DOFIndexedBase*>::iterator it;
DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDofAdmin(iadmin));
std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
for (it = admin->beginDOFIndexed(); it != end; ++it)
(*it)->coarseRestrict(coarsenList, n_neigh);
}
// === And now start to coarsen the patch: remove dof's of the coarsening edge. ===
mesh->freeDof(const_cast<int*>(el->getChild(0)->getDof(3)), VERTEX);
mesh->incrementNumberOfVertices(-1);
if (mesh->getNumberOfDofs(EDGE)) {
int node = mesh->getNode(EDGE) + 2;
mesh->freeDof(const_cast<int*>(el->getChild(0)->getDof(node)), EDGE);
mesh->freeDof(const_cast<int*>(el->getChild(1)->getDof(node)), EDGE);
}
if (coarsenList.getElement(0)->isNewCoordSet())
coarsenList.getElement(0)->eraseNewCoord();
for (int i = 0; i < n_neigh; i++) {
coarsenList.getElement(i)->setNewCoord(NULL);
coarsenTetrahedron(coarsenList, i);
}
// === If the coarsening edge is an interior edge there are n_neigh + 1 edges ===
// === and 2 * n_neigh + 1 faces removed; if it is a boundary edge it is one ===
// === more edge and one more face. ===
if (bound) {
mesh->incrementNumberOfEdges(-(n_neigh + 2));
mesh->incrementNumberOfFaces(-(2 * n_neigh + 1));
} else {
mesh->incrementNumberOfEdges(-(n_neigh + 1));
mesh->incrementNumberOfFaces(-(2 * n_neigh));
}
}
}