diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index 778a7ffa229d555e46876fe8ba7f83d76507ec1f..30945b2e5c1a28b19c5e7a90354668ee51dfde41 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -49,10 +49,10 @@ namespace AMDiS {
{
public:
/// Type of scalars in the underlying matrix
- typedef double value_type;
+ typedef double value_type;
/// Type of underlying matrix
- typedef mtl::compressed2D<value_type> base_matrix_type;
+ typedef mtl::compressed2D<value_type> base_matrix_type;
/// Type of inserter for the base matrix;
typedef mtl::matrix::inserter<base_matrix_type, mtl::operations::update_plus<value_type> > inserter_type;
diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc
index d7291ed419540c2a217f97d02516d10b407f3924..ded7f6e805e5023b98658210ff1844563dfaed34 100644
--- a/AMDiS/src/Debug.cc
+++ b/AMDiS/src/Debug.cc
@@ -62,11 +62,11 @@ namespace AMDiS {
}
- void writeElementIndexMesh(FiniteElemSpace *feSpace, std::string filename)
+ void writeElementIndexMesh(Mesh *mesh, std::string filename)
{
std::map<int, double> vec;
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
int index = elInfo->getElement()->getIndex();
@@ -74,16 +74,15 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
- ElementFileWriter::writeFile(vec, feSpace, filename);
+ ElementFileWriter::writeFile(vec, mesh, filename);
}
- void highlightElementIndexMesh(FiniteElemSpace *feSpace, int idx,
- std::string filename)
+ void highlightElementIndexMesh(Mesh *mesh, int idx, std::string filename)
{
std::map<int, double> vec;
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
int index = elInfo->getElement()->getIndex();
@@ -91,7 +90,7 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
- ElementFileWriter::writeFile(vec, feSpace, filename);
+ ElementFileWriter::writeFile(vec, mesh, filename);
}
diff --git a/AMDiS/src/Debug.h b/AMDiS/src/Debug.h
index 0c724cd2fed8c401d5eaa2fa54f38badc2eb8bc7..3301f76d0269e4d8c2fb4450f0e4befb2bae76ae 100644
--- a/AMDiS/src/Debug.h
+++ b/AMDiS/src/Debug.h
@@ -72,10 +72,9 @@ namespace AMDiS {
* \param[in] feSpace The FE space to be used.
* \param[in] filename Name of the file.
*/
- void writeElementIndexMesh(FiniteElemSpace *feSpace, std::string filename);
+ void writeElementIndexMesh(Mesh *mesh, std::string filename);
- void highlightElementIndexMesh(FiniteElemSpace *feSpace, int idx,
- std::string filename);
+ void highlightElementIndexMesh(Mesh *mesh, int idx, std::string filename);
void colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Element *el);
diff --git a/AMDiS/src/ElementFileWriter.cc b/AMDiS/src/ElementFileWriter.cc
index 8b6fb3737328ad78f8b55184ee7a9189df251e2f..e50e8608373a62bf5fa27e6de82c07d2163f5481 100644
--- a/AMDiS/src/ElementFileWriter.cc
+++ b/AMDiS/src/ElementFileWriter.cc
@@ -7,7 +7,7 @@
namespace AMDiS {
ElementFileWriter::ElementFileWriter(std::string name_,
- const FiniteElemSpace *feSpace_,
+ Mesh *mesh_,
std::map<int, double> &mapvec)
: name(name_),
tecplotExt(".plt"),
@@ -21,8 +21,7 @@ namespace AMDiS {
indexDecimals(3),
tsModulo(1),
timestepNumber(-1),
- mesh(feSpace_->getMesh()),
- feSpace(feSpace_),
+ mesh(mesh_),
vec(mapvec)
{
if (name != "") {
@@ -93,10 +92,10 @@ namespace AMDiS {
void ElementFileWriter::writeFile(std::map<int, double> &vec,
- const FiniteElemSpace *feSpace,
+ Mesh *mesh,
std::string filename)
{
- ElementFileWriter efw("", feSpace, vec);
+ ElementFileWriter efw("", mesh, vec);
efw.writeVtkValues(filename);
}
diff --git a/AMDiS/src/ElementFileWriter.h b/AMDiS/src/ElementFileWriter.h
index 74e16582a6a507d8a95158483969bba22ac5b82b..2b3a4164a3b3c06429fc19faeb3b72260b7ad024 100644
--- a/AMDiS/src/ElementFileWriter.h
+++ b/AMDiS/src/ElementFileWriter.h
@@ -22,7 +22,7 @@ namespace AMDiS {
public:
/// Constructor.
ElementFileWriter(std::string name,
- const FiniteElemSpace *feSpace,
+ Mesh *mesh,
std::map<int, double> &vec);
/// Implementation of FileWriterInterface::writeFiles().
@@ -33,7 +33,7 @@ namespace AMDiS {
/// Simple writing procedure for one element map.
static void writeFile(std::map<int, double> &vec,
- const FiniteElemSpace *feSpace,
+ Mesh *mesh,
std::string filename);
protected:
@@ -101,9 +101,6 @@ namespace AMDiS {
/// Mesh used for output.
Mesh *mesh;
- /// fespace used for output.
- const FiniteElemSpace *feSpace;
-
/// Vector that stores the solution.
std::map<int, double> vec;
};
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index 4a0a2edd87c05b2f1cbdaaee873544b6f0d7fd28..93b33f3f066fd4c7a96a9597dd60ef98832643b2 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/AMDiS/src/ProblemScal.cc
@@ -620,6 +620,7 @@ namespace AMDiS {
createPrecon();
}
+
void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, std::string name)
{
FUNCNAME("ProblemScal::writeResidualMesh()");
@@ -634,10 +635,11 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
- ElementFileWriter fw(name, this->getFeSpace(), vec);
+ ElementFileWriter fw(name, this->getMesh(), vec);
fw.writeFiles(adaptInfo, true);
}
+
void ProblemScal::createPrecon()
{
std::string preconType("no");
@@ -646,15 +648,16 @@ namespace AMDiS {
CreatorInterface<ITL_BasePreconditioner> *preconCreator =
CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
- solver->setLeftPrecon( preconCreator->create(systemMatrix->getBaseMatrix()) );
+ solver->setLeftPrecon(preconCreator->create(systemMatrix->getBaseMatrix()));
preconType= "no";
GET_PARAMETER(0, name + "->solver->right precon", &preconType);
preconCreator = CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
- solver->setRightPrecon( preconCreator->create(systemMatrix->getBaseMatrix()) );
+ solver->setRightPrecon(preconCreator->create(systemMatrix->getBaseMatrix()));
}
+
void ProblemScal::serialize(std::ostream &out)
{
FUNCNAME("ProblemScal::serialize()");
@@ -663,6 +666,7 @@ namespace AMDiS {
solution->serialize(out);
}
+
void ProblemScal::deserialize(std::istream &in)
{
FUNCNAME("ProblemScal::deserialize()");
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index da190afef927db9b9e417785b011b5f824908549..1f9b8f49777d14024aec56cbc3645dedbc5b5e49 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -1538,7 +1538,7 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
- ElementFileWriter::writeFile(vec, this->getFeSpace(comp), name);
+ ElementFileWriter::writeFile(vec, this->getMesh(comp), name);
}
diff --git a/AMDiS/src/RCNeighbourList.cc b/AMDiS/src/RCNeighbourList.cc
index 862a6df11fb82b575c112148afb31433e111ced3..a443191a6e4d0a809a4bc95a1469d6d272846302 100644
--- a/AMDiS/src/RCNeighbourList.cc
+++ b/AMDiS/src/RCNeighbourList.cc
@@ -259,6 +259,7 @@ namespace AMDiS {
}
}
+
void RCNeighbourList::removeDOFParent(int index)
{
Tetrahedron *el = dynamic_cast<Tetrahedron*>(rclist[index]->el);
@@ -297,6 +298,7 @@ namespace AMDiS {
}
}
+
RCNeighbourList *RCNeighbourList::periodicSplit(DegreeOfFreedom *edge[2],
DegreeOfFreedom *nextEdge[2],
int *n_neigh,
diff --git a/AMDiS/src/RefinementManager2d.cc b/AMDiS/src/RefinementManager2d.cc
index c2dc255ab2bc8e4983b833ee89b70b8107472c3a..59a3ddefffc062b75482a27680bde2b6500a0c92 100644
--- a/AMDiS/src/RefinementManager2d.cc
+++ b/AMDiS/src/RefinementManager2d.cc
@@ -12,10 +12,11 @@
#include "Projection.h"
#include "LeafData.h"
#include "VertexVector.h"
+#include "Debug.h"
namespace AMDiS {
- ElInfo* RefinementManager2d::refineFunction(ElInfo* el_info)
+ ElInfo* RefinementManager2d::refineFunction(ElInfo* elInfo)
{
FUNCNAME("RefinementManager::refineFunction()");
@@ -23,41 +24,38 @@ namespace AMDiS {
DegreeOfFreedom *edge[2];
RCNeighbourList* refineList = new RCNeighbourList(2);
- if (el_info->getElement()->getMark() <= 0) {
+ if (elInfo->getElement()->getMark() <= 0) {
delete refineList;
// Element may not be refined.
- return el_info;
+ return elInfo;
}
- refineList->setElement(0, el_info->getElement());
+ refineList->setElement(0, elInfo->getElement());
int n_neigh = 1;
- if (el_info->getProjection(0) &&
- el_info->getProjection(0)->getType() == VOLUME_PROJECTION)
+ if (elInfo->getProjection(0) &&
+ elInfo->getProjection(0)->getType() == VOLUME_PROJECTION)
newCoords = true;
// === Give the refinement edge the right orientation. ===
- if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1, 0)) {
- edge[0] = const_cast<int*>(el_info->getElement()->getDOF(0));
- edge[1] = const_cast<int*>(el_info->getElement()->getDOF(1));
+ if (elInfo->getElement()->getDOF(0, 0) < elInfo->getElement()->getDOF(1, 0)) {
+ edge[0] = const_cast<int*>(elInfo->getElement()->getDOF(0));
+ edge[1] = const_cast<int*>(elInfo->getElement()->getDOF(1));
} else {
- edge[1] = const_cast<int*>(el_info->getElement()->getDOF(0));
- edge[0] = const_cast<int*>(el_info->getElement()->getDOF(1));
+ edge[1] = const_cast<int*>(elInfo->getElement()->getDOF(0));
+ edge[0] = const_cast<int*>(elInfo->getElement()->getDOF(1));
}
// === Get the refinement patch. ===
- if (getRefinePatch(&el_info, edge, 0, refineList, &n_neigh)) {
- // Domain's boundary was reached while looping around the refinement edge.
- getRefinePatch(&el_info, edge, 1, refineList, &n_neigh);
- bound = true;
- }
+ getRefinePatch(&elInfo, edge, 0, refineList, &n_neigh);
+
// === Check for periodic boundary ===
- DegreeOfFreedom *next_edge[2];
+ DegreeOfFreedom *next_edge[2] = {NULL, NULL};
DegreeOfFreedom *first_edge[2] = {edge[0], edge[1]};
DegreeOfFreedom *last_edge[2] = {NULL, NULL};
int n_neigh_periodic;
@@ -121,25 +119,25 @@ namespace AMDiS {
delete refineList;
- return el_info;
+ return elInfo;
}
- void RefinementManager2d::newCoordsFct(ElInfo *el_info)
+ void RefinementManager2d::newCoordsFct(ElInfo *elInfo)
{
- Element *el = el_info->getElement();
+ Element *el = elInfo->getElement();
int dow = Global::getGeo(WORLD);
- Projection *projector = el_info->getProjection(0);
+ Projection *projector = elInfo->getProjection(0);
if (!projector || projector->getType() != VOLUME_PROJECTION)
- projector = el_info->getProjection(2);
+ projector = elInfo->getProjection(2);
if (el->getFirstChild() && projector && (!el->isNewCoordSet())) {
WorldVector<double> *new_coord = new WorldVector<double>;
for (int j = 0; j < dow; j++)
- (*new_coord)[j] = (el_info->getCoord(0)[j] + el_info->getCoord(1)[j]) * 0.5;
+ (*new_coord)[j] = (elInfo->getCoord(0)[j] + elInfo->getCoord(1)[j]) * 0.5;
projector->project(*new_coord);
el->setNewCoord(new_coord);
@@ -235,10 +233,10 @@ namespace AMDiS {
void RefinementManager2d::bisectTriangle(Triangle *el, DegreeOfFreedom* newDOFs[3])
{
FUNCNAME("RefinementManager2d::bisectTriangle()");
+
TEST_EXIT_DBG(mesh)("no mesh!\n");
Triangle *child[2];
-
child[0] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
child[1] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
@@ -302,36 +300,47 @@ namespace AMDiS {
}
- int RefinementManager2d::getRefinePatch(ElInfo **el_info,
+ void RefinementManager2d::getRefinePatch(ElInfo **elInfo,
DegreeOfFreedom *edge[2],
int dir, RCNeighbourList* refineList,
int *n_neigh)
{
FUNCNAME("RefinementManager2d::getRefinePatch()");
- if ((*el_info)->getNeighbour(2) && (*el_info)->getOppVertex(2) != 2) {
+ if ((*elInfo)->getNeighbour(2) && (*elInfo)->getOppVertex(2) != 2) {
// Neighbour is not compatible devisible; refine neighbour first, store the
// opp_vertex to traverse back to el.
- int opp_vertex = (*el_info)->getOppVertex(2);
+ int opp_vertex = (*elInfo)->getOppVertex(2);
- ElInfo *neigh_info = stack->traverseNeighbour2d(*el_info, 2);
+ ElInfo *neigh_info = stack->traverseNeighbour2d(*elInfo, 2);
neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(), 1));
neigh_info = refineFunction(neigh_info);
-
- // Now go back to the original element and refine the patch.
- *el_info = neigh_info = stack->traverseNeighbour2d(neigh_info, opp_vertex);
+
+ // === Now go back to the original element and refine the patch. ===
+
+
+ // In Debug mode we test if traverNeighbour2d returns the expected element.
+#if (DEBUG != 0)
+ int testIndex = neigh_info->getNeighbour(opp_vertex)->getIndex();
+#endif
+
+ *elInfo = neigh_info = stack->traverseNeighbour2d(neigh_info, opp_vertex);
+
+
+ TEST_EXIT_DBG(testIndex == (*elInfo)->getElement()->getIndex())
+ ("Got wrong neighbour! Should be %d, but is %d!\n",
+ testIndex, (*elInfo)->getElement()->getIndex());
+
TEST_EXIT_DBG(neigh_info->getElement() ==
- dynamic_cast<Triangle*>(const_cast<Element*>((*el_info)->getElement())))
+ dynamic_cast<Triangle*>(const_cast<Element*>((*elInfo)->getElement())))
("invalid traverse_neighbour1\n");
}
-
- if (refineList->setElement(1, (*el_info)->getNeighbour(2))) {
- TEST_EXIT_DBG((*el_info)->getOppVertex(2) == 2)
+
+ if (refineList->setElement(1, (*elInfo)->getNeighbour(2))) {
+ TEST_EXIT_DBG((*elInfo)->getOppVertex(2) == 2)
("no compatible ref. edge after recursive refinement of neighbour\n");
*n_neigh = 2;
}
-
- return 0;
}
}
diff --git a/AMDiS/src/RefinementManager2d.h b/AMDiS/src/RefinementManager2d.h
index 3e535eb99e33ae9794cc1c07a1c9f1224f60af2a..ec1e07551abeb0dae07c0ceff6e7921bebebcdf0 100644
--- a/AMDiS/src/RefinementManager2d.h
+++ b/AMDiS/src/RefinementManager2d.h
@@ -54,15 +54,16 @@ namespace AMDiS {
* get_refine_patch returns 1 if the domain's boundary is reached while
* looping around the refinement edge, otherwise 0
*/
- int getRefinePatch(ElInfo **el_info, DegreeOfFreedom *edge[2], int dir,
- RCNeighbourList* refineList, int *n_neigh);
+ void getRefinePatch(ElInfo **elInfo, DegreeOfFreedom *edge[2], int dir,
+ RCNeighbourList* refineList, int *n_neigh);
/// Refines all elements in the patch.
- DegreeOfFreedom refinePatch(DegreeOfFreedom *edge[2], RCNeighbourList* refineList,
+ DegreeOfFreedom refinePatch(DegreeOfFreedom *edge[2],
+ RCNeighbourList* refineList,
int n_neigh, bool bound);
/// Implements RefinementManager::refineFunction.
- ElInfo* refineFunction(ElInfo* el_info);
+ ElInfo* refineFunction(ElInfo* elInfo);
/// Refines one Triangle.
void bisectTriangle(Triangle *el, DegreeOfFreedom* newDofs[3]);
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index 07fb4b5a703860a8e479ad1996fc877185c0f571..99552c989e42b1fd866789fefd0b1c3b570b03b1 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -135,9 +135,8 @@ namespace AMDiS {
currentMacro = traverse_mesh->firstMacroElement();
while (((*currentMacro)->getIndex() % maxThreads != myThreadId) &&
- (currentMacro != traverse_mesh->endOfMacroElements())) {
- currentMacro++;
- }
+ currentMacro != traverse_mesh->endOfMacroElements())
+ currentMacro++;
if (currentMacro == traverse_mesh->endOfMacroElements())
return NULL;
@@ -349,13 +348,8 @@ namespace AMDiS {
int fillIthChild = info_stack[stack_used];
info_stack[stack_used]++;
- if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) {
- if (fillIthChild == 0)
- fillIthChild = 1;
- else
- fillIthChild = 0;
- }
-
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
+ fillIthChild = 1 - fillIthChild;
elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
stack_used++;
@@ -450,7 +444,8 @@ namespace AMDiS {
case 3:
return traverseNeighbour3d(elinfo_old, neighbour);
break;
- default: ERROR_EXIT("invalid dim\n");
+ default:
+ ERROR_EXIT("invalid dim\n");
}
return NULL;
}
@@ -461,11 +456,12 @@ namespace AMDiS {
FUNCNAME("TraverseStackLLtraverseNeighbour3d()");
Element *el2;
- ElInfo *old_elinfo, *elinfo, *elinfo2;
+ ElInfo *elinfo2;
const DegreeOfFreedom *dof;
- int i, nb, opp_vertex, stack2_used, typ = 0;
- int sav_index, sav_neighbour = neighbour;
+ int stack2_used;
+ int sav_neighbour = neighbour;
+ // father.neigh[coarse_nb[i][j]] == child[i-1].neigh[j]
static int coarse_nb[3][3][4] = {{{-2,-2,-2,-2}, {-1,2,3,1}, {-1,3,2,0}},
{{-2,-2,-2,-2}, {-1,2,3,1}, {-1,2,3,0}},
{{-2,-2,-2,-2}, {-1,2,3,1}, {-1,2,3,0}}};
@@ -473,26 +469,29 @@ namespace AMDiS {
TEST_EXIT_DBG(stack_used > 0)("no current element\n");
Parametric *parametric = traverse_mesh->getParametric();
-
if (parametric)
elinfo_old = parametric->removeParametricInfo(elinfo_old);
- TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])
- ("invalid old elinfo\n");
-
+ TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo\n");
TEST_EXIT_DBG(elinfo_stack[stack_used]->getFillFlag().isSet(Mesh::FILL_NEIGH))
("FILL_NEIGH not set");
+
Element *el = elinfo_stack[stack_used]->getElement();
- sav_index = el->getIndex();
+ int sav_index = el->getIndex();
/* first, goto to leaf level, if necessary... */
if ((traverse_fill_flag & Mesh::CALL_LEAF_EL).isAnySet()) {
- if ((el->getChild(0)) && (neighbour < 2)) {
+ if (el->getChild(0) && neighbour < 2) {
if (stack_used >= stack_size - 1)
enlargeTraverseStack();
- i = 1 - neighbour;
+ int i = 1 - neighbour;
+
elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
- info_stack[stack_used] = i + 1;
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
+ info_stack[stack_used] = (i == 0 ? 2 : 1);
+ else
+ info_stack[stack_used] = i + 1;
+
stack_used++;
info_stack[stack_used] = 0;
neighbour = 3;
@@ -502,30 +501,43 @@ namespace AMDiS {
/* save information about current element and its position in the tree */
save_traverse_mel = traverse_mel;
save_stack_used = stack_used;
+ for (int i = stack_used; i <= save_stack_used; i++) {
+ save_info_stack[i] = info_stack[i];
+ *(save_elinfo_stack[i]) = *(elinfo_stack[i]);
+ }
+ ElInfo *old_elinfo = save_elinfo_stack[save_stack_used];
+ int opp_vertex = old_elinfo->getOppVertex(neighbour);
+
- nb = neighbour;
+ // === First phase (see 2D). ===
+
+ int nb = neighbour;
while (stack_used > 1) { /* go up in tree until we can go down again */
stack_used--;
- typ = elinfo_stack[stack_used]->getType();
- nb = coarse_nb[typ][info_stack[stack_used]][nb];
+ int typ = elinfo_stack[stack_used]->getType();
+ int elIsIthChild = info_stack[stack_used];
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE) && elIsIthChild != 0)
+ elIsIthChild = (elIsIthChild == 1 ? 2 : 1);
+
+ TEST_EXIT_DBG(!elinfo_stack[stack_used + 1]->getParent() ||
+ elinfo_stack[stack_used + 1]->getParent()->getChild(elIsIthChild - 1) ==
+ elinfo_stack[stack_used + 1]->getElement())
+ ("Should not happen!\n");
+
+
+ nb = coarse_nb[typ][elIsIthChild][nb];
if (nb == -1)
break;
+
TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
}
- /* save hierarchy information about current element */
-
- for (i = stack_used; i <= save_stack_used; i++) {
- save_info_stack[i] = info_stack[i];
- *(save_elinfo_stack[i]) = *(elinfo_stack[i]);
- }
- old_elinfo = save_elinfo_stack[save_stack_used];
- opp_vertex = old_elinfo->getOppVertex(neighbour);
+ if (nb >= 0) {
+ // Go to macro element neighbour.
- if (nb >= 0) { /* go to macro element neighbour */
- i = traverse_mel->getOppVertex(nb);
+ int i = traverse_mel->getOppVertex(nb);
traverse_mel = traverse_mel->getNeighbour(nb);
if (traverse_mel == NULL)
return NULL;
@@ -541,7 +553,9 @@ namespace AMDiS {
elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
info_stack[stack_used] = 0;
nb = i;
- } else { /* goto other child */
+ } else {
+ // Goto other child.
+
stack2_used = stack_used + 1;
if (save_stack_used > stack2_used)
stack2_used++; /* go down one level in OLD hierarchy */
@@ -551,24 +565,37 @@ namespace AMDiS {
if (stack_used >= stack_size-1)
enlargeTraverseStack();
- i = 2 - info_stack[stack_used];
- info_stack[stack_used] = i+1;
- elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
+
+ int i = 2 - info_stack[stack_used];
+ info_stack[stack_used] = i + 1;
+ int fillIthChild = i;
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
+ fillIthChild = 1 - fillIthChild;
+ elinfo_stack[stack_used+1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
stack_used++;
info_stack[stack_used] = 0;
nb = 0;
}
- elinfo = elinfo_stack[stack_used];
+
+ // === Second phase. ===
+
+ ElInfo *elinfo = elinfo_stack[stack_used];
el = elinfo->getElement();
- while(el->getChild(0)) {
-
- if (nb < 2) { /* go down one level in hierarchy */
- if (stack_used >= stack_size-1)
+ while (el->getChild(0)) {
+ if (nb < 2) {
+ // Go down one level in hierarchy.
+
+ if (stack_used >= stack_size - 1)
enlargeTraverseStack();
- info_stack[stack_used] = 2-nb;
- elinfo_stack[stack_used+1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
+ info_stack[stack_used] = 2 - nb;
+
+ int fillIthChild = 1 - nb;
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
+ fillIthChild = 1 - fillIthChild;
+ elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
+
stack_used++;
info_stack[stack_used] = 0;
elinfo = elinfo_stack[stack_used];
@@ -576,9 +603,12 @@ namespace AMDiS {
nb = 3;
}
- if (save_stack_used > stack2_used) { /* `refine' both el and el2 */
+ if (save_stack_used > stack2_used) {
+ // `refine' both el and el2.
+
TEST_EXIT_DBG(el->getChild(0))("invalid new refinement?\n");
+ int i = 0;
if (el->getDOF(0) == el2->getDOF(0))
i = save_info_stack[stack2_used] - 1;
else if (el->getDOF(1) == el2->getDOF(0))
@@ -593,20 +623,22 @@ namespace AMDiS {
}
}
- if (el->getChild(0) &&
+ if (el->getChild(0) &&
(el->getChild(i)->getDOF(1) == el->getDOF(nb) ||
- traverse_mesh->associated(el->getChild(i)->getDOF(1, 0), el->getDOF(nb, 0))))
- { /* el->child[0] Kuni 22.08.96 */
- nb = 1;
- }
- else {
- nb = 2;
- }
+ traverse_mesh->associated(el->getChild(i)->getDOF(1, 0), el->getDOF(nb, 0))))
+ nb = 1;
+ else
+ nb = 2;
- info_stack[stack_used] = i+1;
- if (stack_used >= stack_size-1)
+ info_stack[stack_used] = i + 1;
+ if (stack_used >= stack_size - 1)
enlargeTraverseStack();
- elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
+
+ int fillIthChild = i;
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
+ fillIthChild = 1 - fillIthChild;
+ elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
+
stack_used++;
info_stack[stack_used] = 0;
@@ -627,9 +659,9 @@ namespace AMDiS {
el2 = elinfo2->getElement();
}
}
+ } else {
+ // Now we're done...
- }
- else { /* now we're done... */
elinfo = elinfo_stack[stack_used];
el = elinfo->getElement();
break;
@@ -639,11 +671,11 @@ namespace AMDiS {
if (elinfo->getNeighbour(opp_vertex) != old_elinfo->getElement()) {
MSG(" looking for neighbour %d of element %d at %p\n",
- neighbour, old_elinfo->getElement()->getIndex(), reinterpret_cast<void*>( old_elinfo->getElement()));
+ neighbour, old_elinfo->getElement()->getIndex(), reinterpret_cast<void*>(old_elinfo->getElement()));
MSG(" originally: neighbour %d of element %d at %p\n",
- sav_neighbour, sav_index, reinterpret_cast<void*>( old_elinfo->getElement()));
+ sav_neighbour, sav_index, reinterpret_cast<void*>(old_elinfo->getElement()));
MSG(" got element %d at %p with opp_vertex %d neigh %d\n",
- elinfo->getElement()->getIndex(), reinterpret_cast<void*>( elinfo->getElement()),
+ elinfo->getElement()->getIndex(), reinterpret_cast<void*>(elinfo->getElement()),
opp_vertex, (elinfo->getNeighbour(opp_vertex))?(elinfo->getNeighbour(opp_vertex))->getIndex():-1);
TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
("didn't succeed !?!?!?\n");
@@ -669,49 +701,56 @@ namespace AMDiS {
{
FUNCNAME("TraverseStack::traverseNeighbour2d()");
- Triangle *el, *el2, *sav_el;
- ElInfo *old_elinfo, *elinfo, *elinfo2;
- int i, nb, opp_vertex, stack2_used;
+ Triangle *el2;
+ ElInfo *elinfo2;
+ int stack2_used;
+ int sav_neighbour = neighbour;
+ // father.neigh[coarse_nb[i][j]] == child[i-1].neigh[j]
static int coarse_nb[3][3] = {{-2,-2,-2}, {2,-1,1}, {-1,2,0}};
- /* father.neigh[coarse_nb[i][j]] == child[i-1].neigh[j] */
-
- int sav_index, sav_neighbour = neighbour;
TEST_EXIT_DBG(stack_used > 0)("no current element");
- TEST_EXIT_DBG(traverse_fill_flag.isSet(Mesh::CALL_LEAF_EL))("invalid traverse_fill_flag");
Parametric *parametric = traverse_mesh->getParametric();
-
if (parametric)
elinfo_old = parametric->removeParametricInfo(elinfo_old);
TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo");
elinfo_stack[stack_used]->testFlag(Mesh::FILL_NEIGH);
- el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo_stack[stack_used]->getElement()));
- sav_index = el->getIndex();
- sav_el = el;
+ Triangle *el =
+ dynamic_cast<Triangle*>(const_cast<Element*>(elinfo_stack[stack_used]->getElement()));
+ int sav_index = el->getIndex();
+ Triangle *sav_el = el;
/* first, goto to leaf level, if necessary... */
- if (!(el->isLeaf()) && (neighbour < 2)) {
- if (stack_used >= static_cast<int>( elinfo_stack.size())-1) enlargeTraverseStack();
- i = 1 - neighbour;
- elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
- info_stack[stack_used] = i + 1;
+ if (!(el->isLeaf()) && neighbour < 2) {
+ if (stack_used >= static_cast<int>(elinfo_stack.size()) - 1)
+ enlargeTraverseStack();
+
+ // If we should search for neighbour 0, take second child, if for
+ // neighbour 1, take the first child.
+ int i = 1 - neighbour;
+
+ elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
+ info_stack[stack_used] = (i == 0 ? 2 : 1);
+ else
+ info_stack[stack_used] = i + 1;
+
stack_used++;
- neighbour = 2;
+ neighbour = 2;
}
/* save information about current element and its position in the tree */
save_traverse_mel = traverse_mel;
- save_stack_used = stack_used;
- for (i=0; i<=stack_used; i++)
+ save_stack_used = stack_used;
+ for (int i = 0; i <= stack_used; i++) {
save_info_stack[i] = info_stack[i];
- for (i=0; i<=stack_used; i++)
(*(save_elinfo_stack[i])) = (*(elinfo_stack[i]));
- old_elinfo = save_elinfo_stack[stack_used];
- opp_vertex = old_elinfo->getOppVertex(neighbour);
+ }
+ ElInfo *old_elinfo = save_elinfo_stack[stack_used];
+ int opp_vertex = old_elinfo->getOppVertex(neighbour);
/****************************************************************************/
@@ -721,58 +760,72 @@ namespace AMDiS {
/* element of the OLD hierarchy branch to the new branch */
/****************************************************************************/
- nb = neighbour;
+ int nb = neighbour;
while (stack_used > 1) {
stack_used--;
- nb = coarse_nb[info_stack[stack_used]][nb];
+ int elIsIthChild = info_stack[stack_used];
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE) && elIsIthChild != 0)
+ elIsIthChild = (elIsIthChild == 1 ? 2 : 1);
+
+ TEST_EXIT_DBG(!elinfo_stack[stack_used + 1]->getParent() ||
+ elinfo_stack[stack_used + 1]->getParent()->getChild(elIsIthChild - 1) ==
+ elinfo_stack[stack_used + 1]->getElement())
+ ("Should not happen!\n");
+
+ nb = coarse_nb[elIsIthChild][nb];
if (nb == -1)
break;
TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb);
}
+
/****************************************************************************/
/* Now, goto neighbouring element at the local hierarchy entry */
/* This is either a macro element neighbour or the other child of parent. */
/* initialize nb for second phase (see below) */
/****************************************************************************/
- if (nb >= 0) { /* go to macro element neighbour */
+ if (nb >= 0) {
+ // Go to macro element neighbour.
- if ((nb < 2) && (save_stack_used > 1)) {
+ if (nb < 2 && save_stack_used > 1)
stack2_used = 2; /* go down one level in OLD hierarchy */
- } else {
+ else
stack2_used = 1;
- }
+
elinfo2 = save_elinfo_stack[stack2_used];
- el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
+ el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
- i = traverse_mel->getOppVertex(nb);
+ int i = traverse_mel->getOppVertex(nb);
traverse_mel = traverse_mel->getNeighbour(nb);
if (traverse_mel == NULL)
return NULL;
nb = i;
stack_used = 1;
- elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>( traverse_mel));
+ elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>(traverse_mel));
info_stack[stack_used] = 0;
- } else { /* goto other child */
+ } else {
+ // Go to other child.
stack2_used = stack_used + 1;
if (save_stack_used > stack2_used)
stack2_used++; /* go down one level in OLD hierarchy */
elinfo2 = save_elinfo_stack[stack2_used];
- el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
-
+ el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
if (stack_used >= stack_size - 1)
enlargeTraverseStack();
- i = 2 - info_stack[stack_used];
+ int i = 2 - info_stack[stack_used];
info_stack[stack_used] = i + 1;
- elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
+ int fillIthChild = i;
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
+ fillIthChild = 1 - fillIthChild;
+ elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
stack_used++;
nb = 1 - i;
}
@@ -783,45 +836,58 @@ namespace AMDiS {
/* new hierarchy branch to the OLD branch. */
/****************************************************************************/
- elinfo = elinfo_stack[stack_used];
- el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
+ ElInfo *elinfo = elinfo_stack[stack_used];
+ el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement()));
while (el->getFirstChild()) {
+ if (nb < 2) {
+ // Go down one level in hierarchy.
- if (nb < 2) { /* go down one level in hierarchy */
- if (stack_used >= stack_size-1)
+ if (stack_used >= stack_size - 1)
enlargeTraverseStack();
- elinfo_stack[stack_used + 1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
+
+ int fillIthChild = 1 - nb;
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
+ fillIthChild = 1 - fillIthChild;
+ elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
+
info_stack[stack_used] = 2 - nb;
stack_used++;
nb = 2;
}
- if (save_stack_used > stack2_used) { /* `refine' both el and el2 */
+ if (save_stack_used > stack2_used) {
+ // `refine' both el and el2.
+
TEST_EXIT_DBG(el->getFirstChild())("invalid new refinement?");
- /* use child i, neighbour of el2->child[nb-1] */
- i = 2 - save_info_stack[stack2_used];
+ // Use child i, neighbour of el2->child[nb - 1].
+ int i = 2 - save_info_stack[stack2_used];
TEST_EXIT_DBG(i < 2)("invalid OLD refinement?");
- info_stack[stack_used] = i+1;
- elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
+ info_stack[stack_used] = i + 1;
+
+ int fillIthChild = i;
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
+ fillIthChild = 1 - fillIthChild;
+ elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
+
stack_used++;
nb = i;
elinfo = elinfo_stack[stack_used];
- el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
+ el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement()));
stack2_used++;
- if (save_stack_used > stack2_used) {
+ if (save_stack_used > stack2_used)
stack2_used++; /* go down one level in OLD hierarchy */
- }
+
elinfo2 = save_elinfo_stack[stack2_used];
- el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
+ el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
+ } else {
+ // Now we're done...
- }
- else { /* now we're done... */
elinfo = elinfo_stack[stack_used];
- el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
+ el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement()));
}
}
@@ -832,7 +898,7 @@ namespace AMDiS {
sav_neighbour, sav_index, sav_el);
MSG(" got element %d at %8X with opp_vertex %d neigh %d\n",
elinfo->getElement()->getIndex(), elinfo->getElement(), opp_vertex,
- elinfo->getNeighbour(opp_vertex)?elinfo->getNeighbour(opp_vertex)->getIndex():-1);
+ elinfo->getNeighbour(opp_vertex) ? elinfo->getNeighbour(opp_vertex)->getIndex() : -1);
TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
("didn't succeed !?!?!?");
}
@@ -845,8 +911,7 @@ namespace AMDiS {
MSG(" got element %d at %8X with opp_vertex %d neigh %d\n",
elinfo->getElement()->getIndex(), elinfo->getElement(), opp_vertex,
elinfo->getNeighbour(opp_vertex)->getIndex());
- MSG("got no leaf element\n");
- WAIT_REALLY;
+ ERROR_EXIT("got no leaf element\n");
}
if (elinfo) {
diff --git a/AMDiS/src/Triangle.cc b/AMDiS/src/Triangle.cc
index e66d46cfb6e096692db62644b8e2ee64626ac4af..703b34589799d8e24425c59ea3ba50a2a9502ae3 100644
--- a/AMDiS/src/Triangle.cc
+++ b/AMDiS/src/Triangle.cc
@@ -81,6 +81,13 @@ namespace AMDiS {
{
FUNCNAME("Triangle::getVertexDofs()");
+ if (bound.subObj == VERTEX) {
+ dofs.push_back(dof[bound.ithObj]);
+ return;
+ }
+
+ TEST_EXIT_DBG(bound.subObj == EDGE)("This should not happen!\n");
+
BoundaryObject nextBound = bound;
switch (bound.ithObj) {
diff --git a/AMDiS/src/parallel/ParallelDomainBase.cc b/AMDiS/src/parallel/ParallelDomainBase.cc
index 9ba19d5d0b2c5434db1dab7075a13342b8cce737..5081cc1085422d7ce6c75fe617768375987b2b23 100644
--- a/AMDiS/src/parallel/ParallelDomainBase.cc
+++ b/AMDiS/src/parallel/ParallelDomainBase.cc
@@ -73,9 +73,12 @@ namespace AMDiS {
TEST_EXIT(mesh)("No mesh has been defined for mesh distribution!\n");
- // If the problem has been already read from a file, we do not need to do anything.
- if (deserialized)
+ // If the problem has been already read from a file, we need only to remove
+ // periodic boundary conditions (if there are some).
+ if (deserialized) {
+ removePeriodicBoundaryConditions();
return;
+ }
// Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
@@ -192,32 +195,7 @@ namespace AMDiS {
// === Remove periodic boundary conditions in sequential problem definition. ===
- // Remove periodic boundaries in boundary manager on matrices and vectors.
- for (unsigned int i = 0; i < probStat.size(); i++) {
- int nComponents = probStat[i]->getNumComponents();
-
- for (int j = 0; j < nComponents; j++) {
- for (int k = 0; k < nComponents; k++) {
- DOFMatrix* mat = probStat[i]->getSystemMatrix(j, k);
- if (mat && mat->getBoundaryManager())
- removeBoundaryCondition(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
- }
-
- if (probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager())
- removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
-
- if (probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager())
- removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
- }
- }
-
- // Remove periodic boundaries on elements in mesh.
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
- while (elInfo) {
- elInfo->getElement()->deleteElementData(PERIODIC);
- elInfo = stack.traverseNext(elInfo);
- }
+ removePeriodicBoundaryConditions();
}
@@ -426,7 +404,38 @@ namespace AMDiS {
}
- void MeshDistributor::removeBoundaryCondition(BoundaryIndexMap& boundaryMap)
+ void MeshDistributor::removePeriodicBoundaryConditions()
+ {
+ // Remove periodic boundaries in boundary manager on matrices and vectors.
+ for (unsigned int i = 0; i < probStat.size(); i++) {
+ int nComponents = probStat[i]->getNumComponents();
+
+ for (int j = 0; j < nComponents; j++) {
+ for (int k = 0; k < nComponents; k++) {
+ DOFMatrix* mat = probStat[i]->getSystemMatrix(j, k);
+ if (mat && mat->getBoundaryManager())
+ removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
+ }
+
+ if (probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager())
+ removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
+
+ if (probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager())
+ removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
+ }
+ }
+
+ // Remove periodic boundaries on elements in mesh.
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
+ while (elInfo) {
+ elInfo->getElement()->deleteElementData(PERIODIC);
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+
+
+ void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
{
BoundaryIndexMap::iterator it = boundaryMap.begin();
while (it != boundaryMap.end()) {
@@ -473,6 +482,10 @@ namespace AMDiS {
if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
allBound[it.getRank()].push_back(*it);
+ for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
+ if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
+ allBound[it.getRank()].push_back(*it);
+
// === Check the boundaries and adapt mesh if necessary. ===
@@ -602,9 +615,8 @@ namespace AMDiS {
if (s1 != -1 && s2 != -1) {
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
- while (elInfo && elInfo->getElement() != el) {
- elInfo = stack.traverseNext(elInfo);
- }
+ while (elInfo && elInfo->getElement() != el)
+ elInfo = stack.traverseNext(elInfo);
meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
return meshChanged;
@@ -1028,10 +1040,10 @@ namespace AMDiS {
// We have found an element that is at an interior, periodic boundary.
AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
b = bound;
- if (mpiRank > otherElementRank)
- b.rankObj.setReverseMode(b.neighObj, feSpace);
- else
- b.neighObj.setReverseMode(b.rankObj, feSpace);
+ if (mpiRank > otherElementRank)
+ b.rankObj.setReverseMode(b.neighObj, feSpace);
+ else
+ b.neighObj.setReverseMode(b.rankObj, feSpace);
} else {
// We have found an element that is at an interior, non-periodic boundary.
if (mpiRank > otherElementRank) {
@@ -1353,6 +1365,11 @@ namespace AMDiS {
AtomicBoundary& b = myIntBoundary.getNewAtomic(it->first);
b.rankObj = rankDofs[it->second[i]];
b.neighObj = stdMpiObjs.getRecvData(it->first)[i];
+
+ if (mpiRank == 3 && it->first == 0) {
+ MSG("ADDED VERTEX BOUNDARY WITH RANK 0!: %d %d %d\n",
+ it->second[i], b.rankObj.el->getIndex(), b.rankObj.ithObj);
+ }
sendObjects[it->first].push_back(b.rankObj);
}
@@ -1641,7 +1658,7 @@ namespace AMDiS {
for (DofContainer::iterator dofIt = it->second.begin();
dofIt != it->second.end(); ++dofIt)
if (oldVertexDof[*dofIt] && vertexDof[*dofIt])
- sendDofs[it->first].push_back(*dofIt);
+ sendDofs[it->first].push_back(*dofIt);
for (RankToDofContainer::iterator it = oldRecvDofs.begin();
@@ -1666,7 +1683,7 @@ namespace AMDiS {
it->rankObj.el->getNonVertexDofs(feSpace, it->rankObj, dofs);
for (int i = 0; i < static_cast<int>(dofs.size()); i++)
- sendDofs[it.getRank()].push_back(dofs[i]);
+ sendDofs[it.getRank()].push_back(dofs[i]);
}
@@ -1994,7 +2011,7 @@ namespace AMDiS {
boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, tmpdofs);
for (unsigned int i = 0; i < tmpdofs.size(); i++)
- dofs.push_back(tmpdofs[i]);
+ dofs.push_back(tmpdofs[i]);
}
// Added the received dofs to the mapping.
@@ -2184,7 +2201,7 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
- ElementFileWriter::writeFile(vec, feSpace, filename);
+ ElementFileWriter::writeFile(vec, mesh, filename);
}
}
diff --git a/AMDiS/src/parallel/ParallelDomainBase.h b/AMDiS/src/parallel/ParallelDomainBase.h
index d21dcd68ac5875831c8f172c2b7c6d70d93c9529..1849a97d7cc23c7efbf40260561cd1af247879a1 100644
--- a/AMDiS/src/parallel/ParallelDomainBase.h
+++ b/AMDiS/src/parallel/ParallelDomainBase.h
@@ -298,8 +298,12 @@ namespace AMDiS {
*/
void writePartitioningMesh(std::string filename);
+ /// Removes all periodic boundary condition information from all matrices and
+ /// vectors of all stationary problems and from the mesh itself.
+ void removePeriodicBoundaryConditions();
+
// Removes all periodic boundaries from a given boundary map.
- void removeBoundaryCondition(BoundaryIndexMap& boundaryMap);
+ void removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap);
bool fitElementToMeshCode(MeshStructure &code,
Element *el,
diff --git a/AMDiS/src/parallel/ParallelDomainDbg.cc b/AMDiS/src/parallel/ParallelDomainDbg.cc
index eeed93bad3dda296bb91daf372eebcea827df0d7..4391379c8ebe95df6878cbfdda7132d868040ee1 100644
--- a/AMDiS/src/parallel/ParallelDomainDbg.cc
+++ b/AMDiS/src/parallel/ParallelDomainDbg.cc
@@ -144,18 +144,19 @@ namespace AMDiS {
MPI::Request::Waitall(requestCounter, request);
+ int foundError = 0;
for (int i = 0; i < pdb.mpiSize; i++) {
if (i == pdb.mpiRank)
continue;
if (recvSize[i] != recvSizeBuffer[i]) {
- std::cout << "Error: MPI rank " << pdb.mpiRank << " expectes to receive "
- << recvSize[i] << " DOFs from rank " << i << ". But this rank sends "
- << recvSizeBuffer[i] << " DOFs!" << std::endl;
-
- ERROR_EXIT("Should not happen!\n");
+ ERROR("MPI rank %d expectes to receive %d DOFs from rank %d. But this rank sends %d DOFs!\n",
+ pdb.mpiRank, recvSize[i], i, recvSizeBuffer[i]);
+ foundError = 1;
}
}
+ mpi::globalAdd(foundError);
+ TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
// === Now we know that the number of send and received DOFs fits together. ===
// === So we can check if also the coordinates of the communicated DOFs are ===
@@ -201,10 +202,15 @@ namespace AMDiS {
debug::printInfoByDof(pdb.feSpace, *(pdb.recvDofs[it->first][i]));
}
- TEST_EXIT(c)("Wrong DOFs in rank %d!\n", pdb.mpiRank);
+ if (!c) {
+ ERROR("Wrong DOFs in rank %d!\n", pdb.mpiRank);
+ foundError = 1;
+ }
}
}
}
+ mpi::globalAdd(foundError);
+ TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
INFO(pdb.info, 8)("Test common dofs needed %.5f seconds\n", TIME_USED(first, clock()));
}