Commit 77dbdc69 authored by Thomas Witkowski's avatar Thomas Witkowski

Some code refactoring on the last work in parallel AMDiS.

parent 8a2a712f
......@@ -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()");
......
......@@ -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);
......
......@@ -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);