add isConnectedAtBoundary function to NurbsPatchGeometry class

dune/iga/geometrykernel/nurbspatchgeometry.hh

@@ -317,6 +317,102 @@ namespace Dune::IGANEW::GeometryKernel {
     /* @brief Get the patch data of the NURBS patch. */
     auto& patchData() { return patchData_; }
 
+    struct ConnectionResult {
+      bool isConnected;
+      int boundary;
+      operator bool() const { return isConnected; }
+    };
+
+    /**
+     * @brief Check if the patch is connected at the given boundary.
+     *
+     * @details This function checks if the patch is connected at the specified boundary. The boundary index is
+     * determined by the vertex numbering for lines, the edge numbering for quadrilaterals, and the face numbering for
+     * cubes. For more details, refer to DUNE Book Chapter 5.5.
+     *
+     * @param boundary The index of the boundary to check (vertex for lines, edge for quadrilaterals, face for cubes).
+     * @param tol Absolute tolerance for comparing control points. Default is 1e-8.
+     *
+     * @return A ConnectionResult object indicating whether the patch is connected at the specified boundary.
+     *         The isConnected member is true if the patch is connected, and the boundary member is the index of the
+     * adjacent boundary if connected, or -1 if not connected.
+     *
+     * @tparam boundary The boundary index to check.
+     * @tparam tol Absolute tolerance for comparing control points. Default is 1e-8.
+     *
+     * @throws Dune::NotImplemented If the dimension is greater than 3 (higher dimensions not implemented yet).
+     */
+    ConnectionResult isConnectedAtBoundary(int boundary, double tol = 1e-8) const {
+      const auto& cps = patchData_.controlPoints;
+      if constexpr (mydimension == 1) {
+        if (auto cp = cps.directGetAll();
+            FloatCmp::eq<FieldVector<ctype, coorddimension>, FloatCmp::CmpStyle::absolute>(cp.front().p, cp.back().p,
+                                                                                           tol))
+          return ConnectionResult(true, !boundary);
+        else
+          return ConnectionResult(false, -1);
+      } else if constexpr (mydimension == 2) {
+        auto edgeRange = [&cps](int b) {
+          switch (b) {
+            case 0:
+              return cps.hyperSurfFront({1});
+            case 1:
+              return cps.hyperSurfBack({1});
+            case 2:
+              return cps.hyperSurfFront({0});
+            case 3:
+              return cps.hyperSurfBack({0});
+            default:
+              assert("No valid boundary index");
+          }
+          __builtin_unreachable();
+        };
+
+        auto theEdgeToCheck = edgeRange(boundary);
+        for (auto i : Dune::range(4)) {
+          if (boundary == i) continue;
+          if (std::ranges::equal(theEdgeToCheck, edgeRange(i), [tol](auto f, auto s) {
+                return FloatCmp::eq<FieldVector<ctype, coorddimension>, FloatCmp::CmpStyle::absolute>(f.p, s.p, tol);
+              }))
+            return ConnectionResult(true, i);
+        }
+        return ConnectionResult(false, -1);
+      } else if constexpr (mydimension == 3) {
+        auto edgeRange = [&cps](int b) {
+          switch (b) {
+            case 0:
+              return cps.hyperSurfFront({1, 2});
+            case 1:
+              return cps.hyperSurfBack({1, 2});
+            case 2:
+              return cps.hyperSurfFront({0, 1});
+            case 3:
+              return cps.hyperSurfBack({0, 1});
+            case 4:
+              return cps.hyperSurfFront({0, 2});
+            case 5:
+              return cps.hyperSurfBack({0, 2});
+            default:
+              assert("No valid boundary index");
+          }
+          __builtin_unreachable();
+        };
+
+        auto theEdgeToCheck = edgeRange(boundary);
+        for (auto i : Dune::range(6)) {
+          if (boundary == i) continue;
+          if (std::ranges::equal(theEdgeToCheck, edgeRange(i), [tol](auto f, auto s) {
+                return FloatCmp::eq<FieldVector<ctype, coorddimension>, FloatCmp::CmpStyle::absolute>(f.p, s.p, tol);
+              }))
+            return ConnectionResult(true, i);
+        }
+        return ConnectionResult(false, -1);
+      } else
+        DUNE_THROW(Dune::NotImplemented,
+                   "Higher dimensions not implemented yet. But should be straight-forward to do so, somebody just has "
+                   "to come up with a resenable indexing of hyper surfaces.");
+    }
+
     friend auto referenceElement(const NURBSPatch& geo) {
       return referenceElement<double, mydimension>(Dune::GeometryTypes::cube(mydimension));
     }

dune/iga/hierarchicpatch/patchgrid.hh

@@ -456,6 +456,7 @@ namespace Dune::IGANEW {
     auto untrimmedElementNumbers(int lvl) const { return patchGeometries_[lvl].numberOfSpans(); }
 
     const auto& trimmer() const { return *trimmer_; }
+    const auto& patchGeometry(int i) const { return patchGeometries_.at(i); }
 
   private:
     PatchGrid() = default;

dune/iga/test/gridtests.cpp

@@ -184,6 +184,21 @@ auto testNurbsGridCylinder() {
 
   TestSuite testSuite;
 
+  auto patchGeo = grid.patchGeometry(0);
+
+  // Since this is only a quarter of a cylinder no glueing should be present
+  auto [glued, edgeIndex] = patchGeo.isConnectedAtBoundary(0);
+  testSuite.check(!glued) << "The edge " << 0 << " should be glued!";
+  testSuite.check(edgeIndex == -1);
+  auto [glued2, edgeIndex2] = patchGeo.isConnectedAtBoundary(0);
+  testSuite.check(!glued2) << "The edge " << 2 << " shouldn't be glued!";
+  testSuite.check(edgeIndex == -1);
+
+  Dune::RefinementIntervals refinementIntervals1(3);
+  SubsamplingVTKWriter<decltype(grid.leafGridView())> vtkWriter(grid.leafGridView(), refinementIntervals1);
+  //  vtkWriter.write("NURBSGridTest-CurveNewFineResample");
+  vtkWriter.write("NURBSGridTest-Zylinder");
+
   testSuite.subTest(thoroughGridCheck(grid));
   return testSuite;
 }
@@ -206,6 +221,33 @@ auto testHierarchicPatch() {
 
   Dune::IGANEW::PatchGrid<2, 3, GridFamily> patch(nurbsPatchData);
 
+  auto glueIndicator = [](int i) {
+    switch (i) {
+      case 0:
+        return 1;
+      case 1:
+        return 0;
+      case 2:
+        return 3;
+      case 3:
+        return 2;
+      default:
+        DUNE_THROW(RangeError, "Invalid index");
+    }
+  };
+  auto patchGeo = patch.patchGeometry(0);
+  for (auto i : Dune::range(4)) {
+    auto [glued, edgeIndex] = patchGeo.isConnectedAtBoundary(i);
+    t.check(glued) << "The edge " << i << " should be glued!";
+    t.check(edgeIndex == glueIndicator(i))
+        << "The edge " << i << " should be glued to edge " << glueIndicator(i) << " but is glued to " << edgeIndex;
+  }
+
+  Dune::RefinementIntervals refinementIntervals1(3);
+  SubsamplingVTKWriter<decltype(patch.leafGridView())> vtkWriter(patch.leafGridView(), refinementIntervals1);
+  //  vtkWriter.write("NURBSGridTest-CurveNewFineResample");
+  vtkWriter.write("NURBSGridTest-Torus");
+
   t.subTest(thoroughGridCheck(patch));
 
   return t;
@@ -326,6 +368,13 @@ auto testNURBSGridCurve() {
 
   TestSuite t;
   IGANEW::PatchGrid<dim, dimworld, GridFamily> grid(patchData);
+
+  auto patchGeo = grid.patchGeometry(0);
+
+  auto [glued, edgeIndex] = patchGeo.isConnectedAtBoundary(0);
+  t.check(glued) << "The edge " << 0 << " should be glued!";
+  t.check(edgeIndex == -1);
+
   grid.globalRefine(3);
   auto gridView        = grid.leafGridView();
   const auto& indexSet = gridView.indexSet();
@@ -389,6 +438,13 @@ auto testNURBSGridSurface() {
   auto gridView        = grid.leafGridView();
   const auto& indexSet = gridView.indexSet();
 
+  auto patchGeo = grid.patchGeometry(0);
+  for (auto i : Dune::range(4)) {
+    auto [glued, edgeIndex] = patchGeo.isConnectedAtBoundary(i);
+    t.check(!glued) << "The edge " << i << " shouldn't be glued!";
+    t.check(edgeIndex == -1);
+  }
+
   Dune::RefinementIntervals refinementIntervals1(subSampling);
   SubsamplingVTKWriter<decltype(gridView)> vtkWriter(gridView, refinementIntervals1);
   vtkWriter.write("NURBSGridTest-Surface");
@@ -451,6 +507,14 @@ auto test3DGrid() {
 
   auto gridView = grid.leafGridView();
   TestSuite t;
+  auto patchGeo = grid.patchGeometry(0);
+
+  for (auto i : Dune::range(6)) {
+    auto [glued, edgeIndex] = patchGeo.isConnectedAtBoundary(i);
+    t.check(!glued) << "The surface " << i << " shouldn't be glued!";
+    t.check(edgeIndex == -1);
+  }
+
   t.subTest(checkUniqueEdges(gridView));
   t.subTest(checkUniqueSurfaces(gridView));
 

dune/iga/trimmer/defaulttrimmer/trimelement.hh

@@ -98,32 +98,31 @@ namespace Dune::IGANEW::DefaultTrim {
             = Util::callFindIntersection(patchTrimData.getCurve(currentCurveIdx), getEdgeIdx(vV2), pt2, corners);
         assertPoint(pt2, curvePoint);
 
-        // If currentT > tParam it can have 2 reasons. 1) the trim id only on one side of the element, and it has to connect back
-        // or 2) if the trimming curve consist only of one curve where it has the same front and back Controlpoint, so sometimes
-        // the wrong tParam (front or back) gets determined
+        // If currentT > tParam it can have 2 reasons. 1) the trim id only on one side of the element, and it has to
+        // connect back or 2) if the trimming curve consist only of one curve where it has the same front and back
+        // Controlpoint, so sometimes the wrong tParam (front or back) gets determined
         if (currentT > tParam) {
           bool success = false;
           if (currentCurveIdx.first > 0 && patchTrimData.loops()[currentCurveIdx.first].size() == 1) {
             auto curve = patchTrimData.getCurve(currentCurveIdx);
-            if (auto cp = curve.patchData().controlPoints.directGetAll(); FloatCmp::eq(cp.front().p, cp.back().p)) {
+            if (curve.isConnectedAtBoundary(0)) {
               auto elementTrimmingCurve = Util::createTrimmingCurveSlice(curve, currentT, curve.domain()[0].back());
               elementTrimData.addEdgeNewNew(elementTrimmingCurve, curvePoint);
-              currentT    = std::numeric_limits<double>::infinity();
-              success = true;
+              currentT = std::numeric_limits<double>::infinity();
+              success  = true;
             }
           } else if (getEdgeIdx(vV1) == getEdgeIdx(vV2)) {
             auto trimmedEdge = Util::createHostGeometry<TrimmingCurve>(foundVertices.back(), curvePoint);
             elementTrimData.addEdgeNewNewOnHost(getEdgeIdx(vV2), trimmedEdge, curvePoint);
-            currentT    = tParam;
-            success = true;
+            currentT = tParam;
+            success  = true;
           }
-          if (not success)
-            throwGridError();
+          if (not success) throwGridError();
         } else {
           auto elementTrimmingCurve
               = Util::createTrimmingCurveSlice(patchTrimData.getCurve(currentCurveIdx), currentT, tParam);
           elementTrimData.addEdgeNewNew(elementTrimmingCurve, curvePoint);
-          currentT    = std::numeric_limits<double>::infinity();
+          currentT = std::numeric_limits<double>::infinity();
         }
 
         foundVertices.push_back(curvePoint);