#include "MeshStructure.h" #include "MeshStructure_ED.h" #include "PartitionElementData.h" #include "Mesh.h" #include "Element.h" #include "Traverse.h" #include "ElInfo.h" #include "RefinementManager.h" namespace AMDiS { const int MeshStructure::unsignedLongSize = sizeof(unsigned long int) * 8; void MeshStructure::insertElement(bool isLeaf) { // overflow? -> next index if (pos >= unsignedLongSize) { code.push_back(currentCode); pos = 0; currentCode = 0; } // insert element in binary code if (!isLeaf) { unsigned long int one = 1; currentCode += (one << pos); } pos++; nElements++; } void MeshStructure::clear() { currentCode = 0; code.resize(0); pos = 0; nElements = 0; currentElement = 0; } void MeshStructure::init(Mesh *mesh) { clear(); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { insertElement(elInfo->getElement()->isLeaf()); elInfo = stack.traverseNext(elInfo); } commit(); } void MeshStructure::init(BoundaryObject &bound) { FUNCNAME("MeshStructure::init()"); Element *el = bound.el; clear(); int s1 = el->getSubObjOfChild(0, bound.subObj, bound.ithObj, bound.elType); int s2 = el->getSubObjOfChild(1, bound.subObj, bound.ithObj, bound.elType); TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n"); if (debugMode) { MSG("addAlondSide(%d, %d, %d, %d)\n", bound.elIndex, bound.ithObj, bound.elType, bound.reverseMode); MSG("Element is leaf: %d\n", el->isLeaf()); MSG("s1 = %d s2 = %d\n", s1, s2); } if (!el->isLeaf()) { if (s1 == -1) addAlongSide(el->getSecondChild(), bound.subObj, s2, el->getChildType(bound.elType), bound.reverseMode); else if (s2 == -1) addAlongSide(el->getFirstChild(), bound.subObj, s1, el->getChildType(bound.elType), bound.reverseMode); else addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode); } commit(); } void MeshStructure::addAlongSide(Element *el, GeoIndex subObj, int ithObj, int elType, bool reverseOrder) { FUNCNAME("MeshStructure::addAlongSide()"); if (debugMode) { MSG("addAlondSide(%d, %d, %d, %d)\n", el->getIndex(), ithObj, elType, reverseOrder); MSG("Element is leaf: %d\n", el->isLeaf()); } insertElement(el->isLeaf()); if (!el->isLeaf()) { int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType); int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType); if (debugMode) { MSG("Child index %d %d\n", el->getFirstChild()->getIndex(), el->getSecondChild()->getIndex()); MSG("s1 = %d s2 = %d\n", s1, s2); MSG(" \n"); } if (!reverseOrder) { if (s1 != -1) addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder); if (s2 != -1) addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder); } else { if (s2 != -1) addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder); if (s1 != -1) addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder); } } } void MeshStructure::reset() { currentIndex = 0; pos = 0; currentElement = 0; if (code.size() > 0) currentCode = code[0]; else currentCode = 0; } bool MeshStructure::nextElement(MeshStructure *insert) { FUNCNAME("MeshStructure::nextElement()"); if (insert) insert->insertElement(isLeafElement()); pos++; currentElement++; if (currentElement >= nElements) return false; if (pos >= unsignedLongSize) { currentIndex++; TEST_EXIT_DBG(currentIndex < static_cast(code.size())) ("End of structure reached!\n"); pos = 0; currentCode = code[currentIndex]; } else { currentCode >>= 1; } return true; } bool MeshStructure::skipBranch(MeshStructure *insert) { FUNCNAME("MeshStructure::skipBranch()"); if (isLeafElement()) { return nextElement(insert); } else { bool cont = nextElement(insert); cont = skipBranch(insert); // left branch TEST_EXIT_DBG(cont)("Invalid structure!\n"); cont = skipBranch(insert); // righ branch return cont; } } void MeshStructure::merge(MeshStructure *structure1, MeshStructure *structure2, MeshStructure *result) { FUNCNAME("MeshStructure::merge()"); result->clear(); structure1->reset(); structure2->reset(); bool cont = true; while (cont) { bool cont1, cont2; if (structure1->isLeafElement() == structure2->isLeafElement()) { cont1 = structure1->nextElement(result); cont2 = structure2->nextElement(); } else { if (structure1->isLeafElement()) { cont1 = structure1->nextElement(); cont2 = structure2->skipBranch(result); } else { cont1 = structure1->skipBranch(result); cont2 = structure2->nextElement(); } } TEST_EXIT_DBG(cont1 == cont2)("Structures don't match!\n"); cont = cont1; } result->commit(); } void MeshStructure::fitMeshToStructure(Mesh *mesh, RefinementManager *manager, bool checkPartition, bool debugMode) { FUNCNAME("MeshStructure::fitMeshToStructure()"); bool cont = true; // decorate leaf data reset(); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { TEST_EXIT(cont)("unexpected structure code end!\n"); Element *element = elInfo->getElement(); if (isLeafElement()) { TEST_EXIT_DBG(element->isLeaf())("mesh finer than code\n"); } if (element->isLeaf() && !isLeafElement()) { MeshStructure *structure = new MeshStructure(); cont = skipBranch(structure); structure->commit(); bool decorate = true; if (checkPartition) { PartitionElementData *partitionData = dynamic_cast (element->getElementData(PARTITION_ED)); TEST_EXIT_DBG(partitionData)("no partition element data\n"); PartitionStatus status = partitionData->getPartitionStatus(); if (debugMode == false && (status == OUT || status == UNDEFINED)) decorate = false; } if (decorate) { MeshStructure_ED *elData = new MeshStructure_ED(element->getElementData()); elData->setStructure(structure); element->setElementData(elData); } else { delete structure; } } else { cont = nextElement(); } elInfo = stack.traverseNext(elInfo); } // refine mesh bool finished = true; do { finished = true; elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { Element *element = elInfo->getElement(); if (element->getElementData(MESH_STRUCTURE) != NULL) { element->setMark(1); finished = false; } else { element->setMark(0); } elInfo = stack.traverseNext(elInfo); } manager->refineMesh(mesh); } while (!finished); } bool MeshStructure::compare(MeshStructure &other) { return (other.getCode() == code); } }