.

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;

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>