diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h index 4cfc027ee3e47e040eb1a65ffb9dca48ce4909e1..2f28c95689b85ed2507409ff06b784d97ab8be9b 100644 --- a/AMDiS/src/Element.h +++ b/AMDiS/src/Element.h @@ -388,13 +388,10 @@ namespace AMDiS { * \param[in] geoPos Must be either EDGE or FACE. Defines whether an * edge or a face (only in 3d) should be traversed. * \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, BoundaryObject bound, - DofContainer& dofs, - bool parentVertices = false) const = 0; + DofContainer& dofs) const = 0; /** \brief * Traverses an edge/face of a given element (this includes also all children of diff --git a/AMDiS/src/Line.h b/AMDiS/src/Line.h index 7698cb6e9bf8631a7ace1d32c118a3d35837d209..e54f3242073c0935d4dbb233493dae5e22298c60 100644 --- a/AMDiS/src/Line.h +++ b/AMDiS/src/Line.h @@ -154,8 +154,7 @@ namespace AMDiS { return "Line"; } - void getVertexDofs(FiniteElemSpace*, BoundaryObject, DofContainer&, - bool parentVertices = false) const + void getVertexDofs(FiniteElemSpace*, BoundaryObject, DofContainer&) const { FUNCNAME("Line::getVertexDofs()"); ERROR_EXIT("Not yet implemented!\n"); diff --git a/AMDiS/src/Tetrahedron.cc b/AMDiS/src/Tetrahedron.cc index 8899c1f2fe343318ade1297d1f959159f8e36932..71eb822de119f10c0fa85cfa17a00dcb62b342ea 100644 --- a/AMDiS/src/Tetrahedron.cc +++ b/AMDiS/src/Tetrahedron.cc @@ -187,8 +187,7 @@ namespace AMDiS { void Tetrahedron::getVertexDofs(FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs, - bool parentVertices) const + DofContainer& dofs) const { FUNCNAME("Tetrahedron::getVertexDofs()"); @@ -197,10 +196,10 @@ namespace AMDiS { dofs.push_back(dof[bound.ithObj]); break; case EDGE: - getVertexDofsAtEdge(feSpace, bound, dofs, parentVertices); + getVertexDofsAtEdge(feSpace, bound, dofs); break; case FACE: - getVertexDofsAtFace(feSpace, bound, dofs, parentVertices); + getVertexDofsAtFace(feSpace, bound, dofs); break; default: ERROR_EXIT("Should not happen!\n"); @@ -210,31 +209,8 @@ namespace AMDiS { void Tetrahedron::getVertexDofsAtFace(FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs, - bool parentVertices) const + DofContainer& dofs) 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]) return; @@ -294,15 +270,10 @@ namespace AMDiS { void Tetrahedron::getVertexDofsAtEdge(FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs, - bool parentVertices) const + DofContainer& dofs) const { FUNCNAME("Tetrahedron::getVertexDofsAtEdge()"); - if (parentVertices) { - ERROR_EXIT("Einbau excludedSubstructures!\n"); - } - if (!child[0]) return; diff --git a/AMDiS/src/Tetrahedron.h b/AMDiS/src/Tetrahedron.h index f20eca12033f91b3e62aaa0e3ccafbd9b1219076..0b312eb46d7fede67b5dd61a306947dc572e8d0a 100644 --- a/AMDiS/src/Tetrahedron.h +++ b/AMDiS/src/Tetrahedron.h @@ -137,18 +137,15 @@ namespace AMDiS { void getVertexDofs(FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs, - bool parentVertices = false) const; + DofContainer& dofs) const; void getVertexDofsAtFace(FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs, - bool parentVertices) const; + DofContainer& dofs) const; void getVertexDofsAtEdge(FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs, - bool parentVertices) const; + DofContainer& dofs) const; void getNonVertexDofs(FiniteElemSpace* feSpace, BoundaryObject bound, diff --git a/AMDiS/src/Triangle.cc b/AMDiS/src/Triangle.cc index 8772ab1ca511a344c617a507d57d026ad05c6009..9ccb8f8ab161d112a2b305e3f0d3432e2beec342 100644 --- a/AMDiS/src/Triangle.cc +++ b/AMDiS/src/Triangle.cc @@ -77,27 +77,10 @@ namespace AMDiS { void Triangle::getVertexDofs(FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs, - bool parentVertices) const + DofContainer& dofs) const { 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; switch (bound.ithObj) { @@ -143,22 +126,6 @@ namespace AMDiS { default: 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"); - } - } } diff --git a/AMDiS/src/Triangle.h b/AMDiS/src/Triangle.h index 8b6a369abf9fd479991f84b1084fb7085e1e3aaa..54d07f681d10e4516bd77b50c3c8d27b73f07730 100644 --- a/AMDiS/src/Triangle.h +++ b/AMDiS/src/Triangle.h @@ -166,8 +166,7 @@ namespace AMDiS { void getVertexDofs(FiniteElemSpace* feSpace, BoundaryObject bound, - DofContainer& dofs, - bool parentVertices = false) const; + DofContainer& dofs) const; void getNonVertexDofs(FiniteElemSpace* feSpace, BoundaryObject bound, diff --git a/AMDiS/src/parallel/GlobalMatrixSolver.cc b/AMDiS/src/parallel/GlobalMatrixSolver.cc index a8370cfb7bc549bb32de20b95939b74fddfec047..2a9b5229ec2b101b298cf28bd6c26932c0eb419d 100644 --- a/AMDiS/src/parallel/GlobalMatrixSolver.cc +++ b/AMDiS/src/parallel/GlobalMatrixSolver.cc @@ -10,7 +10,7 @@ namespace AMDiS { 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; return 0; diff --git a/AMDiS/src/parallel/InteriorBoundary.h b/AMDiS/src/parallel/InteriorBoundary.h index aa19f87598e3a4453c3453ec910625de30a0df8d..59713de7fa52d647a764caa930c5431203e20060 100644 --- a/AMDiS/src/parallel/InteriorBoundary.h +++ b/AMDiS/src/parallel/InteriorBoundary.h @@ -154,6 +154,16 @@ namespace AMDiS { } } + inline AtomicBoundary& operator*() + { + return *vecIt; + } + + inline AtomicBoundary* operator->() + { + return &(*vecIt); + } + void nextRank() { ++mapIt; @@ -168,11 +178,6 @@ namespace AMDiS { return mapIt->first; } - AtomicBoundary& getBound() - { - return *vecIt; - } - protected: inline void nextNonempty() diff --git a/AMDiS/src/parallel/ParallelDomainBase.cc b/AMDiS/src/parallel/ParallelDomainBase.cc index 2cc0f8b7f8e91683e2812002031e2791b297f529..d1e5245e794b5020f32efe881415b9fe5c1e49da 100644 --- a/AMDiS/src/parallel/ParallelDomainBase.cc +++ b/AMDiS/src/parallel/ParallelDomainBase.cc @@ -376,12 +376,12 @@ namespace AMDiS { RankToBoundMap allBound; for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) - if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim)) - allBound[it.getRank()].push_back(it.getBound()); + if (it->rankObj.subObj == INDEX_OF_DIM(dim - 1, dim)) + allBound[it.getRank()].push_back(*it); for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) - if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim)) - allBound[it.getRank()].push_back(it.getBound()); + if (it->rankObj.subObj == INDEX_OF_DIM(dim - 1, dim)) + allBound[it.getRank()].push_back(*it); // === Check the boundaries and adapt mesh if necessary. === @@ -795,17 +795,7 @@ namespace AMDiS { int nNeighbours = mesh->getGeo(NEIGH); int dim = mesh->getDim(); - GeoIndex subObj = CENTER; - switch (dim) { - case 2: - subObj = EDGE; - break; - case 3: - subObj = FACE; - break; - default: - ERROR_EXIT("What is this?\n"); - } + GeoIndex subObj = INDEX_OF_DIM(dim - 1, dim); // === First, traverse the mesh and search for all elements that have an === // === boundary with an element within another partition. === @@ -849,16 +839,9 @@ namespace AMDiS { bound.neighObj.subObj = subObj; 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. int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()]; - // === Add the boundary information object to the corresponding overall === // === boundary. There are three possibilities: if the boundary is a === // === periodic boundary, just add it to \ref periodicBounadry. Here it === @@ -896,24 +879,22 @@ namespace AMDiS { { FUNCNAME("ParallelDomainBase::createSubstructureInteriorBoundaryInfo()"); - // === 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 === - // === the information of macro elements. === + // === Seach for all vertices/edges, which are part of an interior boundary, === + // === but are not part of the interior boundaries that are created based on === + // === the information of macro elements. === int dim = mesh->getDim(); 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. std::map<DegreeOfFreedom, int> dofOwner; - // Maps each DOF in ranks partition of the domain to the element object that - // contains this DOF. + // Maps each DOF in ranks partition to the element object that contains this DOF. std::map<DegreeOfFreedom, BoundaryObject> rankDofs; // Maps each edge in the whole domain to the rank that ownes this edge. std::map<GlobalEdge, int> edgeOwner; - // Maps each edge in ranks partition of the domain to the element object that - // contains this edge. + // Maps each edge in ranks partition to the element object that contains this edge. std::map<GlobalEdge, BoundaryObject> rankEdges; @@ -927,7 +908,19 @@ namespace AMDiS { PartitionElementData *partitionData = dynamic_cast<PartitionElementData*> (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. if (dim == 3) { for (int i = 0; i < 6; i++) { @@ -938,28 +931,18 @@ namespace AMDiS { // Update the owner of the current edge. 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 - // to the set of all rank's edges. + // If the edge is part of an element that is within ranks partition, add it + // to the set of all ranks edges. if (partitionData && partitionData->getPartitionStatus() == IN) 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); } - // === 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. std::set<GlobalEdge> rankBoundaryEdges; @@ -969,12 +952,22 @@ namespace AMDiS { // First, traverse the rank owned elements at the interior boundaries. 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); + 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) { 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 dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)]; GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); @@ -984,51 +977,42 @@ namespace AMDiS { // Otherwise, it is part of the interior boundary and we add it to the set // rankBoundaryEdges. 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 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. 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); + 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) { 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 dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)]; GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); 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 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 === @@ -1486,16 +1470,16 @@ namespace AMDiS { for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) { DofContainer dofs; - it.getBound().rankObj.el->getVertexDofs(feSpace, it.getBound().rankObj, dofs); - it.getBound().rankObj.el->getNonVertexDofs(feSpace, it.getBound().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++) sendDofs[it.getRank()].push_back(dofs[i]); } for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) { DofContainer dofs; - it.getBound().rankObj.el->getNonVertexDofs(feSpace, it.getBound().rankObj, dofs); - it.getBound().rankObj.el->getVertexDofs(feSpace, it.getBound().rankObj, dofs); + it->rankObj.el->getNonVertexDofs(feSpace, it->rankObj, dofs); + it->rankObj.el->getVertexDofs(feSpace, it->rankObj, dofs); for (int i = 0; i < static_cast<int>(dofs.size()); i++) { DofContainer::iterator eraseIt =