Commit 4d67c011 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Several bugfixes. Now, parallelization should also work in 2d.

parent 1f8b7dbf
...@@ -44,7 +44,7 @@ namespace AMDiS { ...@@ -44,7 +44,7 @@ namespace AMDiS {
int myRank = MPI::COMM_WORLD.Get_rank(); int myRank = MPI::COMM_WORLD.Get_rank();
if (rank == -1 || myRank == rank) { if (rank == -1 || myRank == rank) {
DOFVector<double> tmp(feSpace, "tmp"); DOFVector<double> tmp(feSpace, "tmp");
VtkWriter::writeFile(tmp, filename + lexical_cast<std::string>(myRank) + ".vtu"); VtkWriter::writeFile(tmp, filename + lexical_cast<std::string>(myRank) + ".vtu", false);
} }
} }
#endif #endif
......
...@@ -86,41 +86,57 @@ namespace AMDiS { ...@@ -86,41 +86,57 @@ namespace AMDiS {
switch (bound.ithObj) { switch (bound.ithObj) {
case 0: case 0:
{ {
Element *child1 = child[1]; if (child[1] && child[1]->getFirstChild()) {
if (child1 && child1->getFirstChild()) { if (bound.reverseMode) {
nextBound.ithObj = 0; nextBound.ithObj = 1;
child1->getFirstChild()->getVertexDofs(feSpace, nextBound, dofs); child[1]->getSecondChild()->getVertexDofs(feSpace, nextBound, dofs);
dofs.push_back(child[1]->getFirstChild()->getDOF(2));
dofs.push_back(child1->getFirstChild()->getDOF(2)); nextBound.ithObj = 0;
child[1]->getFirstChild()->getVertexDofs(feSpace, nextBound, dofs);
nextBound.ithObj = 1; } else {
child1->getSecondChild()->getVertexDofs(feSpace, nextBound, dofs); nextBound.ithObj = 0;
child[1]->getFirstChild()->getVertexDofs(feSpace, nextBound, dofs);
dofs.push_back(child[1]->getFirstChild()->getDOF(2));
nextBound.ithObj = 1;
child[1]->getSecondChild()->getVertexDofs(feSpace, nextBound, dofs);
}
} }
} }
break; break;
case 1: case 1:
{ {
Element *child0 = child[0]; if (child[0] && child[0]->getFirstChild()) {
if (child0 && child0->getFirstChild()) { if (bound.reverseMode) {
nextBound.ithObj = 0; nextBound.ithObj = 1;
child0->getFirstChild()->getVertexDofs(feSpace, nextBound, dofs); child[0]->getSecondChild()->getVertexDofs(feSpace, nextBound, dofs);
dofs.push_back(child[0]->getFirstChild()->getDOF(2));
dofs.push_back(child[0]->getFirstChild()->getDOF(2)); nextBound.ithObj = 0;
child[0]->getFirstChild()->getVertexDofs(feSpace, nextBound, dofs);
nextBound.ithObj = 1; } else {
child0->getSecondChild()->getVertexDofs(feSpace, nextBound, dofs); nextBound.ithObj = 0;
child[0]->getFirstChild()->getVertexDofs(feSpace, nextBound, dofs);
dofs.push_back(child[0]->getFirstChild()->getDOF(2));
nextBound.ithObj = 1;
child[0]->getSecondChild()->getVertexDofs(feSpace, nextBound, dofs);
}
} }
} }
break; break;
case 2: case 2:
if (child[0]) { if (child[0]) {
nextBound.ithObj = 0; if (bound.reverseMode) {
child[0]->getVertexDofs(feSpace, nextBound, dofs); nextBound.ithObj = 1;
child[1]->getVertexDofs(feSpace, nextBound, dofs);
dofs.push_back(child[0]->getDOF(2)); dofs.push_back(child[0]->getDOF(2));
nextBound.ithObj = 0;
nextBound.ithObj = 1; child[0]->getVertexDofs(feSpace, nextBound, dofs);
child[1]->getVertexDofs(feSpace, nextBound, dofs); } else {
nextBound.ithObj = 0;
child[0]->getVertexDofs(feSpace, nextBound, dofs);
dofs.push_back(child[0]->getDOF(2));
nextBound.ithObj = 1;
child[1]->getVertexDofs(feSpace, nextBound, dofs);
}
} }
break; break;
default: default:
...@@ -159,10 +175,17 @@ namespace AMDiS { ...@@ -159,10 +175,17 @@ namespace AMDiS {
break; break;
case 2: case 2:
if (child[0]) { if (child[0]) {
nextBound.ithObj = 0; if (bound.reverseMode) {
child[0]->getNonVertexDofs(feSpace, nextBound, dofs); nextBound.ithObj = 1;
nextBound.ithObj = 1; child[1]->getNonVertexDofs(feSpace, nextBound, dofs);
child[1]->getNonVertexDofs(feSpace, nextBound, dofs); nextBound.ithObj = 0;
child[0]->getNonVertexDofs(feSpace, nextBound, dofs);
} else {
nextBound.ithObj = 0;
child[0]->getNonVertexDofs(feSpace, nextBound, dofs);
nextBound.ithObj = 1;
child[1]->getNonVertexDofs(feSpace, nextBound, dofs);
}
} else { } else {
addThisEdge = true; addThisEdge = true;
} }
...@@ -173,13 +196,22 @@ namespace AMDiS { ...@@ -173,13 +196,22 @@ namespace AMDiS {
} }
if (addThisEdge) { if (addThisEdge) {
DofContainer addDofs;
ElementDofIterator elDofIter(feSpace, true); ElementDofIterator elDofIter(feSpace, true);
elDofIter.reset(this); elDofIter.reset(this);
do { do {
if (elDofIter.getCurrentPos() == 1 && if (elDofIter.getCurrentPos() == 1 &&
elDofIter.getCurrentElementPos() == bound.ithObj) elDofIter.getCurrentElementPos() == bound.ithObj)
dofs.push_back(elDofIter.getDofPtr()); addDofs.push_back(elDofIter.getDofPtr());
} while(elDofIter.next()); } while(elDofIter.next());
if (bound.reverseMode)
for (int i = addDofs.size() - 1; i >= 0; i--)
dofs.push_back(addDofs[i]);
else
for (unsigned int i = 0; i < addDofs.size(); i++)
dofs.push_back(addDofs[i]);
} }
} }
......
...@@ -128,7 +128,9 @@ namespace AMDiS { ...@@ -128,7 +128,9 @@ namespace AMDiS {
return 0; return 0;
} }
void VtkWriter::writeFile(DOFVector<double> *values, std::string filename) void VtkWriter::writeFile(DOFVector<double> *values,
std::string filename,
bool writeParallel)
{ {
FUNCNAME("VtkWriter::writeFile()"); FUNCNAME("VtkWriter::writeFile()");
...@@ -138,16 +140,18 @@ namespace AMDiS { ...@@ -138,16 +140,18 @@ namespace AMDiS {
VtkWriter writer(&dcList); VtkWriter writer(&dcList);
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
using boost::lexical_cast; if (writeParallel) {
using boost::lexical_cast;
int sPos = filename.find(".vtu");
TEST_EXIT(sPos >= 0)("Failed to find file postfix!\n"); int sPos = filename.find(".vtu");
std::string name = filename.substr(0, sPos); TEST_EXIT(sPos >= 0)("Failed to find file postfix!\n");
std::string name = filename.substr(0, sPos);
if (MPI::COMM_WORLD.Get_rank() == 0)
writer.writeParallelFile(name + ".pvtu", MPI::COMM_WORLD.Get_size(), name, ".vtu"); if (MPI::COMM_WORLD.Get_rank() == 0)
writer.writeParallelFile(name + ".pvtu", MPI::COMM_WORLD.Get_size(), name, ".vtu");
filename = name + "-p" + lexical_cast<std::string>(MPI::COMM_WORLD.Get_rank()) + "-.vtu";
filename = name + "-p" + lexical_cast<std::string>(MPI::COMM_WORLD.Get_rank()) + "-.vtu";
}
#endif #endif
writer.writeFile(filename); writer.writeFile(filename);
......
...@@ -55,12 +55,16 @@ namespace AMDiS { ...@@ -55,12 +55,16 @@ namespace AMDiS {
std::string fnPrefix, std::string fnPostfix); std::string fnPrefix, std::string fnPostfix);
/// May be used to simply write ParaView files. /// May be used to simply write ParaView files.
static void writeFile(DOFVector<double> *values, std::string filename); static void writeFile(DOFVector<double> *values,
std::string filename,
bool writeParallel = true);
/// May be used to simply write ParaView files. /// May be used to simply write ParaView files.
static void writeFile(DOFVector<double> &values, std::string filename) static void writeFile(DOFVector<double> &values,
std::string filename,
bool writeParallel = true)
{ {
writeFile(&values, filename); writeFile(&values, filename, writeParallel);
} }
/// Set a compressing method for file output. /// Set a compressing method for file output.
......
...@@ -12,24 +12,22 @@ namespace AMDiS { ...@@ -12,24 +12,22 @@ namespace AMDiS {
bool otherMode = false; bool otherMode = false;
const BasisFunction *basFcts = feSpace->getBasisFcts();
int nBasFcts = basFcts->getNumber();
std::vector<DegreeOfFreedom> localDofs0(nBasFcts), localDofs1(nBasFcts);
switch (feSpace->getMesh()->getDim()) { switch (feSpace->getMesh()->getDim()) {
case 2: case 2:
if (ithObj == 2 || ithObj != otherBound.ithObj) otherMode = true;
otherMode = true;
break; break;
case 3: case 3:
if (ithObj == 2 || ithObj == 3) { if (ithObj == 2 || ithObj == 3) {
basFcts->getLocalIndices(el, feSpace->getAdmin(), localDofs0); const BasisFunction *basFcts = feSpace->getBasisFcts();
int nBasFcts = basFcts->getNumber();
std::vector<DegreeOfFreedom> localDofs0(nBasFcts), localDofs1(nBasFcts);
basFcts->getLocalIndices(el, feSpace->getAdmin(), localDofs0);
basFcts->getLocalIndices(otherBound.el, feSpace->getAdmin(), localDofs1); basFcts->getLocalIndices(otherBound.el, feSpace->getAdmin(), localDofs1);
otherMode = (localDofs0[0] != localDofs1[0]); otherMode = (localDofs0[0] != localDofs1[0]);
} }
break; break;
default: default:
ERROR_EXIT("This should not happen!\n"); ERROR_EXIT("This should not happen!\n");
} }
......
...@@ -170,6 +170,10 @@ namespace AMDiS { ...@@ -170,6 +170,10 @@ namespace AMDiS {
if (globalRefinement > 0) { if (globalRefinement > 0) {
refineManager->globalRefine(mesh, globalRefinement); refineManager->globalRefine(mesh, globalRefinement);
#if (DEBUG != 0)
debug::writeMesh(feSpace, -1, "grmesh");
#endif
updateLocalGlobalNumbering(); updateLocalGlobalNumbering();
// === Update periodic mapping, if there are periodic boundaries. === // === Update periodic mapping, if there are periodic boundaries. ===
...@@ -177,7 +181,6 @@ namespace AMDiS { ...@@ -177,7 +181,6 @@ namespace AMDiS {
createPeriodicMap(); createPeriodicMap();
} }
/// === Set DOF rank information to all matrices and vectors. === /// === Set DOF rank information to all matrices and vectors. ===
for (int i = 0; i < nComponents; i++) { for (int i = 0; i < nComponents; i++) {
...@@ -910,7 +913,7 @@ namespace AMDiS { ...@@ -910,7 +913,7 @@ namespace AMDiS {
(el->getElementData(PARTITION_ED)); (el->getElementData(PARTITION_ED));
// In 2d and 3d, traverse all vertices of current element. // In 2d and 3d, traverse all vertices of current element.
for (int i = 0; i < 4; i++) { for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
DegreeOfFreedom dof0 = localIndices[i]; DegreeOfFreedom dof0 = localIndices[i];
// Update the owner of the current vertex DOF. // Update the owner of the current vertex DOF.
dofOwner[dof0] = max(dofOwner[dof0], partitionVec[el->getIndex()]); dofOwner[dof0] = max(dofOwner[dof0], partitionVec[el->getIndex()]);
...@@ -955,8 +958,8 @@ namespace AMDiS { ...@@ -955,8 +958,8 @@ namespace AMDiS {
Element *el = it->rankObj.el; Element *el = it->rankObj.el;
basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
for (int j = 0; j < mesh->getGeo(VERTEX); j++) { for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
int dofNo = el->getVertexOfPosition(it->rankObj.subObj, it->rankObj.ithObj, j); int dofNo = el->getVertexOfPosition(it->rankObj.subObj, it->rankObj.ithObj, i);
DegreeOfFreedom dof = localIndices[dofNo]; DegreeOfFreedom dof = localIndices[dofNo];
if (dofOwner[dof] > mpiRank) if (dofOwner[dof] > mpiRank)
...@@ -966,8 +969,8 @@ namespace AMDiS { ...@@ -966,8 +969,8 @@ namespace AMDiS {
} }
if (dim == 3) { if (dim == 3) {
for (int j = 0; j < 3; j++) { for (int i = 0; i < 3; i++) {
int edgeNo = el->getEdgeOfFace(it->rankObj.ithObj, j); int edgeNo = el->getEdgeOfFace(it->rankObj.ithObj, i);
DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)]; DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)];
DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)]; DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)];
GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1));
...@@ -989,8 +992,8 @@ namespace AMDiS { ...@@ -989,8 +992,8 @@ namespace AMDiS {
Element *el = it->rankObj.el; Element *el = it->rankObj.el;
basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
for (int j = 0; j < mesh->getGeo(VERTEX); j++) { for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
int dofNo = el->getVertexOfPosition(it->rankObj.subObj, it->rankObj.ithObj, j); int dofNo = el->getVertexOfPosition(it->rankObj.subObj, it->rankObj.ithObj, i);
DegreeOfFreedom dof = localIndices[dofNo]; DegreeOfFreedom dof = localIndices[dofNo];
...@@ -1001,8 +1004,8 @@ namespace AMDiS { ...@@ -1001,8 +1004,8 @@ namespace AMDiS {
} }
if (dim == 3) { if (dim == 3) {
for (int j = 0; j < 3; j++) { for (int i = 0; i < 3; i++) {
int edgeNo = el->getEdgeOfFace(it->rankObj.ithObj, j); int edgeNo = el->getEdgeOfFace(it->rankObj.ithObj, i);
DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)]; DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)];
DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)]; DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)];
GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1));
...@@ -1015,10 +1018,10 @@ namespace AMDiS { ...@@ -1015,10 +1018,10 @@ namespace AMDiS {
} }
} }
// === Create the new interior boundaries consisting only of edges. This === // === Create the new interior boundaries consisting only of edges. These ===
// === boundaries are created on that ranks, which do not own the boundary === // === boundaries are created on that ranks, which do not own the boundary ===
// === but are on the other side of the edge. Than, theses ranks inform the === // === but are on the other side of the edge. Then, these ranks inform the ===
// === owner of the edge, that they both share them. === // === owner of the edge. ===
// Maps that stores to each rank number the global edges this rank has an // Maps that stores to each rank number the global edges this rank has an
// interior boundary with. // interior boundary with.
...@@ -1030,13 +1033,13 @@ namespace AMDiS { ...@@ -1030,13 +1033,13 @@ namespace AMDiS {
recvObjs[i].clear(); recvObjs[i].clear();
} }
// Traverse all edges in rank's domain. // Traverse all edges in ranks domain.
for (std::map<GlobalEdge, BoundaryObject>::iterator it = rankEdges.begin(); for (std::map<GlobalEdge, BoundaryObject>::iterator it = rankEdges.begin();
it != rankEdges.end(); ++it) { it != rankEdges.end(); ++it) {
// If we have found an edge in rank's domain that has an owner with a higher // If we have found an edge in ranks domain that has an owner with a higher
// rank number, i.e., it must be part of an interior domain, but have not found // rank number, i.e., it must be part of an interior domain, but have not found
// it before to be part of an interior domain, we have found an edge interior // it before to be part of an interior domain, we have found an edge interior
// domain we have to add to the corresponding interior domain object. // boundary.
if (edgeOwner[it->first] > mpiRank && rankBoundaryEdges.count(it->first) == 0) { if (edgeOwner[it->first] > mpiRank && rankBoundaryEdges.count(it->first) == 0) {
std::vector<GlobalEdge> &ownerEdges = recvEdges[edgeOwner[it->first]]; std::vector<GlobalEdge> &ownerEdges = recvEdges[edgeOwner[it->first]];
...@@ -1472,14 +1475,15 @@ namespace AMDiS { ...@@ -1472,14 +1475,15 @@ namespace AMDiS {
DofContainer dofs; DofContainer dofs;
it->rankObj.el->getVertexDofs(feSpace, it->rankObj, dofs); it->rankObj.el->getVertexDofs(feSpace, it->rankObj, dofs);
it->rankObj.el->getNonVertexDofs(feSpace, it->rankObj, dofs); it->rankObj.el->getNonVertexDofs(feSpace, it->rankObj, dofs);
for (int i = 0; i < static_cast<int>(dofs.size()); i++) 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]);
} }
for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) { for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) {
DofContainer dofs; DofContainer dofs;
it->rankObj.el->getNonVertexDofs(feSpace, it->rankObj, dofs);
it->rankObj.el->getVertexDofs(feSpace, it->rankObj, dofs); it->rankObj.el->getVertexDofs(feSpace, it->rankObj, dofs);
it->rankObj.el->getNonVertexDofs(feSpace, it->rankObj, dofs);
for (int i = 0; i < static_cast<int>(dofs.size()); i++) { for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
DofContainer::iterator eraseIt = DofContainer::iterator eraseIt =
...@@ -1488,8 +1492,8 @@ namespace AMDiS { ...@@ -1488,8 +1492,8 @@ namespace AMDiS {
rankDofs.erase(eraseIt); rankDofs.erase(eraseIt);
recvDofs[it.getRank()].push_back(dofs[i]); recvDofs[it.getRank()].push_back(dofs[i]);
} }
} }
nRankDofs = rankDofs.size(); nRankDofs = rankDofs.size();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment