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

Fixed a very special refinement bug occured only in parallel computations.

parent 067ba67d
...@@ -49,10 +49,10 @@ namespace AMDiS { ...@@ -49,10 +49,10 @@ namespace AMDiS {
{ {
public: public:
/// Type of scalars in the underlying matrix /// Type of scalars in the underlying matrix
typedef double value_type; typedef double value_type;
/// Type of underlying matrix /// 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; /// Type of inserter for the base matrix;
typedef mtl::matrix::inserter<base_matrix_type, mtl::operations::update_plus<value_type> > inserter_type; typedef mtl::matrix::inserter<base_matrix_type, mtl::operations::update_plus<value_type> > inserter_type;
......
...@@ -62,11 +62,11 @@ namespace AMDiS { ...@@ -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; std::map<int, double> vec;
TraverseStack stack; 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) { while (elInfo) {
int index = elInfo->getElement()->getIndex(); int index = elInfo->getElement()->getIndex();
...@@ -74,16 +74,15 @@ namespace AMDiS { ...@@ -74,16 +74,15 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
ElementFileWriter::writeFile(vec, feSpace, filename); ElementFileWriter::writeFile(vec, mesh, filename);
} }
void highlightElementIndexMesh(FiniteElemSpace *feSpace, int idx, void highlightElementIndexMesh(Mesh *mesh, int idx, std::string filename)
std::string filename)
{ {
std::map<int, double> vec; std::map<int, double> vec;
TraverseStack stack; 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) { while (elInfo) {
int index = elInfo->getElement()->getIndex(); int index = elInfo->getElement()->getIndex();
...@@ -91,7 +90,7 @@ namespace AMDiS { ...@@ -91,7 +90,7 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
ElementFileWriter::writeFile(vec, feSpace, filename); ElementFileWriter::writeFile(vec, mesh, filename);
} }
......
...@@ -72,10 +72,9 @@ namespace AMDiS { ...@@ -72,10 +72,9 @@ namespace AMDiS {
* \param[in] feSpace The FE space to be used. * \param[in] feSpace The FE space to be used.
* \param[in] filename Name of the file. * \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, void highlightElementIndexMesh(Mesh *mesh, int idx, std::string filename);
std::string filename);
void colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Element *el); void colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Element *el);
......
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
namespace AMDiS { namespace AMDiS {
ElementFileWriter::ElementFileWriter(std::string name_, ElementFileWriter::ElementFileWriter(std::string name_,
const FiniteElemSpace *feSpace_, Mesh *mesh_,
std::map<int, double> &mapvec) std::map<int, double> &mapvec)
: name(name_), : name(name_),
tecplotExt(".plt"), tecplotExt(".plt"),
...@@ -21,8 +21,7 @@ namespace AMDiS { ...@@ -21,8 +21,7 @@ namespace AMDiS {
indexDecimals(3), indexDecimals(3),
tsModulo(1), tsModulo(1),
timestepNumber(-1), timestepNumber(-1),
mesh(feSpace_->getMesh()), mesh(mesh_),
feSpace(feSpace_),
vec(mapvec) vec(mapvec)
{ {
if (name != "") { if (name != "") {
...@@ -93,10 +92,10 @@ namespace AMDiS { ...@@ -93,10 +92,10 @@ namespace AMDiS {
void ElementFileWriter::writeFile(std::map<int, double> &vec, void ElementFileWriter::writeFile(std::map<int, double> &vec,
const FiniteElemSpace *feSpace, Mesh *mesh,
std::string filename) std::string filename)
{ {
ElementFileWriter efw("", feSpace, vec); ElementFileWriter efw("", mesh, vec);
efw.writeVtkValues(filename); efw.writeVtkValues(filename);
} }
......
...@@ -22,7 +22,7 @@ namespace AMDiS { ...@@ -22,7 +22,7 @@ namespace AMDiS {
public: public:
/// Constructor. /// Constructor.
ElementFileWriter(std::string name, ElementFileWriter(std::string name,
const FiniteElemSpace *feSpace, Mesh *mesh,
std::map<int, double> &vec); std::map<int, double> &vec);
/// Implementation of FileWriterInterface::writeFiles(). /// Implementation of FileWriterInterface::writeFiles().
...@@ -33,7 +33,7 @@ namespace AMDiS { ...@@ -33,7 +33,7 @@ namespace AMDiS {
/// Simple writing procedure for one element map. /// Simple writing procedure for one element map.
static void writeFile(std::map<int, double> &vec, static void writeFile(std::map<int, double> &vec,
const FiniteElemSpace *feSpace, Mesh *mesh,
std::string filename); std::string filename);
protected: protected:
...@@ -101,9 +101,6 @@ namespace AMDiS { ...@@ -101,9 +101,6 @@ namespace AMDiS {
/// Mesh used for output. /// Mesh used for output.
Mesh *mesh; Mesh *mesh;
/// fespace used for output.
const FiniteElemSpace *feSpace;
/// Vector that stores the solution. /// Vector that stores the solution.
std::map<int, double> vec; std::map<int, double> vec;
}; };
......
...@@ -620,6 +620,7 @@ namespace AMDiS { ...@@ -620,6 +620,7 @@ namespace AMDiS {
createPrecon(); createPrecon();
} }
void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, std::string name) void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, std::string name)
{ {
FUNCNAME("ProblemScal::writeResidualMesh()"); FUNCNAME("ProblemScal::writeResidualMesh()");
...@@ -634,10 +635,11 @@ namespace AMDiS { ...@@ -634,10 +635,11 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
ElementFileWriter fw(name, this->getFeSpace(), vec); ElementFileWriter fw(name, this->getMesh(), vec);
fw.writeFiles(adaptInfo, true); fw.writeFiles(adaptInfo, true);
} }
void ProblemScal::createPrecon() void ProblemScal::createPrecon()
{ {
std::string preconType("no"); std::string preconType("no");
...@@ -646,15 +648,16 @@ namespace AMDiS { ...@@ -646,15 +648,16 @@ namespace AMDiS {
CreatorInterface<ITL_BasePreconditioner> *preconCreator = CreatorInterface<ITL_BasePreconditioner> *preconCreator =
CreatorMap<ITL_BasePreconditioner>::getCreator(preconType); CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
solver->setLeftPrecon( preconCreator->create(systemMatrix->getBaseMatrix()) ); solver->setLeftPrecon(preconCreator->create(systemMatrix->getBaseMatrix()));
preconType= "no"; preconType= "no";
GET_PARAMETER(0, name + "->solver->right precon", &preconType); GET_PARAMETER(0, name + "->solver->right precon", &preconType);
preconCreator = CreatorMap<ITL_BasePreconditioner>::getCreator(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) void ProblemScal::serialize(std::ostream &out)
{ {
FUNCNAME("ProblemScal::serialize()"); FUNCNAME("ProblemScal::serialize()");
...@@ -663,6 +666,7 @@ namespace AMDiS { ...@@ -663,6 +666,7 @@ namespace AMDiS {
solution->serialize(out); solution->serialize(out);
} }
void ProblemScal::deserialize(std::istream &in) void ProblemScal::deserialize(std::istream &in)
{ {
FUNCNAME("ProblemScal::deserialize()"); FUNCNAME("ProblemScal::deserialize()");
......
...@@ -1538,7 +1538,7 @@ namespace AMDiS { ...@@ -1538,7 +1538,7 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
ElementFileWriter::writeFile(vec, this->getFeSpace(comp), name); ElementFileWriter::writeFile(vec, this->getMesh(comp), name);
} }
......
...@@ -259,6 +259,7 @@ namespace AMDiS { ...@@ -259,6 +259,7 @@ namespace AMDiS {
} }
} }
void RCNeighbourList::removeDOFParent(int index) void RCNeighbourList::removeDOFParent(int index)
{ {
Tetrahedron *el = dynamic_cast<Tetrahedron*>(rclist[index]->el); Tetrahedron *el = dynamic_cast<Tetrahedron*>(rclist[index]->el);
...@@ -297,6 +298,7 @@ namespace AMDiS { ...@@ -297,6 +298,7 @@ namespace AMDiS {
} }
} }
RCNeighbourList *RCNeighbourList::periodicSplit(DegreeOfFreedom *edge[2], RCNeighbourList *RCNeighbourList::periodicSplit(DegreeOfFreedom *edge[2],
DegreeOfFreedom *nextEdge[2], DegreeOfFreedom *nextEdge[2],
int *n_neigh, int *n_neigh,
......
...@@ -12,10 +12,11 @@ ...@@ -12,10 +12,11 @@
#include "Projection.h" #include "Projection.h"
#include "LeafData.h" #include "LeafData.h"
#include "VertexVector.h" #include "VertexVector.h"
#include "Debug.h"
namespace AMDiS { namespace AMDiS {
ElInfo* RefinementManager2d::refineFunction(ElInfo* el_info) ElInfo* RefinementManager2d::refineFunction(ElInfo* elInfo)
{ {
FUNCNAME("RefinementManager::refineFunction()"); FUNCNAME("RefinementManager::refineFunction()");
...@@ -23,41 +24,38 @@ namespace AMDiS { ...@@ -23,41 +24,38 @@ namespace AMDiS {
DegreeOfFreedom *edge[2]; DegreeOfFreedom *edge[2];
RCNeighbourList* refineList = new RCNeighbourList(2); RCNeighbourList* refineList = new RCNeighbourList(2);
if (el_info->getElement()->getMark() <= 0) { if (elInfo->getElement()->getMark() <= 0) {
delete refineList; delete refineList;
// Element may not be refined. // 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; int n_neigh = 1;
if (el_info->getProjection(0) && if (elInfo->getProjection(0) &&
el_info->getProjection(0)->getType() == VOLUME_PROJECTION) elInfo->getProjection(0)->getType() == VOLUME_PROJECTION)
newCoords = true; newCoords = true;
// === Give the refinement edge the right orientation. === // === Give the refinement edge the right orientation. ===
if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1, 0)) { if (elInfo->getElement()->getDOF(0, 0) < elInfo->getElement()->getDOF(1, 0)) {
edge[0] = const_cast<int*>(el_info->getElement()->getDOF(0)); edge[0] = const_cast<int*>(elInfo->getElement()->getDOF(0));
edge[1] = const_cast<int*>(el_info->getElement()->getDOF(1)); edge[1] = const_cast<int*>(elInfo->getElement()->getDOF(1));
} else { } else {
edge[1] = const_cast<int*>(el_info->getElement()->getDOF(0)); edge[1] = const_cast<int*>(elInfo->getElement()->getDOF(0));
edge[0] = const_cast<int*>(el_info->getElement()->getDOF(1)); edge[0] = const_cast<int*>(elInfo->getElement()->getDOF(1));
} }
// === Get the refinement patch. === // === Get the refinement patch. ===
if (getRefinePatch(&el_info, edge, 0, refineList, &n_neigh)) { getRefinePatch(&elInfo, 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;
}
// === Check for periodic boundary === // === Check for periodic boundary ===
DegreeOfFreedom *next_edge[2]; DegreeOfFreedom *next_edge[2] = {NULL, NULL};
DegreeOfFreedom *first_edge[2] = {edge[0], edge[1]}; DegreeOfFreedom *first_edge[2] = {edge[0], edge[1]};
DegreeOfFreedom *last_edge[2] = {NULL, NULL}; DegreeOfFreedom *last_edge[2] = {NULL, NULL};
int n_neigh_periodic; int n_neigh_periodic;
...@@ -121,25 +119,25 @@ namespace AMDiS { ...@@ -121,25 +119,25 @@ namespace AMDiS {
delete refineList; 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); int dow = Global::getGeo(WORLD);
Projection *projector = el_info->getProjection(0); Projection *projector = elInfo->getProjection(0);
if (!projector || projector->getType() != VOLUME_PROJECTION) if (!projector || projector->getType() != VOLUME_PROJECTION)
projector = el_info->getProjection(2); projector = elInfo->getProjection(2);
if (el->getFirstChild() && projector && (!el->isNewCoordSet())) { if (el->getFirstChild() && projector && (!el->isNewCoordSet())) {
WorldVector<double> *new_coord = new WorldVector<double>; WorldVector<double> *new_coord = new WorldVector<double>;
for (int j = 0; j < dow; j++) 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); projector->project(*new_coord);
el->setNewCoord(new_coord); el->setNewCoord(new_coord);
...@@ -235,10 +233,10 @@ namespace AMDiS { ...@@ -235,10 +233,10 @@ namespace AMDiS {
void RefinementManager2d::bisectTriangle(Triangle *el, DegreeOfFreedom* newDOFs[3]) void RefinementManager2d::bisectTriangle(Triangle *el, DegreeOfFreedom* newDOFs[3])
{ {
FUNCNAME("RefinementManager2d::bisectTriangle()"); FUNCNAME("RefinementManager2d::bisectTriangle()");
TEST_EXIT_DBG(mesh)("no mesh!\n"); TEST_EXIT_DBG(mesh)("no mesh!\n");
Triangle *child[2]; Triangle *child[2];
child[0] = dynamic_cast<Triangle*>(mesh->createNewElement(el)); child[0] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
child[1] = dynamic_cast<Triangle*>(mesh->createNewElement(el)); child[1] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
...@@ -302,36 +300,47 @@ namespace AMDiS { ...@@ -302,36 +300,47 @@ namespace AMDiS {
} }
int RefinementManager2d::getRefinePatch(ElInfo **el_info, void RefinementManager2d::getRefinePatch(ElInfo **elInfo,
DegreeOfFreedom *edge[2], DegreeOfFreedom *edge[2],
int dir, RCNeighbourList* refineList, int dir, RCNeighbourList* refineList,
int *n_neigh) int *n_neigh)
{ {
FUNCNAME("RefinementManager2d::getRefinePatch()"); 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 // Neighbour is not compatible devisible; refine neighbour first, store the
// opp_vertex to traverse back to el. // 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->getElement()->setMark(max(neigh_info->getElement()->getMark(), 1));
neigh_info = refineFunction(neigh_info); neigh_info = refineFunction(neigh_info);
// Now go back to the original element and refine the patch. // === Now go back to the original element and refine the patch. ===
*el_info = neigh_info = stack->traverseNeighbour2d(neigh_info, opp_vertex);
// 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() == 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"); ("invalid traverse_neighbour1\n");
} }
if (refineList->setElement(1, (*el_info)->getNeighbour(2))) { if (refineList->setElement(1, (*elInfo)->getNeighbour(2))) {
TEST_EXIT_DBG((*el_info)->getOppVertex(2) == 2) TEST_EXIT_DBG((*elInfo)->getOppVertex(2) == 2)
("no compatible ref. edge after recursive refinement of neighbour\n"); ("no compatible ref. edge after recursive refinement of neighbour\n");
*n_neigh = 2; *n_neigh = 2;
} }
return 0;
} }
} }
...@@ -54,15 +54,16 @@ namespace AMDiS { ...@@ -54,15 +54,16 @@ namespace AMDiS {
* get_refine_patch returns 1 if the domain's boundary is reached while * get_refine_patch returns 1 if the domain's boundary is reached while
* looping around the refinement edge, otherwise 0 * looping around the refinement edge, otherwise 0
*/ */
int getRefinePatch(ElInfo **el_info, DegreeOfFreedom *edge[2], int dir, void getRefinePatch(ElInfo **elInfo, DegreeOfFreedom *edge[2], int dir,
RCNeighbourList* refineList, int *n_neigh); RCNeighbourList* refineList, int *n_neigh);
/// Refines all elements in the patch. /// 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); int n_neigh, bool bound);
/// Implements RefinementManager::refineFunction. /// Implements RefinementManager::refineFunction.
ElInfo* refineFunction(ElInfo* el_info); ElInfo* refineFunction(ElInfo* elInfo);
/// Refines one Triangle. /// Refines one Triangle.
void bisectTriangle(Triangle *el, DegreeOfFreedom* newDofs[3]); void bisectTriangle(Triangle *el, DegreeOfFreedom* newDofs[3]);
......
...@@ -135,9 +135,8 @@ namespace AMDiS { ...@@ -135,9 +135,8 @@ namespace AMDiS {
currentMacro = traverse_mesh->firstMacroElement(); currentMacro = traverse_mesh->firstMacroElement();
while (((*currentMacro)->getIndex() % maxThreads != myThreadId) && while (((*currentMacro)->getIndex() % maxThreads != myThreadId) &&
(currentMacro != traverse_mesh->endOfMacroElements())) { currentMacro != traverse_mesh->endOfMacroElements())
currentMacro++; currentMacro++;
}
if (currentMacro == traverse_mesh->endOfMacroElements()) if (currentMacro == traverse_mesh->endOfMacroElements())
return NULL; return NULL;
...@@ -349,13 +348,8 @@ namespace AMDiS { ...@@ -349,13 +348,8 @@ namespace AMDiS {
int fillIthChild = info_stack[stack_used]; int fillIthChild = info_stack[stack_used];
info_stack[stack_used]++; info_stack[stack_used]++;
if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) { if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE))
if (fillIthChild == 0) fillIthChild = 1 - fillIthChild;
fillIthChild = 1;
else
fillIthChild = 0;
}
elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]); elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);