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

Small bugfix for 2d domain decomposition.

parent b7b61ade
...@@ -388,13 +388,10 @@ namespace AMDiS { ...@@ -388,13 +388,10 @@ namespace AMDiS {
* \param[in] geoPos Must be either EDGE or FACE. Defines whether an * \param[in] geoPos Must be either EDGE or FACE. Defines whether an
* edge or a face (only in 3d) should be traversed. * edge or a face (only in 3d) should be traversed.
* \param[out] dofs List of dofs, where the result is stored. * \param[out] dofs List of dofs, where the result is stored.
* \param[in] parentVertices If true, also the two vertices of the parent
* element are put into the result list.
*/ */
virtual void getVertexDofs(FiniteElemSpace* feSpace, virtual void getVertexDofs(FiniteElemSpace* feSpace,
BoundaryObject bound, BoundaryObject bound,
DofContainer& dofs, DofContainer& dofs) const = 0;
bool parentVertices = false) const = 0;
/** \brief /** \brief
* Traverses an edge/face of a given element (this includes also all children of * Traverses an edge/face of a given element (this includes also all children of
......
...@@ -154,8 +154,7 @@ namespace AMDiS { ...@@ -154,8 +154,7 @@ namespace AMDiS {
return "Line"; return "Line";
} }
void getVertexDofs(FiniteElemSpace*, BoundaryObject, DofContainer&, void getVertexDofs(FiniteElemSpace*, BoundaryObject, DofContainer&) const
bool parentVertices = false) const
{ {
FUNCNAME("Line::getVertexDofs()"); FUNCNAME("Line::getVertexDofs()");
ERROR_EXIT("Not yet implemented!\n"); ERROR_EXIT("Not yet implemented!\n");
......
...@@ -187,8 +187,7 @@ namespace AMDiS { ...@@ -187,8 +187,7 @@ namespace AMDiS {
void Tetrahedron::getVertexDofs(FiniteElemSpace* feSpace, void Tetrahedron::getVertexDofs(FiniteElemSpace* feSpace,
BoundaryObject bound, BoundaryObject bound,
DofContainer& dofs, DofContainer& dofs) const
bool parentVertices) const
{ {
FUNCNAME("Tetrahedron::getVertexDofs()"); FUNCNAME("Tetrahedron::getVertexDofs()");
...@@ -197,10 +196,10 @@ namespace AMDiS { ...@@ -197,10 +196,10 @@ namespace AMDiS {
dofs.push_back(dof[bound.ithObj]); dofs.push_back(dof[bound.ithObj]);
break; break;
case EDGE: case EDGE:
getVertexDofsAtEdge(feSpace, bound, dofs, parentVertices); getVertexDofsAtEdge(feSpace, bound, dofs);
break; break;
case FACE: case FACE:
getVertexDofsAtFace(feSpace, bound, dofs, parentVertices); getVertexDofsAtFace(feSpace, bound, dofs);
break; break;
default: default:
ERROR_EXIT("Should not happen!\n"); ERROR_EXIT("Should not happen!\n");
...@@ -210,31 +209,8 @@ namespace AMDiS { ...@@ -210,31 +209,8 @@ namespace AMDiS {
void Tetrahedron::getVertexDofsAtFace(FiniteElemSpace* feSpace, void Tetrahedron::getVertexDofsAtFace(FiniteElemSpace* feSpace,
BoundaryObject bound, BoundaryObject bound,
DofContainer& dofs, DofContainer& dofs) const
bool parentVertices) const
{ {
if (parentVertices) {
ERROR_EXIT("Einbau excludedSubstructures!\n");
switch (bound.ithObj) {
case 0:
dofs.push_back(dof[1]); dofs.push_back(dof[2]); dofs.push_back(dof[3]);
break;
case 1:
dofs.push_back(dof[0]); dofs.push_back(dof[2]); dofs.push_back(dof[3]);
break;
case 2:
dofs.push_back(dof[0]); dofs.push_back(dof[1]); dofs.push_back(dof[3]);
break;
case 3:
dofs.push_back(dof[0]); dofs.push_back(dof[1]); dofs.push_back(dof[2]);
break;
default:
ERROR_EXIT("Should never happen!\n");
}
}
if (!child[0]) if (!child[0])
return; return;
...@@ -294,15 +270,10 @@ namespace AMDiS { ...@@ -294,15 +270,10 @@ namespace AMDiS {
void Tetrahedron::getVertexDofsAtEdge(FiniteElemSpace* feSpace, void Tetrahedron::getVertexDofsAtEdge(FiniteElemSpace* feSpace,
BoundaryObject bound, BoundaryObject bound,
DofContainer& dofs, DofContainer& dofs) const
bool parentVertices) const
{ {
FUNCNAME("Tetrahedron::getVertexDofsAtEdge()"); FUNCNAME("Tetrahedron::getVertexDofsAtEdge()");
if (parentVertices) {
ERROR_EXIT("Einbau excludedSubstructures!\n");
}
if (!child[0]) if (!child[0])
return; return;
......
...@@ -137,18 +137,15 @@ namespace AMDiS { ...@@ -137,18 +137,15 @@ namespace AMDiS {
void getVertexDofs(FiniteElemSpace* feSpace, void getVertexDofs(FiniteElemSpace* feSpace,
BoundaryObject bound, BoundaryObject bound,
DofContainer& dofs, DofContainer& dofs) const;
bool parentVertices = false) const;
void getVertexDofsAtFace(FiniteElemSpace* feSpace, void getVertexDofsAtFace(FiniteElemSpace* feSpace,
BoundaryObject bound, BoundaryObject bound,
DofContainer& dofs, DofContainer& dofs) const;
bool parentVertices) const;
void getVertexDofsAtEdge(FiniteElemSpace* feSpace, void getVertexDofsAtEdge(FiniteElemSpace* feSpace,
BoundaryObject bound, BoundaryObject bound,
DofContainer& dofs, DofContainer& dofs) const;
bool parentVertices) const;
void getNonVertexDofs(FiniteElemSpace* feSpace, void getNonVertexDofs(FiniteElemSpace* feSpace,
BoundaryObject bound, BoundaryObject bound,
......
...@@ -77,27 +77,10 @@ namespace AMDiS { ...@@ -77,27 +77,10 @@ namespace AMDiS {
void Triangle::getVertexDofs(FiniteElemSpace* feSpace, void Triangle::getVertexDofs(FiniteElemSpace* feSpace,
BoundaryObject bound, BoundaryObject bound,
DofContainer& dofs, DofContainer& dofs) const
bool parentVertices) const
{ {
FUNCNAME("Triangle::getVertexDofs()"); FUNCNAME("Triangle::getVertexDofs()");
if (parentVertices) {
switch (bound.ithObj) {
case 0:
dofs.push_back(dof[1]);
break;
case 1:
dofs.push_back(dof[0]);
break;
case 2:
dofs.push_back(dof[0]);
break;
default:
ERROR_EXIT("Should never happen!\n");
}
}
BoundaryObject nextBound = bound; BoundaryObject nextBound = bound;
switch (bound.ithObj) { switch (bound.ithObj) {
...@@ -143,22 +126,6 @@ namespace AMDiS { ...@@ -143,22 +126,6 @@ namespace AMDiS {
default: default:
ERROR_EXIT("Should never happen!\n"); ERROR_EXIT("Should never happen!\n");
} }
if (parentVertices) {
switch (bound.ithObj) {
case 0:
dofs.push_back(dof[2]);
break;
case 1:
dofs.push_back(dof[2]);
break;
case 2:
dofs.push_back(dof[1]);
break;
default:
ERROR_EXIT("Should never happen!\n");
}
}
} }
......
...@@ -166,8 +166,7 @@ namespace AMDiS { ...@@ -166,8 +166,7 @@ namespace AMDiS {
void getVertexDofs(FiniteElemSpace* feSpace, void getVertexDofs(FiniteElemSpace* feSpace,
BoundaryObject bound, BoundaryObject bound,
DofContainer& dofs, DofContainer& dofs) const;
bool parentVertices = false) const;
void getNonVertexDofs(FiniteElemSpace* feSpace, void getNonVertexDofs(FiniteElemSpace* feSpace,
BoundaryObject bound, BoundaryObject bound,
......
...@@ -10,7 +10,7 @@ namespace AMDiS { ...@@ -10,7 +10,7 @@ namespace AMDiS {
PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *) PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
{ {
if (iter % 10 == 0 && MPI::COMM_WORLD.Get_rank() == 0) if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
std::cout << "[0] Petsc-Iteration " << iter << ": " << rnorm << std::endl; std::cout << "[0] Petsc-Iteration " << iter << ": " << rnorm << std::endl;
return 0; return 0;
......
...@@ -154,6 +154,16 @@ namespace AMDiS { ...@@ -154,6 +154,16 @@ namespace AMDiS {
} }
} }
inline AtomicBoundary& operator*()
{
return *vecIt;
}
inline AtomicBoundary* operator->()
{
return &(*vecIt);
}
void nextRank() void nextRank()
{ {
++mapIt; ++mapIt;
...@@ -168,11 +178,6 @@ namespace AMDiS { ...@@ -168,11 +178,6 @@ namespace AMDiS {
return mapIt->first; return mapIt->first;
} }
AtomicBoundary& getBound()
{
return *vecIt;
}
protected: protected:
inline void nextNonempty() inline void nextNonempty()
......
...@@ -376,12 +376,12 @@ namespace AMDiS { ...@@ -376,12 +376,12 @@ namespace AMDiS {
RankToBoundMap allBound; RankToBoundMap allBound;
for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim)) if (it->rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
allBound[it.getRank()].push_back(it.getBound()); allBound[it.getRank()].push_back(*it);
for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim)) if (it->rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
allBound[it.getRank()].push_back(it.getBound()); allBound[it.getRank()].push_back(*it);
// === Check the boundaries and adapt mesh if necessary. === // === Check the boundaries and adapt mesh if necessary. ===
...@@ -795,17 +795,7 @@ namespace AMDiS { ...@@ -795,17 +795,7 @@ namespace AMDiS {
int nNeighbours = mesh->getGeo(NEIGH); int nNeighbours = mesh->getGeo(NEIGH);
int dim = mesh->getDim(); int dim = mesh->getDim();
GeoIndex subObj = CENTER; GeoIndex subObj = INDEX_OF_DIM(dim - 1, dim);
switch (dim) {
case 2:
subObj = EDGE;
break;
case 3:
subObj = FACE;
break;
default:
ERROR_EXIT("What is this?\n");
}
// === First, traverse the mesh and search for all elements that have an === // === First, traverse the mesh and search for all elements that have an ===
// === boundary with an element within another partition. === // === boundary with an element within another partition. ===
...@@ -849,16 +839,9 @@ namespace AMDiS { ...@@ -849,16 +839,9 @@ namespace AMDiS {
bound.neighObj.subObj = subObj; bound.neighObj.subObj = subObj;
bound.neighObj.ithObj = elInfo->getSideOfNeighbour(i); bound.neighObj.ithObj = elInfo->getSideOfNeighbour(i);
if (dim == 2) {
// i == 2 => getSideOfNeighbour(i) == 2
TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
("Should not happen!\n");
}
// Get rank number of the neighbouring element. // Get rank number of the neighbouring element.
int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()]; int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()];
// === Add the boundary information object to the corresponding overall === // === Add the boundary information object to the corresponding overall ===
// === boundary. There are three possibilities: if the boundary is a === // === boundary. There are three possibilities: if the boundary is a ===
// === periodic boundary, just add it to \ref periodicBounadry. Here it === // === periodic boundary, just add it to \ref periodicBounadry. Here it ===
...@@ -896,24 +879,22 @@ namespace AMDiS { ...@@ -896,24 +879,22 @@ namespace AMDiS {
{ {
FUNCNAME("ParallelDomainBase::createSubstructureInteriorBoundaryInfo()"); FUNCNAME("ParallelDomainBase::createSubstructureInteriorBoundaryInfo()");
// === Seach for all vertices/edges, which are part of an interior boundary, === // === Seach for all vertices/edges, which are part of an interior boundary, ===
// === but are not a part of the interior boundaries that are created based on === // === but are not part of the interior boundaries that are created based on ===
// === the information of macro elements. === // === the information of macro elements. ===
int dim = mesh->getDim(); int dim = mesh->getDim();
const BasisFunction *basFcts = feSpace->getBasisFcts(); const BasisFunction *basFcts = feSpace->getBasisFcts();
std::vector<DegreeOfFreedom> localIndices(basFcts->getNumber(), 0); std::vector<DegreeOfFreedom> localIndices(basFcts->getNumber());
// Maps each DOF in the whole domain to the rank that ownes this DOF. // Maps each DOF in the whole domain to the rank that ownes this DOF.
std::map<DegreeOfFreedom, int> dofOwner; std::map<DegreeOfFreedom, int> dofOwner;
// Maps each DOF in ranks partition of the domain to the element object that // Maps each DOF in ranks partition to the element object that contains this DOF.
// contains this DOF.
std::map<DegreeOfFreedom, BoundaryObject> rankDofs; std::map<DegreeOfFreedom, BoundaryObject> rankDofs;
// Maps each edge in the whole domain to the rank that ownes this edge. // Maps each edge in the whole domain to the rank that ownes this edge.
std::map<GlobalEdge, int> edgeOwner; std::map<GlobalEdge, int> edgeOwner;
// Maps each edge in ranks partition of the domain to the element object that // Maps each edge in ranks partition to the element object that contains this edge.
// contains this edge.
std::map<GlobalEdge, BoundaryObject> rankEdges; std::map<GlobalEdge, BoundaryObject> rankEdges;
...@@ -927,7 +908,19 @@ namespace AMDiS { ...@@ -927,7 +908,19 @@ namespace AMDiS {
PartitionElementData *partitionData = dynamic_cast<PartitionElementData*> PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
(el->getElementData(PARTITION_ED)); (el->getElementData(PARTITION_ED));
// In 2d and 3d, traverse all vertices of current element.
for (int i = 0; i < 4; i++) {
DegreeOfFreedom dof0 = localIndices[i];
// Update the owner of the current vertex DOF.
dofOwner[dof0] = max(dofOwner[dof0], partitionVec[el->getIndex()]);
// If the DOF is part of an element that is within ranks partition, add it
// to the set of all ranks vertex DOFs.
if (partitionData && partitionData->getPartitionStatus() == IN)
rankDofs[dof0] = BoundaryObject(el, elInfo->getType(), VERTEX, i);
}
// In 3d, traverse all edges of current element. // In 3d, traverse all edges of current element.
if (dim == 3) { if (dim == 3) {
for (int i = 0; i < 6; i++) { for (int i = 0; i < 6; i++) {
...@@ -938,28 +931,18 @@ namespace AMDiS { ...@@ -938,28 +931,18 @@ namespace AMDiS {
// Update the owner of the current edge. // Update the owner of the current edge.
edgeOwner[edge] = max(edgeOwner[edge], partitionVec[el->getIndex()]); edgeOwner[edge] = max(edgeOwner[edge], partitionVec[el->getIndex()]);
// If the edge is part of an element that is part of rank's domain, add it // If the edge is part of an element that is within ranks partition, add it
// to the set of all rank's edges. // to the set of all ranks edges.
if (partitionData && partitionData->getPartitionStatus() == IN) if (partitionData && partitionData->getPartitionStatus() == IN)
rankEdges[edge] = BoundaryObject(el, elInfo->getType(), EDGE, i); rankEdges[edge] = BoundaryObject(el, elInfo->getType(), EDGE, i);
} }
} }
// In 2d and 3d, traverse all vertices of current element.
for (int i = 0; i < 4; i++) {
DegreeOfFreedom dof0 = localIndices[i];
dofOwner[dof0] = max(dofOwner[dof0], partitionVec[el->getIndex()]);
if (partitionData && partitionData->getPartitionStatus() == IN)
rankDofs[dof0] = BoundaryObject(el, elInfo->getType(), VERTEX, i);
}
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
// === Create a set of all edges and vertices at rank's interior boundaries. === // === Create a set of all edges and vertices at ranks interior boundaries. ===
// Stores all edges at rank's interior boundaries. // Stores all edges at rank's interior boundaries.
std::set<GlobalEdge> rankBoundaryEdges; std::set<GlobalEdge> rankBoundaryEdges;
...@@ -969,12 +952,22 @@ namespace AMDiS { ...@@ -969,12 +952,22 @@ namespace AMDiS {
// First, traverse the rank owned elements at the interior boundaries. // First, traverse the rank owned elements at the interior boundaries.
for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) { for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) {
Element *el = it.getBound().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++) {
int dofNo = el->getVertexOfPosition(it->rankObj.subObj, it->rankObj.ithObj, j);
DegreeOfFreedom dof = localIndices[dofNo];
if (dofOwner[dof] > mpiRank)
it->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
else
rankBoundaryDofs.insert(dof);
}
if (dim == 3) { if (dim == 3) {
for (int j = 0; j < 3; j++) { for (int j = 0; j < 3; j++) {
int edgeNo = el->getEdgeOfFace(it.getBound().rankObj.ithObj, j); int edgeNo = el->getEdgeOfFace(it->rankObj.ithObj, j);
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));
...@@ -984,51 +977,42 @@ namespace AMDiS { ...@@ -984,51 +977,42 @@ namespace AMDiS {
// Otherwise, it is part of the interior boundary and we add it to the set // Otherwise, it is part of the interior boundary and we add it to the set
// rankBoundaryEdges. // rankBoundaryEdges.
if (edgeOwner[edge] > mpiRank) if (edgeOwner[edge] > mpiRank)
it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo)); it->rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
else else
rankBoundaryEdges.insert(edge); rankBoundaryEdges.insert(edge);
} }
} }
for (int j = 0; j < 3; j++) {
int dofNo = el->getVertexOfPosition(FACE, it.getBound().rankObj.ithObj, j);
DegreeOfFreedom dof = localIndices[dofNo];
if (dofOwner[dof] > mpiRank)
it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
else
rankBoundaryDofs.insert(dof);
}
} }
// Now, do the same with all other elements at the interior boundaries. // Now, do the same with all other elements at the interior boundaries.
for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) { for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) {
Element *el = it.getBound().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++) {
int dofNo = el->getVertexOfPosition(it->rankObj.subObj, it->rankObj.ithObj, j);
DegreeOfFreedom dof = localIndices[dofNo];
if (dofOwner[dof] > it.getRank())
it->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
else
rankBoundaryDofs.insert(dof);
}
if (dim == 3) { if (dim == 3) {
for (int j = 0; j < 3; j++) { for (int j = 0; j < 3; j++) {
int edgeNo = el->getEdgeOfFace(it.getBound().rankObj.ithObj, j); int edgeNo = el->getEdgeOfFace(it->rankObj.ithObj, j);
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));
if (edgeOwner[edge] > it.getRank()) if (edgeOwner[edge] > it.getRank())
it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo)); it->rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
else else
rankBoundaryEdges.insert(edge); rankBoundaryEdges.insert(edge);
} }
} }
for (int j = 0; j < 3; j++) {
int dofNo = el->getVertexOfPosition(FACE, it.getBound().rankObj.ithObj, j);
DegreeOfFreedom dof = localIndices[dofNo];
if (dofOwner[dof] > it.getRank())
it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
else
rankBoundaryDofs.insert(dof);
}
} }
// === Create the new interior boundaries consisting only of edges. This === // === Create the new interior boundaries consisting only of edges. This ===
...@@ -1486,16 +1470,16 @@ namespace AMDiS { ...@@ -1486,16 +1470,16 @@ namespace AMDiS {