diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index c540d1b6a9cf03107ac88e96b994fb363a7d006a..425d54d4d5ef43dae11fd426d4f1fe4ef0758667 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -134,7 +134,7 @@ namespace AMDiS {
return;
}
-
+
// Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
// only for macro meshes, so it will not work yet if the mesh is already refined
// in some way.
@@ -540,10 +540,11 @@ namespace AMDiS {
void MeshDistributor::checkMeshChange()
{
- FUNCNAME("MeshDistributor::checkMeshChange()");
+ FUNCNAME("MeshDistributor::checkMeshChange()");
double first = MPI::Wtime();
+
// === If mesh has not been changed on all ranks, return. ===
int recvAllValues = 0;
@@ -557,6 +558,8 @@ namespace AMDiS {
// === adapted to the new mesh structure. ===
do {
+ bool meshChanged = false;
+
// To check the interior boundaries, the ownership of the boundaries is not
// important. Therefore, we add all boundaries to one boundary container.
RankToBoundMap allBound;
@@ -573,7 +576,20 @@ namespace AMDiS {
for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it) {
if (it.getRank() == mpiRank) {
- WARNING("Na, du weisst schon!\n");
+ if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) ||
+ (mesh->getDim() == 3 && it->rankObj.subObj == FACE)) {
+ MeshStructure elCode;
+ elCode.init(it->rankObj);
+
+ MeshManipulation meshManipulation(feSpace);
+ meshChanged |=
+ !(meshManipulation.startFitElementToMeshCode(elCode,
+ it->neighObj.el,
+ it->neighObj.subObj,
+ it->neighObj.ithObj,
+ it->neighObj.elType,
+ it->neighObj.reverseMode));
+ }
} else {
if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) ||
(mesh->getDim() == 3 && it->rankObj.subObj == FACE))
@@ -587,7 +603,7 @@ namespace AMDiS {
MSG("Run checkAndAdaptBoundary ...\n");
#endif
- bool meshChanged = checkAndAdaptBoundary(allBound);
+ meshChanged |= checkAndAdaptBoundary(allBound);
// === Check on all ranks if at least one rank's mesh has changed. ===
@@ -671,15 +687,17 @@ namespace AMDiS {
if (elCode.getCode() != recvCodes[i].getCode()) {
TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
- bool b = startFitElementToMeshCode(recvCodes[i],
- boundIt->rankObj.el,
- boundIt->rankObj.subObj,
- boundIt->rankObj.ithObj,
- boundIt->rankObj.elType,
- boundIt->rankObj.reverseMode);
+ MeshManipulation meshManipulation(feSpace);
+ bool b =
+ meshManipulation.startFitElementToMeshCode(recvCodes[i],
+ boundIt->rankObj.el,
+ boundIt->rankObj.subObj,
+ boundIt->rankObj.ithObj,
+ boundIt->rankObj.elType,
+ boundIt->rankObj.reverseMode);
if (b)
- meshFitTogether = false;
+ meshFitTogether = false;
}
i++;
@@ -690,225 +708,6 @@ namespace AMDiS {
}
- bool MeshDistributor::startFitElementToMeshCode(MeshStructure &code,
- Element *el,
- GeoIndex subObj,
- int ithObj,
- int elType,
- bool reverseMode)
- {
- FUNCNAME("MeshDistributor::startFitElementToMeshCode()");
-
- TEST_EXIT_DBG(el)("No element given!\n");
-
- // If the code is empty, the element does not matter and the function can
- // return without chaning the mesh.
- if (code.empty())
- return false;
-
- // s0 and s1 are the number of the edge/face in both child of the element,
- // which contain the edge/face the function has to traverse through. If the
- // edge/face is not contained in one of the children, s0 or s1 is -1.
- int s0 = el->getSubObjOfChild(0, subObj, ithObj, elType);
- int s1 = el->getSubObjOfChild(1, subObj, ithObj, elType);
-
- TEST_EXIT_DBG(s0 != -1 || s1 != -1)("This should not happen!\n");
-
- bool meshChanged = false;
- Flag traverseFlag =
- Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;
-
- // Test for reverse mode, in which the left and right children of elements
- // are flipped.
- if (reverseMode)
- traverseFlag |= Mesh::CALL_REVERSE_MODE;
-
-
- // === If the edge/face is contained in both children. ===
-
- if (s0 != -1 && s1 != -1) {
- // Create traverse stack and traverse within the mesh until the element,
- // which should be fitted to the mesh structure code, is reached.
- 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");
-
- meshChanged = fitElementToMeshCode(code, stack, subObj, ithObj, reverseMode);
- return meshChanged;
- }
-
-
- // === The edge/face is contained in only on of the both children. ===
-
- if (el->isLeaf()) {
-
- // If element is leaf and code contains only one leaf element, we are finished.
- if (code.getNumElements() == 1 && code.isLeafElement())
- return false;
-
- // Create traverse stack and traverse the mesh to the element.
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
- while (elInfo && elInfo->getElement() != el)
- elInfo = stack.traverseNext(elInfo);
-
- TEST_EXIT_DBG(elInfo)("This should not happen!\n");
-
- // Code is not leaf, therefore refine the element.
- el->setMark(1);
- refineManager->setMesh(el->getMesh());
- refineManager->setStack(&stack);
- refineManager->refineFunction(elInfo);
- meshChanged = true;
- }
-
- Element *child0 = el->getFirstChild();
- Element *child1 = el->getSecondChild();
- if (reverseMode) {
- swap(s0, s1);
- swap(child0, child1);
- }
-
- // === We know that the edge/face is contained in only one of the children. ===
- // === Therefore, traverse the mesh to this children and fit this element ===
- // === To the mesh structure code. ===
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
-
- if (s0 != -1) {
- while (elInfo && elInfo->getElement() != child0)
- elInfo = stack.traverseNext(elInfo);
-
- meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
- } else {
- while (elInfo && elInfo->getElement() != child1)
- elInfo = stack.traverseNext(elInfo);
-
- meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
- }
-
-
- return meshChanged;
- }
-
-
- bool MeshDistributor::fitElementToMeshCode(MeshStructure &code,
- TraverseStack &stack,
- GeoIndex subObj,
- int ithObj,
- bool reverseMode)
- {
- FUNCNAME("MeshDistributor::fitElementToMeshCode()");
-
-
- // === Test if there are more elements in stack to check with the code. ===
-
- ElInfo *elInfo = stack.getElInfo();
- if (!elInfo)
- return false;
-
-
- // === Test if code contains a leaf element. If this is the case, the ===
- // === current element is finished. Traverse the mesh to the next ===
- // === coarser element. ===
-
- if (code.isLeafElement()) {
- int level = elInfo->getLevel();
-
- do {
- elInfo = stack.traverseNext(elInfo);
- } while (elInfo && elInfo->getLevel() > level);
-
- return false;
- }
-
-
- bool meshChanged = false;
- Element *el = elInfo->getElement();
-
-
- // === If element is leaf (and the code is not), refine the element. ===
-
- if (el->isLeaf()) {
- TEST_EXIT_DBG(elInfo->getLevel() < 255)("This should not happen!\n");
-
- el->setMark(1);
- refineManager->setMesh(el->getMesh());
- refineManager->setStack(&stack);
- refineManager->refineFunction(elInfo);
- meshChanged = true;
- }
-
-
- // === Continue fitting the mesh structure code to the children of the ===
- // === current element. ===
-
- int s0 = el->getSubObjOfChild(0, subObj, ithObj, elInfo->getType());
- int s1 = el->getSubObjOfChild(1, subObj, ithObj, elInfo->getType());
- Element *child0 = el->getFirstChild();
- Element *child1 = el->getSecondChild();
- if (reverseMode) {
- swap(s0, s1);
- swap(child0, child1);
- }
-
-
- // === Traverse left child. ===
-
- if (s0 != -1) {
- // The edge/face is contained in the left child, therefore fit this
- // child to the mesh structure code.
-
- stack.traverseNext(elInfo);
- code.nextElement();
- meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
- elInfo = stack.getElInfo();
- } else {
- // The edge/face is not contained in the left child. Hence we need
- // to traverse through all subelements of the left child until we get
- // the second child of the current element.
-
- do {
- elInfo = stack.traverseNext(elInfo);
- } while (elInfo && elInfo->getElement() != child1);
-
- TEST_EXIT_DBG(elInfo != NULL)("This should not happen!\n");
- }
-
- TEST_EXIT_DBG(elInfo->getElement() == child1)
- ("Got wrong child with idx = %d! Searched for child idx = %d\n",
- elInfo->getElement()->getIndex(), child1->getIndex());
-
-
- // === Traverse right child. ===
-
- if (s1 != -1) {
- // The edge/face is contained in the right child, therefore fit this
- // child to the mesh structure code.
-
- code.nextElement();
- meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
- } else {
- // The edge/face is not contained in the right child. Hence we need
- // to traverse through all subelements of the right child until we have
- // finished traversing the current element with all its subelements.
-
- int level = elInfo->getLevel();
-
- do {
- elInfo = stack.traverseNext(elInfo);
- } while (elInfo && elInfo->getLevel() > level);
- }
-
-
- return meshChanged;
- }
-
-
void MeshDistributor::serialize(ostream &out, DofContainer &data)
{
int vecSize = data.size();
@@ -1428,10 +1227,7 @@ namespace AMDiS {
mesh->getDofIndexCoords(it->first.first, feSpace, c0);
mesh->getDofIndexCoords(it->first.second, feSpace, c1);
- // MSG("CREATE BOUNDARY FOR DOF MAP: %d (%.3f %.3f %.3f)<-> %d (%.3f %.3f %.3f)\n", it->first.first,
- // c0[0], c0[1], c0[2], it->first.second, c1[0], c1[1], c1[2]);
ElementObjectData& perDofEl0 = elObjects.getElementsInRank(it->first.first)[mpiRank];
- // MSG("DATA: %d %d %d\n", perDofEl0.elIndex, VERTEX, perDofEl0.ithObject);
for (map<int, ElementObjectData>::iterator elIt = elObjects.getElementsInRank(it->first.second).begin();
elIt != elObjects.getElementsInRank(it->first.second).end(); ++elIt) {
@@ -1611,26 +1407,6 @@ namespace AMDiS {
stdMpi.recv(recvBounds.boundary);
stdMpi.startCommunication<int>(MPI_INT);
- /*
- for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
- rankIt != periodicBoundary.boundary.end(); ++rankIt) {
-
- for (unsigned int j = 0; j < rankIt->second.size(); j++) {
-
- MSG("%-ith boundary with rank %d\n", j, rankIt->first);
- MSG(" bound: el = %d subobj = %d ithobj = %d\n",
- periodicBoundary.boundary[rankIt->first][j].rankObj.elIndex,
- periodicBoundary.boundary[rankIt->first][j].rankObj.subObj,
- periodicBoundary.boundary[rankIt->first][j].rankObj.ithObj);
- MSG(" neigh: el = %d subobj = %d ithobj = %d\n",
- periodicBoundary.boundary[rankIt->first][j].neighObj.elIndex,
- periodicBoundary.boundary[rankIt->first][j].neighObj.subObj,
- periodicBoundary.boundary[rankIt->first][j].neighObj.ithObj);
- }
- }
-
- */
-
for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
rankIt != periodicBoundary.boundary.end(); ++rankIt) {
@@ -1874,7 +1650,9 @@ namespace AMDiS {
bound.neighObj.el->getVertexDofs(feSpace, bound.neighObj, dofs1);
bound.neighObj.el->getNonVertexDofs(feSpace, bound.neighObj, dofs1);
- TEST_EXIT_DBG(dofs0.size() == dofs1.size())("Should not happen!\n");
+ TEST_EXIT_DBG(dofs0.size() == dofs1.size())
+ ("Number of DOFs does not fit together: %d %d\n",
+ dofs0.size(), dofs1.size());
BoundaryType type = bound.type;
@@ -1969,6 +1747,9 @@ namespace AMDiS {
for (unsigned int i = 0; i < dofs.size(); i++) {
DegreeOfFreedom globalDof = mapLocalGlobalDofs[*dofs[i]];
+ TEST_EXIT_DBG(periodicDofAssociations.count(globalDof))
+ ("Should hot happen!\n");
+
for (std::set<BoundaryType>::iterator perAscIt = periodicDofAssociations[globalDof].begin();
perAscIt != periodicDofAssociations[globalDof].end(); ++perAscIt)
if (*perAscIt >= -3) {
@@ -2038,6 +1819,19 @@ namespace AMDiS {
}
+ DegreeOfFreedom MeshDistributor::mapGlobalToLocal(DegreeOfFreedom dof)
+ {
+ FUNCNAME("MeshDistributor::mapGlobalToLocal()");
+
+ for (DofMapping::iterator it = mapLocalGlobalDofs.begin();
+ it != mapLocalGlobalDofs.end(); ++it)
+ if (it->second == dof)
+ return it->first;
+
+ return -1;
+ }
+
+
void MeshDistributor::serialize(ostream &out)
{
FUNCNAME("MeshDistributor::serialize()");
diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h
index d5d64372769d0a20e45ee336d8b075f353995b23..d6926e26e0c499b6f8add4541ef6be15db343a1f 100644
--- a/AMDiS/src/parallel/MeshDistributor.h
+++ b/AMDiS/src/parallel/MeshDistributor.h
@@ -147,6 +147,8 @@ namespace AMDiS {
return mapLocalGlobalDofs[dof];
}
+ DegreeOfFreedom mapGlobalToLocal(DegreeOfFreedom dof);
+
/// Maps a local dof to its local index.
inline DegreeOfFreedom mapLocalToDofIndex(DegreeOfFreedom dof)
{
@@ -160,7 +162,7 @@ namespace AMDiS {
}
/// Returns for a global dof index its periodic mapping for a given boundary type.
- inline int getPeriodicMapping(BoundaryType type, int globalDofIndex)
+ inline int getPeriodicMapping(int globalDofIndex, BoundaryType type)
{
FUNCNAME("MeshDistributor::getPeriodicMapping()");
@@ -175,13 +177,17 @@ namespace AMDiS {
/// associations, i.e., the boundary types the DOF is associated to, for this DOF.
inline std::set<BoundaryType>& getPerDofAssociations(int globalDofIndex)
{
+ TEST_EXIT_DBG(periodicDofAssociations.count(globalDofIndex))
+ ("Should not happen!\n");
+
return periodicDofAssociations[globalDofIndex];
}
/// Returns true, if the DOF (global index) is a periodic DOF.
inline bool isPeriodicDof(int globalDofIndex)
{
- return (periodicDofAssociations.count(globalDofIndex) > 0);
+ return (periodicDofAssociations.count(globalDofIndex) > 0 &&
+ periodicDofAssociations[globalDofIndex].size() > 0);
}
/// Returns true, if the DOF (global index) is a periodic DOF for the given
@@ -326,50 +332,6 @@ namespace AMDiS {
// Removes all periodic boundaries from a given boundary map.
void removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap);
- /** \brief
- * Starts the procedure to fit a given edge/face of an element with a mesh
- * structure code. This functions prepares some data structures and call
- * then \ref fitElementToMeshCode, that mainly refines the element such that
- * it fits to the mesh structure code.
- *
- * \param[in] code The mesh structure code to which the edge/face of
- * an element must be fitted.
- * \param[in] el Pointer to the element.
- * \param[in] subObj Defines whether an edge or a face must be fitted.
- * \param[in] ithObj Defines which edge/face must be fitted.
- * \param[in] elType Element type of the element (only important in 3D).
- * \param[in] reverseMode Defines, whether the mesh structure code is given
- * in reverse mode, i.e., left and right children where
- * changed when the code was created.
- */
- bool startFitElementToMeshCode(MeshStructure &code,
- Element *el,
- GeoIndex subObj,
- int ithObj,
- int elType,
- bool reverseMode);
-
- /** \brief
- * Recursively fits a given mesh structure code to an edge/face of an element.
- * This function is always initialy called from \ref startFitElementToMeshCode.
- *
- * \param[in] code The mesh structure code which is used to fit an
- * edge/face of an element.
- * \param[in] stack A traverse stack object. The upper most element in this
- * stack must be used for fitting the mesh structure code
- * at the current position.
- * \param[in] subObj Defines whether an edge or a face must be fitted.
- * \param[in] ithObj Defines which edge/face must be fitted.
- * \param[in] reverseMode Defines, whether the mesh structure code is given
- * in reverse mode, i.e., left and right children where
- * changed when the code was created.
- */
- bool fitElementToMeshCode(MeshStructure &code,
- TraverseStack &stack,
- GeoIndex subObj,
- int ithObj,
- bool reverseMode);
-
/// Writes a vector of dof pointers to an output stream.
void serialize(ostream &out, DofContainer &data);
diff --git a/AMDiS/src/parallel/MeshManipulation.cc b/AMDiS/src/parallel/MeshManipulation.cc
index d5040cd79f0bc87fce2e7ec6eafc1213ed9fc823..166e2f0a76a828c793e7e6808180282a7a7638d8 100644
--- a/AMDiS/src/parallel/MeshManipulation.cc
+++ b/AMDiS/src/parallel/MeshManipulation.cc
@@ -12,12 +12,36 @@
#include "parallel/MeshManipulation.h"
#include "Mesh.h"
+#include "MeshStructure.h"
#include "BasisFunction.h"
#include "Traverse.h"
#include "Debug.h"
namespace AMDiS {
+ MeshManipulation::MeshManipulation(FiniteElemSpace *space)
+ : feSpace(space),
+ mesh(space->getMesh())
+ {
+ switch (mesh->getDim()) {
+ case 2:
+ refineManager = new RefinementManager2d();
+ break;
+ case 3:
+ refineManager = new RefinementManager3d();
+ break;
+ default:
+ ERROR_EXIT("invalid dim!\n");
+ }
+ }
+
+
+ MeshManipulation::~MeshManipulation()
+ {
+ delete refineManager;
+ }
+
+
void MeshManipulation::deleteDoubleDofs(std::set<MacroElement*>& newMacroEl,
ElementObjects &objects)
{
@@ -176,5 +200,223 @@ namespace AMDiS {
mesh->freeDof(const_cast<DegreeOfFreedom*>(it->first), VERTEX);
}
+
+ bool MeshManipulation::startFitElementToMeshCode(MeshStructure &code,
+ Element *el,
+ GeoIndex subObj,
+ int ithObj,
+ int elType,
+ bool reverseMode)
+ {
+ FUNCNAME("MeshManipulation::startFitElementToMeshCode()");
+
+ TEST_EXIT_DBG(el)("No element given!\n");
+
+ // If the code is empty, the element does not matter and the function can
+ // return without chaning the mesh.
+ if (code.empty())
+ return false;
+
+ // s0 and s1 are the number of the edge/face in both child of the element,
+ // which contain the edge/face the function has to traverse through. If the
+ // edge/face is not contained in one of the children, s0 or s1 is -1.
+ int s0 = el->getSubObjOfChild(0, subObj, ithObj, elType);
+ int s1 = el->getSubObjOfChild(1, subObj, ithObj, elType);
+
+ TEST_EXIT_DBG(s0 != -1 || s1 != -1)("This should not happen!\n");
+
+ bool meshChanged = false;
+ Flag traverseFlag =
+ Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;
+
+ // Test for reverse mode, in which the left and right children of elements
+ // are flipped.
+ if (reverseMode)
+ traverseFlag |= Mesh::CALL_REVERSE_MODE;
+
+
+ // === If the edge/face is contained in both children. ===
+
+ if (s0 != -1 && s1 != -1) {
+ // Create traverse stack and traverse within the mesh until the element,
+ // which should be fitted to the mesh structure code, is reached.
+ 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");
+
+ meshChanged = fitElementToMeshCode(code, stack, subObj, ithObj, reverseMode);
+ return meshChanged;
+ }
+
+
+ // === The edge/face is contained in only on of the both children. ===
+
+ if (el->isLeaf()) {
+
+ // If element is leaf and code contains only one leaf element, we are finished.
+ if (code.getNumElements() == 1 && code.isLeafElement())
+ return false;
+
+ // Create traverse stack and traverse the mesh to the element.
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
+ while (elInfo && elInfo->getElement() != el)
+ elInfo = stack.traverseNext(elInfo);
+
+ TEST_EXIT_DBG(elInfo)("This should not happen!\n");
+
+ // Code is not leaf, therefore refine the element.
+ el->setMark(1);
+ refineManager->setMesh(el->getMesh());
+ refineManager->setStack(&stack);
+ refineManager->refineFunction(elInfo);
+ meshChanged = true;
+ }
+
+ Element *child0 = el->getFirstChild();
+ Element *child1 = el->getSecondChild();
+ if (reverseMode) {
+ swap(s0, s1);
+ swap(child0, child1);
+ }
+
+ // === We know that the edge/face is contained in only one of the children. ===
+ // === Therefore, traverse the mesh to this children and fit this element ===
+ // === To the mesh structure code. ===
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
+
+ if (s0 != -1) {
+ while (elInfo && elInfo->getElement() != child0)
+ elInfo = stack.traverseNext(elInfo);
+
+ meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
+ } else {
+ while (elInfo && elInfo->getElement() != child1)
+ elInfo = stack.traverseNext(elInfo);
+
+ meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
+ }
+
+
+ return meshChanged;
+ }
+
+
+ bool MeshManipulation::fitElementToMeshCode(MeshStructure &code,
+ TraverseStack &stack,
+ GeoIndex subObj,
+ int ithObj,
+ bool reverseMode)
+ {
+ FUNCNAME("MeshManipulation::fitElementToMeshCode()");
+
+
+ // === Test if there are more elements in stack to check with the code. ===
+
+ ElInfo *elInfo = stack.getElInfo();
+ if (!elInfo)
+ return false;
+
+
+ // === Test if code contains a leaf element. If this is the case, the ===
+ // === current element is finished. Traverse the mesh to the next ===
+ // === coarser element. ===
+
+ if (code.isLeafElement()) {
+ int level = elInfo->getLevel();
+
+ do {
+ elInfo = stack.traverseNext(elInfo);
+ } while (elInfo && elInfo->getLevel() > level);
+
+ return false;
+ }
+
+
+ bool meshChanged = false;
+ Element *el = elInfo->getElement();
+
+
+ // === If element is leaf (and the code is not), refine the element. ===
+
+ if (el->isLeaf()) {
+ TEST_EXIT_DBG(elInfo->getLevel() < 255)("This should not happen!\n");
+
+ el->setMark(1);
+ refineManager->setMesh(el->getMesh());
+ refineManager->setStack(&stack);
+ refineManager->refineFunction(elInfo);
+ meshChanged = true;
+ }
+
+
+ // === Continue fitting the mesh structure code to the children of the ===
+ // === current element. ===
+
+ int s0 = el->getSubObjOfChild(0, subObj, ithObj, elInfo->getType());
+ int s1 = el->getSubObjOfChild(1, subObj, ithObj, elInfo->getType());
+ Element *child0 = el->getFirstChild();
+ Element *child1 = el->getSecondChild();
+ if (reverseMode) {
+ swap(s0, s1);
+ swap(child0, child1);
+ }
+
+
+ // === Traverse left child. ===
+
+ if (s0 != -1) {
+ // The edge/face is contained in the left child, therefore fit this
+ // child to the mesh structure code.
+
+ stack.traverseNext(elInfo);
+ code.nextElement();
+ meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
+ elInfo = stack.getElInfo();
+ } else {
+ // The edge/face is not contained in the left child. Hence we need
+ // to traverse through all subelements of the left child until we get
+ // the second child of the current element.
+
+ do {
+ elInfo = stack.traverseNext(elInfo);
+ } while (elInfo && elInfo->getElement() != child1);
+
+ TEST_EXIT_DBG(elInfo != NULL)("This should not happen!\n");
+ }
+
+ TEST_EXIT_DBG(elInfo->getElement() == child1)
+ ("Got wrong child with idx = %d! Searched for child idx = %d\n",
+ elInfo->getElement()->getIndex(), child1->getIndex());
+
+
+ // === Traverse right child. ===
+
+ if (s1 != -1) {
+ // The edge/face is contained in the right child, therefore fit this
+ // child to the mesh structure code.
+
+ code.nextElement();
+ meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
+ } else {
+ // The edge/face is not contained in the right child. Hence we need
+ // to traverse through all subelements of the right child until we have
+ // finished traversing the current element with all its subelements.
+
+ int level = elInfo->getLevel();
+
+ do {
+ elInfo = stack.traverseNext(elInfo);
+ } while (elInfo && elInfo->getLevel() > level);
+ }
+
+
+ return meshChanged;
+ }
}
diff --git a/AMDiS/src/parallel/MeshManipulation.h b/AMDiS/src/parallel/MeshManipulation.h
index 2602f9591a6ef770ff8660a8006f6fffd20c2a7c..cac397290809f5322a293c6a1351136461b3f988 100644
--- a/AMDiS/src/parallel/MeshManipulation.h
+++ b/AMDiS/src/parallel/MeshManipulation.h
@@ -33,17 +33,65 @@ namespace AMDiS {
class MeshManipulation {
public:
- MeshManipulation(FiniteElemSpace *space)
- : feSpace(space),
- mesh(space->getMesh())
- {}
+ MeshManipulation(FiniteElemSpace *space);
+
+ ~MeshManipulation();
void deleteDoubleDofs(std::set<MacroElement*>& newMacroEl, ElementObjects &elObj);
+ /** \brief
+ * Starts the procedure to fit a given edge/face of an element with a mesh
+ * structure code. This functions prepares some data structures and call
+ * then \ref fitElementToMeshCode, that mainly refines the element such that
+ * it fits to the mesh structure code.
+ *
+ * \param[in] code The mesh structure code to which the edge/face of
+ * an element must be fitted.
+ * \param[in] el Pointer to the element.
+ * \param[in] subObj Defines whether an edge or a face must be fitted.
+ * \param[in] ithObj Defines which edge/face must be fitted.
+ * \param[in] elType Element type of the element (only important in 3D).
+ * \param[in] reverseMode Defines, whether the mesh structure code is given
+ * in reverse mode, i.e., left and right children where
+ * changed when the code was created.
+ */
+ bool startFitElementToMeshCode(MeshStructure &code,
+ Element *el,
+ GeoIndex subObj,
+ int ithObj,
+ int elType,
+ bool reverseMode);
+
+ private:
+ /** \brief
+ * Recursively fits a given mesh structure code to an edge/face of an element.
+ * This function is always initialy called from \ref startFitElementToMeshCode.
+ *
+ * \param[in] code The mesh structure code which is used to fit an
+ * edge/face of an element.
+ * \param[in] stack A traverse stack object. The upper most element in this
+ * stack must be used for fitting the mesh structure code
+ * at the current position.
+ * \param[in] subObj Defines whether an edge or a face must be fitted.
+ * \param[in] ithObj Defines which edge/face must be fitted.
+ * \param[in] reverseMode Defines, whether the mesh structure code is given
+ * in reverse mode, i.e., left and right children where
+ * changed when the code was created.
+ */
+ bool fitElementToMeshCode(MeshStructure &code,
+ TraverseStack &stack,
+ GeoIndex subObj,
+ int ithObj,
+ bool reverseMode);
+
+
+
private:
FiniteElemSpace *feSpace;
Mesh *mesh;
+
+ RefinementManager *refineManager;
};
}
diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc
index 4013113c317eb479ca01944cb833a7b6df480317..b312884a8f52da7bc66dde01a2daba50cbe1b0cd 100644
--- a/AMDiS/src/parallel/ParallelDebug.cc
+++ b/AMDiS/src/parallel/ParallelDebug.cc
@@ -143,7 +143,32 @@ namespace AMDiS {
{
FUNCNAME("ParallelDebug::testPeriodicBoundary()");
+ DegreeOfFreedom testDof = pdb.mapGlobalToLocal(968);
+ if (testDof != -1) {
+ WorldVector<double> c;
+ pdb.mesh->getDofIndexCoords(testDof, pdb.feSpace, c);
+
+ MSG("FOUND: LOCAL %d GLOBAL %d COORDS %f %f %f\n",
+ testDof, 968, c[0], c[1], c[2]);
+ }
+
+ testDof = pdb.mapGlobalToLocal(973);
+ if (testDof != -1) {
+ WorldVector<double> c;
+ pdb.mesh->getDofIndexCoords(testDof, pdb.feSpace, c);
+ MSG("FOUND: LOCAL %d GLOBAL %d COORDS %f %f %f\n",
+ testDof, 973, c[0], c[1], c[2]);
+ }
+
+ testDof = pdb.mapGlobalToLocal(1795);
+ if (testDof != -1) {
+ WorldVector<double> c;
+ pdb.mesh->getDofIndexCoords(testDof, pdb.feSpace, c);
+
+ MSG("FOUND: LOCAL %d GLOBAL %d COORDS %f %f %f\n",
+ testDof, 1795, c[0], c[1], c[2]);
+ }
// === 1. check: All periodic DOFs should have at least a correct number ===
// === of periodic associations. ===
diff --git a/AMDiS/src/parallel/PetscSolver.cc b/AMDiS/src/parallel/PetscSolver.cc
index 6a2a23d1a02f9697651c37ffa54897a8821fb9df..f83cb4da9a6472b17916b591f4d9fc9e94027fff 100644
--- a/AMDiS/src/parallel/PetscSolver.cc
+++ b/AMDiS/src/parallel/PetscSolver.cc
@@ -73,6 +73,7 @@ namespace AMDiS {
std::vector<int> globalCols;
+
// === Traverse all rows of the dof matrix and insert row wise the values ===
// === to the petsc matrix. ===
@@ -104,22 +105,28 @@ namespace AMDiS {
cols.push_back(globalColDof * dispMult + dispAddCol);
values.push_back(value(*icursor));
} else {
- std::set<int>& perColAsc = meshDistributor->getPerDofAssociations(globalColDof);
+ std::set<int>& perColAsc =
+ meshDistributor->getPerDofAssociations(globalColDof);
+
std::set<int> perAsc;
- for (std::set<int>::iterator it = perColAsc.begin(); it != perColAsc.end(); ++it)
+ for (std::set<int>::iterator it = perColAsc.begin();
+ it != perColAsc.end(); ++it)
if (*it >= -3)
perAsc.insert(*it);
-
- double scaledValue = value(*icursor) * std::pow(0.5, static_cast<double>(perAsc.size()));
+
+ double scaledValue =
+ value(*icursor) * std::pow(0.5, static_cast<double>(perAsc.size()));
std::vector<int> newCols;
newCols.push_back(globalColDof);
- for (std::set<int>::iterator it = perAsc.begin(); it != perAsc.end(); ++it) {
+ for (std::set<int>::iterator it = perAsc.begin();
+ it != perAsc.end(); ++it) {
int nCols = static_cast<int>(newCols.size());
for (int i = 0; i < nCols; i++) {
- TEST_EXIT(meshDistributor->isPeriodicDof(*it, newCols[i]))
+ TEST_EXIT(meshDistributor->isPeriodicDof(newCols[i], *it))
("Should not happen: %d %d\n", *it, newCols[i]);
- newCols.push_back(meshDistributor->getPeriodicMapping(*it, newCols[i]));
+
+ newCols.push_back(meshDistributor->getPeriodicMapping(newCols[i], *it));
}
}
@@ -139,18 +146,28 @@ namespace AMDiS {
// Global index of the current column index.
int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor));
- std::set<int>& perColAsc = meshDistributor->getPerDofAssociations(globalColDof);
std::set<int> perAsc;
- for (std::set<int>::iterator it = perColAsc.begin(); it != perColAsc.end(); ++it)
- if (*it >= -3)
- perAsc.insert(*it);
- std::set<int>& perRowAsc = meshDistributor->getPerDofAssociations(globalRowDof);
- for (std::set<int>::iterator it = perRowAsc.begin(); it != perRowAsc.end(); ++it)
- if (*it >= -3)
- perAsc.insert(*it);
+ if (meshDistributor->isPeriodicDof(globalColDof)) {
+ std::set<int>& perColAsc =
+ meshDistributor->getPerDofAssociations(globalColDof);
+ for (std::set<int>::iterator it = perColAsc.begin();
+ it != perColAsc.end(); ++it)
+ if (*it >= -3)
+ perAsc.insert(*it);
+ }
+
+ if (meshDistributor->isPeriodicDof(globalRowDof)) {
+ std::set<int>& perRowAsc =
+ meshDistributor->getPerDofAssociations(globalRowDof);
+ for (std::set<int>::iterator it = perRowAsc.begin();
+ it != perRowAsc.end(); ++it)
+ if (*it >= -3)
+ perAsc.insert(*it);
+ }
- double scaledValue = value(*icursor) * std::pow(0.5, static_cast<double>(perAsc.size()));
+ double scaledValue =
+ value(*icursor) * std::pow(0.5, static_cast<double>(perAsc.size()));
std::vector<std::pair<int, int> > entry;
entry.push_back(std::make_pair(globalRowDof, globalColDof));
@@ -159,13 +176,13 @@ namespace AMDiS {
for (int i = 0; i < nEntry; i++) {
int perRowDof = 0;
if (meshDistributor->getPeriodicMapping()[*it].count(entry[i].first))
- perRowDof = meshDistributor->getPeriodicMapping(*it, entry[i].first);
+ perRowDof = meshDistributor->getPeriodicMapping(entry[i].first, *it);
else
perRowDof = entry[i].first;
int perColDof;
if (meshDistributor->getPeriodicMapping()[*it].count(entry[i].second))
- perColDof = meshDistributor->getPeriodicMapping(*it, entry[i].second);
+ perColDof = meshDistributor->getPeriodicMapping(entry[i].second, *it);
else
perColDof = entry[i].second;
@@ -207,7 +224,7 @@ namespace AMDiS {
VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
for (std::set<int>::iterator perIt = perAsc.begin(); perIt != perAsc.end(); ++perIt) {
- int mappedDof = meshDistributor->getPeriodicMapping(*perIt, globalRowDof);
+ int mappedDof = meshDistributor->getPeriodicMapping(globalRowDof, *perIt);
int mappedIndex = mappedDof * dispMult + dispAdd;
VecSetValues(petscVec, 1, &mappedIndex, &value, ADD_VALUES);
}