Skip to content

GitLab

Commit fbc5eb23 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

BLUB

parent 32dfc015
Hide whitespace changes
Inline Side-by-side
Showing with 154 additions and 150 deletions
AMDiS/src/MeshStructure.cc
View file @ fbc5eb23
......@@ -62,12 +62,16 @@ namespace AMDiS {
clear();
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
while (elInfo) {
if (macroElIndex == -1 ||
elInfo->getMacroElement()->getIndex() == macroElIndex)
insertElement(elInfo->getElement()->isLeaf());
ElInfo *elInfo;
if (macroElIndex == -1)
elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
else
elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1,
Mesh::CALL_EVERY_EL_PREORDER);
while (elInfo) {
insertElement(elInfo->getElement()->isLeaf());
elInfo = stack.traverseNext(elInfo);
}
......
AMDiS/src/Traverse.h
View file @ fbc5eb23
......@@ -77,25 +77,19 @@ namespace AMDiS {
}
public:
/** \brief
* Returns the first ElInfo of a non recursive traversal which fullfills the
* selection criterion specified by level and fill_flag and initializes the
* stack. After a call of traverseFirst the next Elements are traversed via
* \ref traverseNext().
*/
/// Returns the first ElInfo of a non recursive traversal which fullfills the
/// selection criterion specified by level and fill_flag and initializes the
/// stack. After a call of traverseFirst the next Elements are traversed via
/// \ref traverseNext().
ElInfo* traverseFirst(Mesh *mesh, int level, Flag fill_flag);
/** \brief
* Works in the same way as \ref traverseFirst defined above, but limits the
* traversal to one macro mesh only.
*/
/// Works in the same way as \ref traverseFirst defined above, but limits the
/// traversal to one macro mesh only.
ElInfo* traverseFirstOneMacro(Mesh *mesh, int macroIndex, int level,
Flag fill_flag);
/** \brief
* Returns the next ElInfo in a traversal initiated by \ref traverseFirst()
* If NULL is returned, the traversal is finished.
*/
/// Returns the next ElInfo in a traversal initiated by \ref traverseFirst()
/// If NULL is returned, the traversal is finished.
ElInfo* traverseNext(ElInfo* elinfo_old);
/// Returns the neighbour-th neighbour of elInfoOld
......@@ -176,20 +170,16 @@ namespace AMDiS {
/// Used by \ref traverseFirst() \ref traverseNext()
ElInfo* traverseElementLevel();
/** \brief
* Avoids copy of a traverse stack. If copy should be possible
* the operator must be implemented (deep copy not flat copy!)
*/
/// Avoids copy of a traverse stack. If copy should be possible
/// the operator must be implemented (deep copy not flat copy!)
void operator=(const TraverseStack& /*rhs*/)
{
FUNCNAME("TraverseStack::operator=()");
ERROR_EXIT("not implemented");
}
/** \brief
* Avoids copy of a traverse stack. If copy should be possible
* the operator must be implemented (deep copy not flat copy!)
*/
/// Avoids copy of a traverse stack. If copy should be possible
/// the operator must be implemented (deep copy not flat copy!)
TraverseStack(const TraverseStack&)
{
FUNCNAME("TraverseStack::TraverseStack()");
......@@ -203,16 +193,12 @@ namespace AMDiS {
/// Mesh which is currently traversed
Mesh* traverse_mesh;
/** \brief
* Traverse level. Used if CALL_EL_LEVEL or CALL_LEAF_EL_LEVEL are set in
* \ref traverse_fill_flag
*/
/// Traverse level. Used if CALL_EL_LEVEL or CALL_LEAF_EL_LEVEL are set in
/// \ref traverse_fill_flag
int traverse_level;
/** \brief
* Flags for traversal. Specifies which elements are visited and which
* information are filled to the ElInfo objects.
*/
/// Flags for traversal. Specifies which elements are visited and which
/// information are filled to the ElInfo objects.
Flag traverse_fill_flag;
/// If -1, the whole mesh is traversed. Otherwise, only the macro element
......@@ -231,13 +217,11 @@ namespace AMDiS {
///
std::vector<ElInfo*> elinfo_stack;
/** \brief
* Stores for all level, which children of this element was already
* visited. If 0, no children were visited (this status occurs only
* in intermediate computations). If 1, only the first was vistied.
* If 2, both children of the element on this level were already
* visited.
*/
/// Stores for all level, which children of this element was already
/// visited. If 0, no children were visited (this status occurs only
/// in intermediate computations). If 1, only the first was vistied.
/// If 2, both children of the element on this level were already
/// visited.
std::vector<int> info_stack;
///
......@@ -255,18 +239,14 @@ namespace AMDiS {
///
int id;
/** \brief
* Is used for parallel mesh traverse. The thread with the id
* myThreadId is only allowed to access coarse elements, which id
* satisfies: elId % maxThreads = myThreadId.
*/
/// Is used for parallel mesh traverse. The thread with the id
/// myThreadId is only allowed to access coarse elements, which id
/// satisfies: elId % maxThreads = myThreadId.
int myThreadId;
/** \brief
* Is used for parallel mesh traverse. The thread with the id
* myThreadId is only allowed to access coarse elements, which id
* satisfies: elId % maxThreads = myThreadId.
*/
/// Is used for parallel mesh traverse. The thread with the id
/// myThreadId is only allowed to access coarse elements, which id
/// satisfies: elId % maxThreads = myThreadId.
int maxThreads;
};
......@@ -278,11 +258,9 @@ namespace AMDiS {
class Traverse
{
public:
/** \brief
* Creates a Traverse object for Mesh m with fill flags f and level l.
* el_fct is a pointer to a the function that will be called for every
* visited element.
*/
/// Creates a Traverse object for Mesh m with fill flags f and level l.
/// el_fct is a pointer to a the function that will be called for every
/// visited element.
Traverse(Mesh* m,
const Flag f,
const unsigned char l,
......
AMDiS/src/parallel/InteriorBoundary.cc
View file @ fbc5eb23
......@@ -48,8 +48,6 @@ namespace AMDiS {
for (int geoPos = 0; geoPos < mesh->getDim(); geoPos++) {
GeoIndex geoIndex = INDEX_OF_DIM(geoPos, mesh->getDim());
MSG("---> start with geo index = %d\n", geoIndex);
while (elObjDb.iterate(geoIndex)) {
map<int, ElementObjectData>& objData = elObjDb.getIterateData();
......
AMDiS/src/parallel/MeshDistributor.cc
View file @ fbc5eb23
......@@ -908,6 +908,7 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::checkMeshChange()");
MPI::COMM_WORLD.Barrier();
double first = MPI::Wtime();
int skip = 0;
......@@ -927,24 +928,36 @@ namespace AMDiS {
if (skip == 0) {
int iterationCounter = 0;
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;
for (InteriorBoundary::iterator it(intBoundary.getOwn()); !it.end(); ++it)
// To check the interior boundaries, the ownership of the boundaries is not
// important. Therefore, we add all boundaries to one boundary container.
RankToBoundMap allBound;
for (InteriorBoundary::iterator it(intBoundary.getOwn()); !it.end(); ++it)
if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) ||
(mesh->getDim() == 3 && it->rankObj.subObj == FACE))
allBound[it.getRank()].push_back(*it);
for (InteriorBoundary::iterator it(intBoundary.getOther());
!it.end(); ++it)
for (InteriorBoundary::iterator it(intBoundary.getOther());
!it.end(); ++it)
if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) ||
(mesh->getDim() == 3 && it->rankObj.subObj == FACE))
allBound[it.getRank()].push_back(*it);
for (InteriorBoundary::iterator it(intBoundary.getPeriodic());
!it.end(); ++it)
if (it.getRank() != mpiRank)
if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) ||
(mesh->getDim() == 3 && it->rankObj.subObj == FACE))
allBound[it.getRank()].push_back(*it);
allBound[it.getRank()].push_back(*it);
do {
bool meshChanged = false;
// === Check the boundaries and adapt mesh if necessary. ===
MSG_DBG("Run checkAndAdaptBoundary ...\n");
// Check for periodic boundaries within rank's subdomain.
for (InteriorBoundary::iterator it(intBoundary.getPeriodic());
!it.end(); ++it) {
if (it.getRank() == mpiRank) {
......@@ -956,16 +969,8 @@ namespace AMDiS {
MeshManipulation mm(mesh);
meshChanged |= mm.fitElementToMeshCode(elCode, it->neighObj);
}
} else {
if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) ||
(mesh->getDim() == 3 && it->rankObj.subObj == FACE))
allBound[it.getRank()].push_back(*it);
}
}
// === Check the boundaries and adapt mesh if necessary. ===
MSG_DBG("Run checkAndAdaptBoundary ...\n");
meshChanged |= checkAndAdaptBoundary(allBound);
......@@ -982,11 +987,14 @@ namespace AMDiS {
MSG("Number of iteration to adapt mesh: %d\n", iterationCounter);
}
MPI::COMM_WORLD.Barrier();
INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n",
MPI::Wtime() - first);
#if (DEBUG != 0)
debug::writeMesh(feSpaces[0], -1, debugOutputDir + "mesh");
#endif
// Because the mesh has been changed, update the DOF numbering and mappings.
updateLocalGlobalNumbering();
......@@ -998,16 +1006,12 @@ namespace AMDiS {
ParallelDebug::testPeriodicBoundary(*this);
#endif
// === The mesh has changed, so check if it is required to repartition ===
// === the mesh. ===
nMeshChangesAfterLastRepartitioning++;
INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n",
MPI::Wtime() - first);
if (repartitioningFailed > 0) {
MSG_DBG("Repartitioning not tried because it has failed in the past!\n");
......@@ -1088,7 +1092,7 @@ namespace AMDiS {
FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
// === Create mesh structure codes for all ranks boundary elements. ===
map<int, MeshCodeVec> sendCodes;
for (RankToBoundMap::iterator it = allBound.begin();
......@@ -1109,10 +1113,12 @@ namespace AMDiS {
stdMpi.startCommunication();
// === Compare received mesh structure codes. ===
bool meshChanged = false;
int cCounter = 0;
for (RankToBoundMap::iterator it = allBound.begin();
it != allBound.end(); ++it) {
......@@ -1272,9 +1278,9 @@ namespace AMDiS {
if (repartitioning == 0)
return;
MPI::COMM_WORLD.Barrier();
double timePoint = MPI::Wtime();
#if (DEBUG != 0)
ParallelDebug::testDoubleDofs(mesh);
int writePartMesh = 1;
......@@ -1500,10 +1506,11 @@ namespace AMDiS {
mesh->dofCompress();
partitioner->createPartitionMap(partitionMap);
createInteriorBoundary(false);
updateLocalGlobalNumbering();
// === After mesh adaption, set values of interchange vectors. ===
for (map<int, vector<int> >::iterator it = partitioner->getRecvElements().begin();
......@@ -1521,12 +1528,10 @@ namespace AMDiS {
}
}
// === Update periodic mapping, if there are periodic boundaries. ===
createPeriodicMap();
#if (DEBUG != 0)
MSG("AMDiS runs in debug mode, so make some test ...\n");
......@@ -1540,12 +1545,10 @@ namespace AMDiS {
MSG("Debug mode tests finished!\n");
#endif
// === In 3D we have to make some test, if the resulting mesh is valid. ===
// === If it is not valid, there is no possiblity yet to fix this ===
// === problem, just exit with an error message. ===
// In 3D we have to make some test, if the resulting mesh is valid.
check3dValidMesh();
MPI::COMM_WORLD.Barrier();
MSG("Mesh repartitioning needed %.5f seconds\n", MPI::Wtime() - timePoint);
}
......@@ -1651,6 +1654,9 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::updateLocalGlobalNumbering()");
MPI::COMM_WORLD.Barrier();
double first = MPI::Wtime();
mesh->dofCompress();
#if (DEBUG != 0)
......@@ -1735,6 +1741,9 @@ namespace AMDiS {
ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], dofMap,
debugOutputDir + "mpi-dbg", "dat");
#endif
MPI::COMM_WORLD.Barrier();
MSG("Rebuild of parallel data structures needed %.5f seconds\n", MPI::Wtime() - first);
}
......@@ -1783,8 +1792,15 @@ namespace AMDiS {
TEST_EXIT(levelData.getLevelNumber() == 1)
("Periodic DOF map does not support multi level domain decomposition!\n");
MPI::COMM_WORLD.Barrier();
double first = MPI::Wtime();
for (unsigned int i = 0; i < feSpaces.size(); i++)
createPeriodicMap(feSpaces[i]);
MPI::COMM_WORLD.Barrier();
INFO(info, 8)("Creation of periodic mapping needed %.5f seconds\n",
MPI::Wtime() - first);
}
......
AMDiS/src/parallel/MeshManipulation.cc
View file @ fbc5eb23
......@@ -298,14 +298,16 @@ namespace AMDiS {
// 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(mesh, -1, traverseFlag);
while (elInfo && elInfo->getElement() != boundEl.el)
elInfo = stack.traverseNext(elInfo);
ElInfo *elInfo =
stack.traverseFirstOneMacro(mesh, boundEl.elIndex, -1, traverseFlag);
TEST_EXIT_DBG(elInfo->getElement() == boundEl.el)("This should not happen!\n");
TEST_EXIT_DBG(elInfo->getElement() == boundEl.el)
("This should not happen!\n");
return fitElementToMeshCode(code, stack, boundEl.subObj,
boundEl.ithObj, boundEl.reverseMode);
pcode = &code;
pstack = &stack;
rMode = boundEl.reverseMode;
return fitElementToMeshCode(boundEl.subObj, boundEl.ithObj);
}
......@@ -319,9 +321,8 @@ namespace AMDiS {
// Create traverse stack and traverse the mesh to the element.
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
while (elInfo && elInfo->getElement() != boundEl.el)
elInfo = stack.traverseNext(elInfo);
ElInfo *elInfo =
stack.traverseFirstOneMacro(mesh, boundEl.elIndex, -1, traverseFlag);
TEST_EXIT_DBG(elInfo)("This should not happen!\n");
......@@ -351,14 +352,18 @@ namespace AMDiS {
while (elInfo && elInfo->getElement() != child0)
elInfo = stack.traverseNext(elInfo);
meshChanged |=
fitElementToMeshCode(code, stack, boundEl.subObj, s0, boundEl.reverseMode);
pcode = &code;
pstack = &stack;
rMode = boundEl.reverseMode;
meshChanged |= fitElementToMeshCode(boundEl.subObj, s0);
} else {
while (elInfo && elInfo->getElement() != child1)
elInfo = stack.traverseNext(elInfo);
meshChanged |=
fitElementToMeshCode(code, stack, boundEl.subObj, s1, boundEl.reverseMode);
pcode = &code;
pstack = &stack;
rMode = boundEl.reverseMode;
meshChanged |= fitElementToMeshCode(boundEl.subObj, s1);
}
......@@ -366,18 +371,13 @@ namespace AMDiS {
}
bool MeshManipulation::fitElementToMeshCode(MeshStructure &code,
TraverseStack &stack,
GeoIndex subObj,
int ithObj,
bool reverseMode)
bool MeshManipulation::fitElementToMeshCode(GeoIndex subObj, int ithObj)
{
FUNCNAME("MeshManipulation::fitElementToMeshCode()");
// === Test if there are more elements in stack to check with the code. ===
ElInfo *elInfo = stack.getElInfo();
ElInfo *elInfo = pstack->getElInfo();
if (!elInfo)
return false;
......@@ -386,11 +386,11 @@ namespace AMDiS {
// === current element is finished. Traverse the mesh to the next ===
// === coarser element. ===
if (code.isLeafElement()) {
if (pcode->isLeafElement()) {
int level = elInfo->getLevel();
do {
elInfo = stack.traverseNext(elInfo);
elInfo = pstack->traverseNext(elInfo);
} while (elInfo && elInfo->getLevel() > level);
return false;
......@@ -423,19 +423,18 @@ namespace AMDiS {
// that should be adapted. Thus, we can ommit the refinement.
if (subObj == FACE) {
if (s0 != -1 && s1 == -1 || s0 == -1 && s1 != -1) {
if (ithObj <= 1 && code.lookAhead() == 0) {
if (ithObj <= 1 && pcode->lookAhead() == 0) {
elementRefined = false;
code.nextElement();
stack.traverseNext(elInfo);
pcode->nextElement();
pstack->traverseNext(elInfo);
}
}
}
if (elementRefined) {
MeshStructure t = code;
el->setMark(1);
refineManager->setMesh(el->getMesh());
refineManager->setStack(&stack);
refineManager->setStack(pstack);
refineManager->refineFunction(elInfo);
meshChanged = true;
}
......@@ -448,7 +447,7 @@ namespace AMDiS {
Element *child0 = el->getFirstChild();
Element *child1 = el->getSecondChild();
if (reverseMode) {
if (rMode) {
swap(s0, s1);
swap(child0, child1);
}
......@@ -460,17 +459,17 @@ namespace AMDiS {
// 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();
pstack->traverseNext(elInfo);
pcode->nextElement();
meshChanged |= fitElementToMeshCode(subObj, s0);
elInfo = pstack->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);
elInfo = pstack->traverseNext(elInfo);
} while (elInfo && elInfo->getElement() != child1);
TEST_EXIT_DBG(elInfo != NULL)("This should not happen!\n");
......@@ -487,8 +486,8 @@ namespace AMDiS {
// 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);
pcode->nextElement();
meshChanged |= fitElementToMeshCode(subObj, s1);
} 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
......@@ -497,7 +496,7 @@ namespace AMDiS {
int level = elInfo->getLevel();
do {
elInfo = stack.traverseNext(elInfo);
elInfo = pstack->traverseNext(elInfo);
} while (elInfo && elInfo->getLevel() > level);
}
......
AMDiS/src/parallel/MeshManipulation.h
View file @ fbc5eb23
......@@ -42,10 +42,10 @@ namespace AMDiS {
ElementObjectDatabase &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.
* Starts the procedure to fit a given edge/face of a macro 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.
......@@ -62,25 +62,13 @@ namespace AMDiS {
* 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.
*
* \return Returns true, if the mesh was changed, i.e. refined, to fit the
* element to the structure code.
*/
bool fitElementToMeshCode(MeshStructure &code,
TraverseStack &stack,
GeoIndex subObj,
int ithObj,
bool reverseMode);
bool fitElementToMeshCode(GeoIndex subObj, int ithObj);
......@@ -88,6 +76,12 @@ namespace AMDiS {
Mesh *mesh;
RefinementManager *refineManager;
MeshStructure *pcode;
TraverseStack *pstack;
bool rMode;
};
}
......
AMDiS/src/parallel/PetscSolverFeti.cc
View file @ fbc5eb23
......@@ -423,6 +423,7 @@ namespace AMDiS {
subdomain->setDofMapping(&localDofMap, &primalDofMap);
if (printTimings) {
MPI::COMM_WORLD.Barrier();
timeCounter = MPI::Wtime() - timeCounter;
MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n",
timeCounter);
......@@ -737,9 +738,11 @@ namespace AMDiS {
MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
if (printTimings)