diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index aa662b5e41d01164eb6b4ed050e9395376b63c2b..5c7dc0bc214b65493dedeb18b0d8ad5985188096 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -185,7 +185,7 @@ int main(int argc, char *argv[]) {
const auto & coarseContactNetwork = *contactNetwork.level(0);
const auto & fineContactNetwork = *contactNetwork.level(1);
SupportPatchFactory<typename ContactNetwork::LevelContactNetwork> supportPatchFactory(coarseContactNetwork, fineContactNetwork);
- size_t bodyID = 1;
+ size_t bodyID = 2;
size_t patchDepth = 0;
std::cout << std::endl;
diff --git a/src/spatial-solving/preconditioners/supportpatchfactory.hh b/src/spatial-solving/preconditioners/supportpatchfactory.hh
index a7274816d1b90f26900b9ca9e90fafa6ef1a6760..a4a91be78eca6725613289d0b81b9d05c6217abf 100644
--- a/src/spatial-solving/preconditioners/supportpatchfactory.hh
+++ b/src/spatial-solving/preconditioners/supportpatchfactory.hh
@@ -293,6 +293,15 @@ class SupportPatchFactory
return false;
};
+ auto isRPatchIntersection = [this](std::vector<std::set<size_t>>& dofs, const std::set<size_t>& patchVertices) {
+ for (auto& dof : dofs[0]) {
+ if (patchVertices.count(dof)) {
+ return true;
+ }
+ }
+ return false;
+ };
+
// construct coarse patch
std::vector<BodyElement> patchElements;
patchElements.emplace_back(seedElement);
@@ -346,19 +355,28 @@ class SupportPatchFactory
if (coarseIntersectionMapper_.count(faceIdx)) {
const auto& intersections = coarseIntersectionMapper_[faceIdx];
for (size_t i=0; i<intersections.size(); i++) {
+ const auto& intersection = intersections[i];
+ const auto& rIs = intersection.get(coarseContactNetwork_);
+
BodyElement neighbor;
- neighbor.set(coarseContactNetwork_, intersections[i]);
+ neighbor.set(coarseContactNetwork_, intersection);
size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);
if (visitedElements.count(neighborIdx))
continue;
- std::cout << "elementID: " << coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element) << std::endl;
+ std::vector<std::set<size_t>> dofs;
+ remoteIntersectionDofs(coarseIndices_, rIs, bodyIdx, neighbor.bodyID, dofs, intersection.flip);
- patchVertices.insert(neighborFaceDofs_[faceIdx].begin(), neighborFaceDofs_[faceIdx].end());
- patchElements.emplace_back(neighbor);
- patchSeeds.push(neighbor);
+ if (isRPatchIntersection(dofs, patchVertices)) {
+ patchVertices.insert(dofs[1].begin(), dofs[1].end());
+ patchElements.emplace_back(neighbor);
+ patchSeeds.push(neighbor);
+ } else {
+ newPatchVertices.insert(dofs[1].begin(), dofs[1].end());
+ nextSeeds.emplace_back(neighbor);
+ }
print(patchVertices, "patchVertices: ");
}
@@ -418,6 +436,8 @@ class SupportPatchFactory
if (visited.count(elemIdx))
continue;
+ std::cout << "coarse element: " << elemIdx << std::endl;
+
visited.insert(elemIdx);
const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID)->gridView().grid();
@@ -505,11 +525,7 @@ class SupportPatchFactory
const size_t outsideID = intersection.bodyID;
std::vector<std::set<size_t>> rIsDofs;
- if (intersection.flip) {
- remoteIntersectionDofs(fineIndices_, rIs, outsideID, bodyID, rIsDofs);
- } else {
- remoteIntersectionDofs(fineIndices_, rIs, bodyID, outsideID, rIsDofs);
- }
+ remoteIntersectionDofs(fineIndices_, rIs, bodyID, outsideID, rIsDofs, intersection.flip);
if (rIs.neighbor()) {
Element outsideFather;
@@ -566,7 +582,9 @@ class SupportPatchFactory
}
template <class RIntersection>
- void remoteIntersectionDofs(const NetworkIndexMapper& indices, const RIntersection& is, const size_t insideID, const size_t outsideID, std::vector<std::set<size_t>>& dofs) {
+ void remoteIntersectionDofs(const NetworkIndexMapper& indices, const RIntersection& is, const size_t insideID, const size_t outsideID, std::vector<std::set<size_t>>& dofs, bool flip = false) {
+ assert(!flip || (flip && is.neighbor()));
+
dofs.resize(2);
dofs[0].clear();
dofs[1].clear();
@@ -574,28 +592,59 @@ class SupportPatchFactory
const auto& isGeo = is.geometry();
const auto& isRefElement = Dune::ReferenceElements<ctype, dim-1>::general(is.type());
- const auto& inside = is.inside();
- const auto& insideGeo = inside.geometry();
- const auto& insideRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
- const int dimElement = Dune::ReferenceElements<ctype, dim>::dimension;
+ if (!flip) {
+ const auto& inside = is.inside();
+ const auto& insideGeo = inside.geometry();
+ const auto& insideRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+ const int dimElement = Dune::ReferenceElements<ctype, dim>::dimension;
+
+ for (int i=0; i<insideRefElement.size(is.indexInInside(), 1, dimElement); i++) {
+ size_t idxInElement = insideRefElement.subEntity(is.indexInInside(), 1, i, dimElement);
+
+ const auto& localCorner = isGeo.local(insideGeo.corner(idxInElement));
+
+ if (isRefElement.checkInside(localCorner)) {
+ dofs[0].insert(indices.vertexIndex(insideID, inside, idxInElement));
+ }
+ }
+
+ if (is.neighbor()) {
+ const auto& outside = is.outside();
+ const auto& outsideGeo = outside.geometry();
+ const auto& outsideRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
+
+ for (int i=0; i<outsideRefElement.size(is.indexInOutside(), 1, dimElement); i++) {
+ size_t idxInElement = outsideRefElement.subEntity(is.indexInOutside(), 1, i, dimElement);
- for (int i=0; i<insideRefElement.size(is.indexInInside(), 1, dimElement); i++) {
- size_t idxInElement = insideRefElement.subEntity(is.indexInInside(), 1, i, dimElement);
+ const auto& localCorner = isGeo.local(outsideGeo.corner(idxInElement));
- const auto& localCorner = isGeo.local(insideGeo.corner(idxInElement));
+ if (isRefElement.checkInside(localCorner)) {
+ dofs[1].insert(indices.vertexIndex(outsideID, outside, idxInElement));
+ }
+ }
+ }
+ } else {
+ const auto& inside = is.outside();
+ const auto& insideGeo = inside.geometry();
+ const auto& insideRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+ const int dimElement = Dune::ReferenceElements<ctype, dim>::dimension;
- if (isRefElement.checkInside(localCorner)) {
- dofs[0].insert(indices.vertexIndex(insideID, inside, idxInElement));
+ for (int i=0; i<insideRefElement.size(is.indexInOutside(), 1, dimElement); i++) {
+ size_t idxInElement = insideRefElement.subEntity(is.indexInOutside(), 1, i, dimElement);
+
+ const auto& localCorner = isGeo.local(insideGeo.corner(idxInElement));
+
+ if (isRefElement.checkInside(localCorner)) {
+ dofs[0].insert(indices.vertexIndex(insideID, inside, idxInElement));
+ }
}
- }
- if (is.neighbor()) {
- const auto& outside = is.outside();
+ const auto& outside = is.inside();
const auto& outsideGeo = outside.geometry();
const auto& outsideRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
- for (int i=0; i<outsideRefElement.size(is.indexInOutside(), 1, dimElement); i++) {
- size_t idxInElement = outsideRefElement.subEntity(is.indexInOutside(), 1, i, dimElement);
+ for (int i=0; i<outsideRefElement.size(is.indexInInside(), 1, dimElement); i++) {
+ size_t idxInElement = outsideRefElement.subEntity(is.indexInInside(), 1, i, dimElement);
const auto& localCorner = isGeo.local(outsideGeo.corner(idxInElement));
@@ -604,6 +653,7 @@ class SupportPatchFactory
}
}
}
+
}
template <class Intersection>