support for vertices inside element

dune/iga/geometrykernel/slicecurve.hh

@@ -67,9 +67,9 @@ namespace Dune::IGANEW {
     static_assert(GeoCurve::mydimension == 1);
 
     auto domain = geoCurve.domain();
-    if (domain[0].front() == t[0] and domain[0].back() == t[1]) return geoCurve;
-    if (domain[0].front() == t[0]) return std::get<0>(splitCurve(geoCurve, t[1]));
-    if (domain[0].back() == t[1]) return std::get<1>(splitCurve(geoCurve, t[0]));
+    if (FloatCmp::eq(domain[0].front(), t[0]) and FloatCmp::eq(domain[0].back(), t[1])) return geoCurve;
+    if (FloatCmp::eq(domain[0].front(),t[0])) return std::get<0>(splitCurve(geoCurve, t[1]));
+    if (FloatCmp::eq(domain[0].back(),t[1])) return std::get<1>(splitCurve(geoCurve, t[0]));
 
     auto tmpCurve = std::get<1>(splitCurve(geoCurve, t[0]));
     return std::get<0>(splitCurve(tmpCurve, t[1]));

dune/iga/test/testibrareader.cpp

@@ -26,12 +26,12 @@ auto testIbraReader() {
   using GridFactory = Dune::GridFactory<PatchGrid>;
 
   auto gridFactory = GridFactory();
-  gridFactory.insertTrimParameters(GridFactory::TrimParameterType{120});
+  gridFactory.insertTrimParameters(GridFactory::TrimParameterType{100});
 
   // gridFactory.insertJson("auxiliaryfiles/element_trim.ibra", true, {1, 1});
-  gridFactory.insertJson("auxiliaryfiles/trim_2edges.ibra", true, {0, 0});
+  gridFactory.insertJson("auxiliaryfiles/trim_2edges.ibra", true, {2, 2});
   // gridFactory.insertJson("auxiliaryfiles/trim_multi.ibra", true, {1, 1});
-  gridFactory.insertJson("auxiliaryfiles/element_trim_xb.ibra", true, {1, 1});
+  // gridFactory.insertJson("auxiliaryfiles/element_trim_xb.ibra", true, {1, 1});
   // gridFactory.insertJson("auxiliaryfiles/pipe_trim.ibra", true, {1, 1});
 
   auto grid = gridFactory.createGrid();

dune/iga/trimmer/defaulttrimmer/createlevel.hh

@@ -158,6 +158,9 @@ namespace Dune::IGANEW::DefaultTrim {
       int edgeIndex                   = 0;
       int vertexIndex                 = 0;
 
+
+      auto figure = matplot::figure(true);
+
       for (const auto& ele : elements(gv)) {
         ElementTrimFlag eleTrimFlag{ElementTrimFlag::full};
         ElementTrimData eleTrimData{eleTrimFlag};
@@ -167,7 +170,7 @@ namespace Dune::IGANEW::DefaultTrim {
 
           // Testing Purposes
           auto resName = "ele_" + std::to_string(unTrimmedElementIndex + trimmedElementIndex);
-          eleTrimData.drawResult(resName, ele.geometry());
+          eleTrimData.drawResult(resName, ele.geometry(), false);
         }
 
         auto hostId = globalIdSetParameterSpace.id(ele);
@@ -248,6 +251,7 @@ namespace Dune::IGANEW::DefaultTrim {
           }
         }
       }
+      matplot::save("grid", "gif");
 
       // save numbers of untrimmed and trimmed elements per level
       entityContainer.numberOfTrimmedElements.push_back(trimmedElementIndex);

dune/iga/trimmer/defaulttrimmer/elementtrimdata.hh

@@ -63,17 +63,21 @@ namespace Dune::IGANEW::DefaultTrim {
     }
 
     /* this is for testing purposes only */
-    void drawResult(const std::string& filename, auto eleGeometry) {
-      if (static_cast<int>(flag_) != 2) return;
+    void drawResult(const std::string& filename, auto eleGeometry, bool newFig = true) {
+      if (static_cast<int>(flag_) == 1) return;
+
+      if (newFig)
+        auto figure = matplot::figure(true);
 
-      auto figure = matplot::figure(true);
       matplot::hold("on");
 
       constexpr std::array<std::array<int, 2>, 4> edgeLookUp{std::array{0, 1}, {1, 3}, {3, 2}, {2, 0}};
       constexpr std::array idxLookUp = {0, 1, 3, 2};
 
-      auto plotEllipse = [](const Vertex& v) {
-        constexpr auto w = 0.025;
+      const double scale = eleGeometry.volume() / 5;
+
+      auto plotEllipse = [&](const Vertex& v) {
+        const auto w = 0.025 * scale;
         const auto c     = matplot::ellipse(v[0] - (w / 2), v[1] - (w / 2), w, w);
         c->color("blue");
       };
@@ -116,7 +120,8 @@ namespace Dune::IGANEW::DefaultTrim {
       }
 
       matplot::axis(matplot::equal);
-      matplot::save(filename, "gif");
+      if (newFig)
+        matplot::save(filename, "gif");
     }
 
     // Getter

dune/iga/trimmer/defaulttrimmer/patchtrimdata.hh

@@ -8,6 +8,7 @@ namespace Dune::IGANEW::DefaultTrim {
   struct PatchTrimDataImpl {
     using TrimmingCurve = typename GridImp::GridFamily::TrimmerTraits::TrimmingCurve;
     using ParameterType = typename GridImp::GridFamily::Trimmer::ParameterType;
+    using ctype = typename GridImp::GridFamily::ctype;
 
     struct BoundaryLoop {
       void insertTrimCurve(const TrimmingCurve& curve) { curves_.push_back(curve); }
@@ -18,6 +19,12 @@ namespace Dune::IGANEW::DefaultTrim {
       std::vector<TrimmingCurve> curves_;
     };
 
+    struct PointInPatch {
+      FieldVector<ctype, 2> pt;
+      size_t curveIdxI;
+      size_t curveIdxJ;
+    };
+
     class CurveManager {
       friend PatchTrimDataImpl;
 
@@ -88,15 +95,38 @@ namespace Dune::IGANEW::DefaultTrim {
     }
     auto getSplitter() const -> typename CurveManager::idx_t { return manager_.splitter_; }
 
-    template <typename ParameterSpaceGrid>
-    void prepare(const ParameterType& parameters, const std::unique_ptr<ParameterSpaceGrid>& parameterSpaceGrid) {
+    auto getPointsInPatch(size_t loopIndex) const -> const std::vector<PointInPatch>& {
+      return pointsInPatch_.at(loopIndex);
+    }
+
+    void prepare(const ParameterType& parameters, const std::array<std::vector<ctype>, 2>& coordinates) {
       manager_.splitter_ = parameters.splitter;
       std::ranges::for_each(loops_, [&](const auto& loop) { manager_.addLoop(loop); });
 
+      std::array domainU{coordinates[0].front(), coordinates[0].back()};
+      std::array domainV{coordinates[1].front(), coordinates[1].back()};
+
+      // Determine points of the outer boundary that are not on the edges
+      auto isInsidePatch = [&](const FieldVector<ctype, 2>& pt) {
+        const bool isInsideU = FloatCmp::gt(pt[0], domainU.front()) and FloatCmp::lt(pt[0], domainU.back());
+        const bool isInsideV = FloatCmp::gt(pt[1], domainV.front()) and FloatCmp::lt(pt[1], domainV.back());
+
+        return isInsideU and isInsideV;
+      };
+
+      pointsInPatch_.push_back({});
+      for (size_t i = 0; const auto& curve : loops_.front().curves()) {
+        if (const auto pt = curve.corner(1); isInsidePatch(pt)) {
+          pointsInPatch_.front().emplace_back(pt, i, (i + 1) % loops_.front().curves().size());
+        }
+        ++i;
+      }
+
       finished_ = true;
     }
 
   private:
+    std::vector<std::vector<PointInPatch>> pointsInPatch_;
     bool finished_ = false;
     std::vector<BoundaryLoop> loops_;
     CurveManager manager_;

dune/iga/trimmer/defaulttrimmer/trimelement.hh

@@ -29,14 +29,19 @@ namespace Dune::IGANEW::DefaultTrim {
     if (flag != ElementTrimFlag::trimmed) return elementTrimData;
 
     /*
-     * Visitors to get  vertex properties
+     * Visitors to get vertex properties
      */
 
     auto nextEntity  = [&](const int i) { return (i + 1) % result.vertices_.size(); };
     auto isNewVertex = [](const auto& vV) { return std::holds_alternative<Util::ClippingResult::NewVertex>(vV); };
+    auto isHostVertex = [](const auto& vV) { return std::holds_alternative<Util::ClippingResult::HostVertex>(vV); };
+    auto isInsideVertex = [](const auto& vV) { return std::holds_alternative<Util::ClippingResult::InsideVertex>(vV); };
     auto getPt       = [](auto&& vV) { return vV.pt; };
     auto getHostIdx  = [](const auto& vV) { return std::get<Util::ClippingResult::HostVertex>(vV).hostIdx; };
     auto getEdgeIdx  = [](const auto& vV) { return std::get<Util::ClippingResult::NewVertex>(vV).onEdgeIdx; };
+    auto getCurveI = [](const auto& vV) {return std::get<Util::ClippingResult::InsideVertex>(vV).curveIdxI; };
+    auto getCurveJ = [](const auto& vV) {return std::get<Util::ClippingResult::InsideVertex>(vV).curveIdxJ; };
+    auto getLoopIdx = [](const auto& vV) {return std::get<Util::ClippingResult::InsideVertex>(vV).loopIdx; };
     auto getTrimmingCurveIdx = [&](auto& vV) -> std::pair<size_t, size_t> {
       return patchTrimData.getIndices(std::get<Util::ClippingResult::NewVertex>(vV).trimmingCurveZ);
     };
@@ -45,6 +50,9 @@ namespace Dune::IGANEW::DefaultTrim {
     // Here the actual code is starting
     ///
 
+    // Major todo: Create fallback to straight line if alog is not able to find correct Trimming CurveIdx
+    // Also todo: Make a second algo that just connects the points as lines (as template?)
+
     // State
     std::vector<FieldVector<ScalarType, dim>> foundVertices;
     bool isOnNewEdge          = false;
@@ -59,14 +67,14 @@ namespace Dune::IGANEW::DefaultTrim {
       auto pt2 = std::visit(getPt, vV2);
 
       // First case edge is completly untrimmed
-      if (!isNewVertex(vV1) and !isNewVertex(vV2)) {
+      if (isHostVertex(vV1) and isHostVertex(vV2)) {
         elementTrimData.addEdge(Util::giveEdgeIdx(getHostIdx(vV1), getHostIdx(vV2)));
         foundVertices.push_back({pt2.x, pt2.y});
         continue;
       }
 
       // Second case edge begins on a hostVertes and ends on a newVertex
-      if (!isNewVertex(vV1) and isNewVertex(vV2)) {
+      if (isHostVertex(vV1) and isNewVertex(vV2)) {
         currentCurveIdx = getTrimmingCurveIdx(vV2);
         auto [tParam, curvePoint]
             = Util::callFindIntersection(patchTrimData.getCurve(currentCurveIdx), getEdgeIdx(vV2), pt2, corners);
@@ -84,7 +92,7 @@ namespace Dune::IGANEW::DefaultTrim {
         continue;
       }
       // Third case newVertex - newVertex
-      if (isNewVertex(vV1), isNewVertex(vV2)) {
+      if (isNewVertex(vV1) and isNewVertex(vV2)) {
         if (foundVertices.empty()) {
           currentCurveIdx = getTrimmingCurveIdx(vV1);
           FieldVector<ScalarType, dim> curvePoint;
@@ -119,7 +127,7 @@ namespace Dune::IGANEW::DefaultTrim {
         continue;
       }
       // Fourth case edge begins on a newVertex and ends in a HostVertex
-      if (isNewVertex(vV1), !isNewVertex(vV2)) {
+      if (isNewVertex(vV1), isHostVertex(vV2)) {
         auto v2 = Dune::FieldVector<ScalarType, dim>{pt2.x, pt2.y};
 
         FieldVector<ScalarType, dim> p;
@@ -135,6 +143,38 @@ namespace Dune::IGANEW::DefaultTrim {
         foundVertices.push_back({pt2.x, pt2.y});
         continue;
       }
+      // Additional cases to cover inside vertices
+      if (isNewVertex(vV1), isInsideVertex(vV2)) {
+        if (foundVertices.empty()) {
+          currentCurveIdx = getTrimmingCurveIdx(vV1);
+          FieldVector<ScalarType, dim> curvePoint;
+          std::tie(currentT, curvePoint)
+              = Util::callFindIntersection(patchTrimData.getCurve(currentCurveIdx), getEdgeIdx(vV1), pt1, corners);
+          std::cout << "Found: " << curvePoint << " From Clipping: " << pt1.x << " " << pt1.y << " t: " << currentT
+                    << std::endl;
+          foundVertices.push_back(curvePoint);
+        }
+        const auto& curve = patchTrimData.loops()[getLoopIdx(vV2)].curves()[getCurveI(vV2)];
+        double tParam = curve.domain().front().back();
+        auto elementTrimmingCurve = Util::createTrimmingCurveSlice(curve, currentT, tParam);
+        auto curvePoint = curve.corner(1);
+        elementTrimData.addEdgeNewNew(elementTrimmingCurve, curvePoint);
+        foundVertices.push_back(curvePoint);
+        continue;
+      }
+      if (isInsideVertex(vV1), isNewVertex(vV2)) {
+        const auto& curve = patchTrimData.loops()[getLoopIdx(vV1)].curves()[getCurveJ(vV1)];
+        currentT = curve.domain().front().front();
+        auto [tParam, curvePoint]
+            = Util::callFindIntersection(curve, getEdgeIdx(vV2), pt2, corners);
+        std::cout << "Found: " << curvePoint << " From Clipping: " << pt2.x << " " << pt2.y << " t: " << tParam
+                  << std::endl;
+        auto elementTrimmingCurve
+              = Util::createTrimmingCurveSlice(curve, currentT, tParam);
+        elementTrimData.addEdgeNewNew(elementTrimmingCurve, curvePoint);
+        foundVertices.push_back(curvePoint);
+        continue;
+      }
     }
 
     return elementTrimData;

dune/iga/trimmer/defaulttrimmer/trimmer.hh

@@ -498,7 +498,8 @@ namespace Dune::IGANEW {
       GridImp* grid_;
 
       void setup() {
-        if (trimData_.has_value()) trimData_->prepare(parameters_, untrimmedParameterSpaceGrid_);
+        if (trimData_.has_value())
+          trimData_->prepare(parameters_, grid_->tensorProductCoordinates(0));
       }
       /**
        * @brief Change the parameters to the trimmer.

dune/iga/trimmer/defaulttrimmer/trimmingutils/clipelementrectangle.hh

@@ -87,6 +87,13 @@ namespace Dune::IGANEW::DefaultTrim::Util {
       if (auto [isHost, idx] = isHostVertex(pt, eleRect); isHost) result.addOriginalVertex(idx);
     }
 
+    // Now we have to check for points inside the elements
+    for (const auto& pt : patchTrimData.getPointsInPatch(0)) {
+      if (const PointD cp{pt.pt[0], pt.pt[1]}; PointInPolygon(cp, eleRect) == PointInPolygonResult::IsInside) {
+        result.addInsideVertex(cp, pt.curveIdxI, pt.curveIdxJ, 0);
+      }
+    }
+
     // Add startpoints of trimming curves to the vertices if they are on one of the element edges
     // \todo not exactly sure but i think this only has to check the outer boundary loop
     for (auto cI = 0; cI < patchTrimData.loops().front().curves().size(); ++cI) {

dune/iga/trimmer/defaulttrimmer/trimmingutils/cliputils.hh

@@ -61,8 +61,14 @@ namespace Dune::IGANEW::DefaultTrim::Util {
       Clipper2Lib::PointD pt{};
       size_t trimmingCurveZ{};
     };
+    struct InsideVertex {
+      Clipper2Lib::PointD pt{};
+      size_t curveIdxI{};
+      size_t curveIdxJ{};
+      size_t loopIdx{};
+    };
 
-    using VertexVariant = std::variant<HostVertex, NewVertex>;
+    using VertexVariant = std::variant<HostVertex, NewVertex, InsideVertex>;
 
     auto isAlreadyThere(const auto& pt) {
       const auto it = std::ranges::find_if(vertices_, [&](const VertexVariant& vertexVariant) {
@@ -95,6 +101,9 @@ namespace Dune::IGANEW::DefaultTrim::Util {
     void addNewVertex(const int edgeIdx, const Clipper2Lib::PointD& pt, const size_t trimmingCurveZ) {
       if (not isAlreadyThere(pt)) vertices_.emplace_back(NewVertex(edgeIdx, pt, trimmingCurveZ));
     }
+    void addInsideVertex(const Clipper2Lib::PointD& pt, const size_t curveIndexI, const size_t curveIndexJ, size_t loopIdx) {
+      vertices_.emplace_back(InsideVertex(pt, curveIndexI, curveIndexJ, loopIdx));
+    }
 
     void finish() {
       // Sort the points in counter clockwise manner such that the first point is in the lower left
@@ -137,18 +146,20 @@ namespace Dune::IGANEW::DefaultTrim::Util {
         std::ranges::rotate(vertices_, it);
       else
         std::cout << "Warning, no HostVertex" << std::endl;
-      // DUNE_THROW(Dune::NotImplemented, "Algorithm needs at least one HostVertex to work");
     }
 
     void report() const {
       std::cout << "Vertices found\n";
       struct Visitor {
-        void operator()(const ClippingResult::HostVertex& v) const {
+        void operator()(const HostVertex& v) const {
           std::cout << "Pt: " << v.pt << " Host Idx: " << v.hostIdx << std::endl;
         }
-        void operator()(const ClippingResult::NewVertex& v) const {
+        void operator()(const NewVertex& v) const {
           std::cout << "Edge: " << v.onEdgeIdx << " Pt: " << v.pt << " On TC: " << v.trimmingCurveZ << std::endl;
         }
+        void operator()(const InsideVertex& v) const {
+          std::cout << " Pt: " << v.pt << " On TC: " << v.curveIdxI << " and " << v.curveIdxJ << std::endl;
+        }
       };
       for (auto& vV : vertices_)
         std::visit(Visitor{}, vV);