Commit 165d9815 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Several bugfixes in adaptivity within parallel code.

parent 29e49b7e
...@@ -28,19 +28,14 @@ namespace AMDiS { ...@@ -28,19 +28,14 @@ namespace AMDiS {
TEST_EXIT_DBG(child[0]->getMark() < 0 && child[1]->getMark() < 0) TEST_EXIT_DBG(child[0]->getMark() < 0 && child[1]->getMark() < 0)
("element %d with children[%d,%d] must not be coarsend!\n", ("element %d with children[%d,%d] must not be coarsend!\n",
el->getIndex(), child[0]->getIndex(), child[1]->getIndex()); el->getIndex(), child[0]->getIndex(), child[1]->getIndex());
if (mesh->getNumberOfDOFs(EDGE)) { // remove dof from common edge of child[0] and child[1]
/****************************************************************************/ if (mesh->getNumberOfDOFs(EDGE))
/* remove dof from common edge of child[0] and child[1] */
/****************************************************************************/
mesh->freeDOF(const_cast<int*>(child[0]->getDOF(4)), EDGE); mesh->freeDOF(const_cast<int*>(child[0]->getDOF(4)), EDGE);
}
// remove dof from the barycenters of child[0] and child[1]
if (mesh->getNumberOfDOFs(CENTER)) { if (mesh->getNumberOfDOFs(CENTER)) {
/****************************************************************************/ int node = mesh->getNode(CENTER);
/* remove dof from the barycenters of child[0] and child[1] */
/****************************************************************************/
int node = mesh->getNode(CENTER);
mesh->freeDOF(const_cast<int*>(child[0]->getDOF(node)), CENTER); mesh->freeDOF(const_cast<int*>(child[0]->getDOF(node)), CENTER);
mesh->freeDOF(const_cast<int*>(child[1]->getDOF(node)), CENTER); mesh->freeDOF(const_cast<int*>(child[1]->getDOF(node)), CENTER);
...@@ -70,8 +65,10 @@ namespace AMDiS { ...@@ -70,8 +65,10 @@ namespace AMDiS {
int n_neigh, int n_neigh,
int bound) int bound)
{ {
Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>(coarsenList->getElement(0))); Triangle *el =
Triangle *neigh = dynamic_cast<Triangle*>(const_cast<Element*>(coarsenList->getElement(1))); dynamic_cast<Triangle*>(const_cast<Element*>(coarsenList->getElement(0)));
Triangle *neigh =
dynamic_cast<Triangle*>(const_cast<Element*>(coarsenList->getElement(1)));
DegreeOfFreedom *dof[3]; DegreeOfFreedom *dof[3];
dof[0] = const_cast<int*>(el->getChild(0)->getDOF(2)); dof[0] = const_cast<int*>(el->getChild(0)->getDOF(2));
...@@ -82,34 +79,28 @@ namespace AMDiS { ...@@ -82,34 +79,28 @@ namespace AMDiS {
dof[1] = dof[2] = 0; dof[1] = dof[2] = 0;
} }
int node = mesh->getNode(EDGE);
if (mesh->getNumberOfDOFs(EDGE)) { if (mesh->getNumberOfDOFs(EDGE)) {
/****************************************************************************/ int node = mesh->getNode(EDGE);
/* get new dof on el at the midpoint of the coarsening edge */ // get new dof on el at the midpoint of the coarsening edge
/****************************************************************************/
if (!el->getDOF(node + 2)) { if (!el->getDOF(node + 2)) {
el->setDOF(node + 2, mesh->getDOF(EDGE)); el->setDOF(node + 2, mesh->getDOF(EDGE));
if (neigh) { if (neigh)
neigh->setDOF(node + 2, const_cast<int*>(el->getDOF(node+2))); neigh->setDOF(node + 2, const_cast<int*>(el->getDOF(node + 2)));
}
} }
} }
if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(CENTER)) { if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(CENTER))
coarsenList->addDOFParents(n_neigh); coarsenList->addDOFParents(n_neigh);
}
// restrict dof vectors to the parents on the patch
/****************************************************************************/
/* restrict dof vectors to the parents on the patch */
/****************************************************************************/
int nrAdmin = mesh->getNumberOfDOFAdmin(); int nrAdmin = mesh->getNumberOfDOFAdmin();
for (int iadmin = 0; iadmin < nrAdmin; iadmin++) { for (int iadmin = 0; iadmin < nrAdmin; iadmin++) {
std::list<DOFIndexedBase*>::iterator it;
DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin)); DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin));
std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed(); for (std::list<DOFIndexedBase*>::iterator it = admin->beginDOFIndexed();
for (it = admin->beginDOFIndexed(); it != end; ++it) { it != admin->endDOFIndexed(); ++it)
(*it)->coarseRestrict(*coarsenList, n_neigh); (*it)->coarseRestrict(*coarsenList, n_neigh);
}
} }
coarsenTriangle(el); coarsenTriangle(el);
...@@ -117,9 +108,8 @@ namespace AMDiS { ...@@ -117,9 +108,8 @@ namespace AMDiS {
if (neigh) if (neigh)
coarsenTriangle(neigh); coarsenTriangle(neigh);
/****************************************************************************/ // now, remove those dofs in the coarcening edge
/* now, remove those dofs in the coarcening edge */
/****************************************************************************/
mesh->freeDOF(dof[0], VERTEX); mesh->freeDOF(dof[0], VERTEX);
if (mesh->getNumberOfDOFs(EDGE)) { if (mesh->getNumberOfDOFs(EDGE)) {
mesh->freeDOF(dof[1], EDGE); mesh->freeDOF(dof[1], EDGE);
...@@ -145,24 +135,18 @@ namespace AMDiS { ...@@ -145,24 +135,18 @@ namespace AMDiS {
return 0; // single leaves don't get coarsened return 0; // single leaves don't get coarsened
if (el->getChild(0)->getMark() >= 0 || el->getChild(1)->getMark() >= 0) { if (el->getChild(0)->getMark() >= 0 || el->getChild(1)->getMark() >= 0) {
/****************************************************************************/ // one of the children must not be coarsend; return :(
/* one of the children must not be coarsend; return :-( */
/****************************************************************************/
el->setMark(0); el->setMark(0);
return 0; return 0;
} }
if (!el->getChild(0)->isLeaf() || !el->getChild(1)->isLeaf()) { if (!el->getChild(0)->isLeaf() || !el->getChild(1)->isLeaf()) {
/****************************************************************************/ // one of the children is not a leaf element; try again later on
/* one of the children is not a leaf element; try again later on */
/****************************************************************************/
doMore = true; doMore = true;
return 0; return 0;
} }
/****************************************************************************/ // give the refinement edge the right orientation
/* give the refinement edge the right orientation */
/****************************************************************************/
if (el->getDOF(0,0) < el->getDOF(1,0)) { if (el->getDOF(0,0) < el->getDOF(1,0)) {
edge[0] = const_cast<int*>(el->getDOF(0)); edge[0] = const_cast<int*>(el->getDOF(0));
...@@ -180,9 +164,7 @@ namespace AMDiS { ...@@ -180,9 +164,7 @@ namespace AMDiS {
coarse_list.setCoarsePatch(1, el_info->getOppVertex(2) == 2); coarse_list.setCoarsePatch(1, el_info->getOppVertex(2) == 2);
} }
/****************************************************************************/ // check wether we can coarsen the patch or not
/* check wether we can coarsen the patch or not */
/****************************************************************************/
// ========================================================================== // ==========================================================================
// === check for periodic boundary ========================================== // === check for periodic boundary ==========================================
......
...@@ -9,6 +9,9 @@ namespace AMDiS { ...@@ -9,6 +9,9 @@ namespace AMDiS {
template<> template<>
void DOFVector<double>::coarseRestrict(RCNeighbourList& list, int n) void DOFVector<double>::coarseRestrict(RCNeighbourList& list, int n)
{ {
MSG("Coarsen operation \"%d\" in vector \"%s\"\n",
coarsenOperation, this->name.c_str());
switch (coarsenOperation) { switch (coarsenOperation) {
case NO_OPERATION: case NO_OPERATION:
return; return;
...@@ -26,7 +29,9 @@ namespace AMDiS { ...@@ -26,7 +29,9 @@ namespace AMDiS {
(const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n); (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
break; break;
default: default:
ERROR_EXIT("invalid coarsen operation\n"); std::cout << "ERROR WITH PTR = " << this << std::endl;
ERROR_EXIT("Invalid coarsen operation \"%d\" in vector \"%s\"\n",
coarsenOperation, this->name.c_str());
} }
} }
......
...@@ -53,7 +53,11 @@ namespace AMDiS { ...@@ -53,7 +53,11 @@ namespace AMDiS {
elementVector(3), elementVector(3),
boundaryManager(NULL), boundaryManager(NULL),
nBasFcts(0) nBasFcts(0)
{} {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
rankDofs = NULL;
#endif
}
DOFVectorBase(const FiniteElemSpace *f, std::string n); DOFVectorBase(const FiniteElemSpace *f, std::string n);
...@@ -213,12 +217,13 @@ namespace AMDiS { ...@@ -213,12 +217,13 @@ namespace AMDiS {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
inline void setRankDofs(std::map<DegreeOfFreedom, bool> &dofmap) inline void setRankDofs(std::map<DegreeOfFreedom, bool> &dofmap)
{ {
// rankDofs = &dofmap; rankDofs = &dofmap;
} }
inline bool isRankDof(DegreeOfFreedom dof) inline bool isRankDof(DegreeOfFreedom dof)
{ {
TEST_EXIT_DBG(rankDofs)("No rank dofs set!\n"); TEST_EXIT_DBG(rankDofs)("No rank dofs set!\n");
return (*rankDofs)[dof]; return (*rankDofs)[dof];
} }
#endif #endif
...@@ -267,7 +272,6 @@ namespace AMDiS { ...@@ -267,7 +272,6 @@ namespace AMDiS {
int dim; int dim;
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
public:
std::map<DegreeOfFreedom, bool> *rankDofs; std::map<DegreeOfFreedom, bool> *rankDofs;
#endif #endif
}; };
......
...@@ -630,6 +630,7 @@ namespace AMDiS { ...@@ -630,6 +630,7 @@ namespace AMDiS {
coarsenOperation = rhs.coarsenOperation; coarsenOperation = rhs.coarsenOperation;
this->operators = rhs.operators; this->operators = rhs.operators;
this->operatorFactor = rhs.operatorFactor; this->operatorFactor = rhs.operatorFactor;
if (rhs.boundaryManager) { if (rhs.boundaryManager) {
if (this->boundaryManager) if (this->boundaryManager)
delete this->boundaryManager; delete this->boundaryManager;
...@@ -639,6 +640,10 @@ namespace AMDiS { ...@@ -639,6 +640,10 @@ namespace AMDiS {
this->boundaryManager = NULL; this->boundaryManager = NULL;
} }
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
this->rankDofs = rhs.rankDofs;
#endif
return *this; return *this;
} }
......
...@@ -54,7 +54,7 @@ namespace AMDiS { ...@@ -54,7 +54,7 @@ namespace AMDiS {
for (int i = 0; i < nBasFcts; i++) { for (int i = 0; i < nBasFcts; i++) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
// if (vector->isRankDof(dofIndices[i])) if (vector->isRankDof(dofIndices[i]))
#endif #endif
if (localBound[i] == boundaryType) { if (localBound[i] == boundaryType) {
if (f) { if (f) {
......
...@@ -580,17 +580,16 @@ namespace AMDiS { ...@@ -580,17 +580,16 @@ namespace AMDiS {
int s1 = el->getSideOfChild(0, ithSide, elType); int s1 = el->getSideOfChild(0, ithSide, elType);
int s2 = el->getSideOfChild(1, ithSide, elType); int s2 = el->getSideOfChild(1, ithSide, elType);
code.nextElement(); if (s1 != -1) {
code.nextElement();
if (s1 != -1)
fitElementToMeshCode(refineManager, code, fitElementToMeshCode(refineManager, code,
el->getFirstChild(), s1, el->getChildType(elType)); el->getFirstChild(), s1, el->getChildType(elType));
}
code.nextElement(); if (s2 != -1) {
code.nextElement();
if (s2 != -1)
fitElementToMeshCode(refineManager, code, fitElementToMeshCode(refineManager, code,
el->getSecondChild(), s2, el->getChildType(elType)); el->getSecondChild(), s2, el->getChildType(elType));
}
} }
} }
...@@ -53,18 +53,19 @@ namespace AMDiS { ...@@ -53,18 +53,19 @@ namespace AMDiS {
commit(); commit();
} }
void MeshStructure::init(Element *el, int ithSide, int elType) void MeshStructure::init(Element *el, int ithSide, int elType, bool reverseOrder)
{ {
FUNCNAME("MeshStructure::init()"); FUNCNAME("MeshStructure::init()");
clear(); clear();
addAlongSide(el, ithSide, elType); addAlongSide(el, ithSide, elType, reverseOrder);
commit(); commit();
} }
void MeshStructure::addAlongSide(Element *el, int ithSide, int elType) void MeshStructure::addAlongSide(Element *el, int ithSide, int elType,
bool reverseOrder)
{ {
insertElement(el->isLeaf()); insertElement(el->isLeaf());
...@@ -72,10 +73,17 @@ namespace AMDiS { ...@@ -72,10 +73,17 @@ namespace AMDiS {
int s1 = el->getSideOfChild(0, ithSide, elType); int s1 = el->getSideOfChild(0, ithSide, elType);
int s2 = el->getSideOfChild(1, ithSide, elType); int s2 = el->getSideOfChild(1, ithSide, elType);
if (s1 != -1) if (!reverseOrder) {
addAlongSide(el->getFirstChild(), s1, el->getChildType(elType)); if (s1 != -1)
if (s2 != -1) addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), reverseOrder);
addAlongSide(el->getSecondChild(), s2, el->getChildType(elType)); if (s2 != -1)
addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), reverseOrder);
} else {
if (s2 != -1)
addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), reverseOrder);
if (s1 != -1)
addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), reverseOrder);
}
} }
} }
......
...@@ -44,7 +44,7 @@ namespace AMDiS { ...@@ -44,7 +44,7 @@ namespace AMDiS {
/// Creates a mesh structure code from a mesh object by traversing it in preorder. /// Creates a mesh structure code from a mesh object by traversing it in preorder.
void init(Mesh *mesh); void init(Mesh *mesh);
void init(Element *el, int ithSide, int elType); void init(Element *el, int ithSide, int elType, bool reverseOrder);
void init(const std::vector<unsigned long int>& initCode, int n) void init(const std::vector<unsigned long int>& initCode, int n)
{ {
...@@ -137,7 +137,7 @@ namespace AMDiS { ...@@ -137,7 +137,7 @@ namespace AMDiS {
/// Insert a new element to the structure code. Is used by the init function. /// Insert a new element to the structure code. Is used by the init function.
void insertElement(bool isLeaf); void insertElement(bool isLeaf);
void addAlongSide(Element *el, int ithSide, int elType); void addAlongSide(Element *el, int ithSide, int elType, bool reverseOrder);
/// Merges two mesh structure codes to one structure code. /// Merges two mesh structure codes to one structure code.
void merge(MeshStructure *structure1, void merge(MeshStructure *structure1,
......
...@@ -741,6 +741,7 @@ namespace AMDiS { ...@@ -741,6 +741,7 @@ namespace AMDiS {
if (mesh->getChangeIndex() == lastMeshChangeIndex) if (mesh->getChangeIndex() == lastMeshChangeIndex)
return; return;
MSG("MESH CHANGED!\n");
// === Create mesh structure codes for all ranks boundary elements. === // === Create mesh structure codes for all ranks boundary elements. ===
...@@ -755,7 +756,7 @@ namespace AMDiS { ...@@ -755,7 +756,7 @@ namespace AMDiS {
MeshStructure elCode; MeshStructure elCode;
elCode.init(boundIt->rankObj.el, boundIt->rankObj.ithObj, elCode.init(boundIt->rankObj.el, boundIt->rankObj.ithObj,
boundIt->rankObj.elType); boundIt->rankObj.elType, false);
sendCodes[it->first].push_back(elCode); sendCodes[it->first].push_back(elCode);
} }
} }
...@@ -763,11 +764,7 @@ namespace AMDiS { ...@@ -763,11 +764,7 @@ namespace AMDiS {
StdMpi<MeshCodeVec, MeshCodeVec> stdMpi(mpiComm); StdMpi<MeshCodeVec, MeshCodeVec> stdMpi(mpiComm);
stdMpi.prepareCommunication(true); stdMpi.prepareCommunication(true);
stdMpi.send(sendCodes); stdMpi.send(sendCodes);
stdMpi.recv(otherIntBoundary.boundary);
for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
it != otherIntBoundary.boundary.end(); ++it)
stdMpi.recv(it->first);
stdMpi.startCommunication<unsigned long int>(); stdMpi.startCommunication<unsigned long int>();
...@@ -775,6 +772,19 @@ namespace AMDiS { ...@@ -775,6 +772,19 @@ namespace AMDiS {
bool meshFitTogether = true; bool meshFitTogether = true;
if (mpiRank == 0) {
DOFVector<double> ddd(feSpace, "tmp");
ddd.set(0.0);
VtkWriter::writeFile(&ddd, "testold0.vtu");
}
if (mpiRank == 1) {
DOFVector<double> ddd(feSpace, "tmp");
ddd.set(0.0);
VtkWriter::writeFile(&ddd, "testold1.vtu");
}
for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin(); for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
it != otherIntBoundary.boundary.end(); ++it) { it != otherIntBoundary.boundary.end(); ++it) {
...@@ -784,15 +794,21 @@ namespace AMDiS { ...@@ -784,15 +794,21 @@ namespace AMDiS {
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin(); for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end(); ++boundIt) { boundIt != it->second.end(); ++boundIt) {
std::cout << "[" << mpiRank << "] GET CODE: ";
recvCodes[i].print();
MeshStructure elCode; MeshStructure elCode;
elCode.init(boundIt->rankObj.el, boundIt->rankObj.ithObj, elCode.init(boundIt->rankObj.el, boundIt->rankObj.ithObj,
boundIt->rankObj.elType); boundIt->rankObj.elType, true);
std::cout << "[" << mpiRank << "] CALC CODE: ";
elCode.print();
if (elCode.getCode() != recvCodes[i].getCode()) { if (elCode.getCode() != recvCodes[i].getCode()) {
TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n"); TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
// recvCodes[i].reset(); recvCodes[i].reset();
// fitElementToMeshCode(refineManager, recvCodes[i], boundIt->rankObj.el, fitElementToMeshCode(refineManager, recvCodes[i], boundIt->rankObj.el,
// boundIt->rankObj.ithObj, boundIt->rankObj.elType); boundIt->rankObj.ithObj, boundIt->rankObj.elType);
meshFitTogether = false; meshFitTogether = false;
} }
...@@ -801,9 +817,27 @@ namespace AMDiS { ...@@ -801,9 +817,27 @@ namespace AMDiS {
} }
} }
MSG("MESH FIT ALGO FINISHED!\n");
if (mpiRank == 0) {
DOFVector<double> ddd(feSpace, "tmp");
ddd.set(0.0);
VtkWriter::writeFile(&ddd, "test0.vtu");
}
if (mpiRank == 1) {
DOFVector<double> ddd(feSpace, "tmp");
ddd.set(0.0);
VtkWriter::writeFile(&ddd, "test1.vtu");
}
if (!meshFitTogether) { if (!meshFitTogether) {
MSG("MESH STRUCTURE CHANGED!\n");
std::cout << "MESH HAS BEEN CHANGED!" << std::endl; std::cout << "MESH HAS BEEN CHANGED!" << std::endl;
exit(0);
exit(0);
} }
updateLocalGlobalNumbering(); updateLocalGlobalNumbering();
......
...@@ -83,6 +83,10 @@ namespace AMDiS { ...@@ -83,6 +83,10 @@ namespace AMDiS {
if (probVec->getRHS()->getDOFVector(i)->getBoundaryManager()) if (probVec->getRHS()->getDOFVector(i)->getBoundaryManager())
removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probVec->getRHS()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap())); removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probVec->getRHS()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));
std::cout << "TEST SOL " << probVec->getSolution()->getDOFVector(i) << ": " << probVec->getSolution()->getDOFVector(i)->getCoarsenOperation() << std::endl;
std::cout << "TEST RHS " << probVec->getRHS()->getDOFVector(i) << ": " << probVec->getRHS()->getDOFVector(i)->getCoarsenOperation() << std::endl;
} }
// Remove periodic boundaries on elements in mesh. // Remove periodic boundaries on elements in mesh.
...@@ -91,7 +95,7 @@ namespace AMDiS { ...@@ -91,7 +95,7 @@ namespace AMDiS {
while (elInfo) { while (elInfo) {
elInfo->getElement()->deleteElementData(PERIODIC); elInfo->getElement()->deleteElementData(PERIODIC);
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
} }
void ParallelDomainVec::exitParallelization(AdaptInfo *adaptInfo) void ParallelDomainVec::exitParallelization(AdaptInfo *adaptInfo)
......
...@@ -338,7 +338,8 @@ namespace AMDiS { ...@@ -338,7 +338,8 @@ namespace AMDiS {
std::string numberedName = "rhs[" + boost::lexical_cast<std::string>(i) + "]"; std::string numberedName = "rhs[" + boost::lexical_cast<std::string>(i) + "]";
rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], numberedName)); rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], numberedName));
numberedName = name + boost::lexical_cast<std::