Intersections maybe working

dune/iga/geometrykernel/nurbspatchgeometry.hh

@@ -80,6 +80,15 @@ public:
   /* @brief Default constructor for NURBSPatch.*/
   NURBSPatch() = default;
 
+  bool operator==(const NURBSPatch & other) const {
+    bool isSame = true;
+    for (auto i : Dune::range(corners())) {
+      if (not FloatCmp::eq(corner(i), other.corner(i), 1e-8))
+        isSame = false;
+    }
+    return isSame;
+  }
+
   /**
    * @brief Get a local view of the NURBS patch.
    * @tparam codim Codimension of the patch.

dune/iga/hierarchicpatch/patchgrid.hh

@@ -491,6 +491,9 @@ public:
   const auto& patchGeometry(int i) const {
     return patchGeometries_.at(i);
   }
+  const auto& patchGeometryAtBack() const {
+    return patchGeometries_.back();
+  }
 
 private:
   PatchGrid() = default;

dune/iga/hierarchicpatch/patchgridentity.hh

@@ -149,7 +149,7 @@ public:
   }
 
   bool isTrimmed() const {
-    return localEntity_.entityInfo_.trimmed;
+    return localEntity_.entityInfo().trimmed;
   }
 
 private:

dune/iga/hierarchicpatch/patchgridintersections.hh

@@ -125,7 +125,7 @@ namespace Impl {
       // @todo trim this will be wrong as soon as the intersection geometry has a special geoemtry
       auto geo = typename Geometry::Implementation(
           parameterSpaceIntersection_.geometry(),
-          patchGrid_->patchGeometries_.back().template localView<1, GridImp::Trimmer>());
+          patchGridGeometry().template localView<1, typename GridImp::Trimmer>());
       return Geometry(geo);
     }
 
@@ -166,6 +166,12 @@ namespace Impl {
     }
 
   private:
+    auto& patchGridGeometry() const {
+      if constexpr (type_ == IntersectionType::Leaf)
+        return  patchGrid_->patchGeometryAtBack();
+      else
+        return patchGrid_->patchGeometry(inside().level());
+    }
     const GridImp* patchGrid_{};
     ParameterSpaceIntersection parameterSpaceIntersection_{};
   };
@@ -184,7 +190,7 @@ class PatchGridLeafIntersection
 {
   friend typename GridImp::Traits::LeafIntersectionIterator;
 
-  friend struct HostGridAccess<typename std::remove_const<GridImp>::type>;
+  friend struct HostGridAccess<std::remove_const_t<GridImp>>;
 
   using Implementation                 = Impl::PatchGridIntersectionImpl<GridImp, Impl::IntersectionType::Leaf>;
   using ParameterSpaceLeafIntersection = typename Implementation::ParameterSpaceIntersection;

dune/iga/test/gridteststrimmed.cpp

@@ -120,6 +120,7 @@ auto myGridCheck(G& grid) {
   static_assert(G::dimension == 2);
 
   auto testGV = [&]<typename GV>(const GV& gv) {
+    auto& indexSet = gv.indexSet();
     for (int eleIdx = 0; const auto& ele : elements(gv)) {
       std::cout << "Element " << eleIdx << std::endl;
       const int numCorners  = ele.subEntities(2);
@@ -139,6 +140,9 @@ auto myGridCheck(G& grid) {
       // Intersections
       int intersectionCount{0};
       for (const auto& intersection : intersections(gv, ele)) {
+        auto geometry = intersection.geometry();
+        std::cout << "Intersection Idx: " << intersectionCount << ", C: " << geometry.corner(0) << ", "
+                  << geometry.corner(1) << std::endl;
         ++intersectionCount;
       }
       const int numEdges = ele.subEntities(1);
@@ -148,6 +152,12 @@ auto myGridCheck(G& grid) {
 
       t.check(intersectionCount == numEdges) << "There should be as many edges as intersections";
 
+      for (const auto i : Dune::range(numEdges)) {
+        auto edge = ele.template subEntity<1>(i);
+        std::cout << "Edge Idx " << indexSet.index(edge) << ", C: " << edge.geometry().corner(0) << ", "
+                  << edge.geometry().corner(1) << std::endl;
+      }
+
       ++eleIdx;
       std::cout << std::endl;
     }
@@ -200,15 +210,15 @@ auto thoroughGridCheck(auto& grid) {
   }
   t.subTest(gvTest(grid.leafGridView()));
 
-  // gridcheck(grid);
-  t.subTest(myGridCheck(grid));
+  gridcheck(grid);
+  //t.subTest(myGridCheck(grid));
 
-  // try {
-  //   checkIntersectionIterator(grid);
-  //   checkLeafIntersections(grid);
-  // } catch (const Dune::NotImplemented& e) {
-  //   std::cout << e.what() << std::endl;
-  // }
+  try {
+    checkIntersectionIterator(grid);
+    checkLeafIntersections(grid);
+  } catch (const Dune::NotImplemented& e) {
+    std::cout << e.what() << std::endl;
+  }
 
   return t;
 }
@@ -226,7 +236,7 @@ auto testPlate() {
 
   auto gridFactory = GridFactory();
   gridFactory.insertTrimParameters(GridFactory::TrimParameterType{100});
-  gridFactory.insertJson("auxiliaryfiles/element_trim.ibra", true, {1, 1});
+  gridFactory.insertJson("auxiliaryfiles/element_trim.ibra", true, {0, 0});
 
   auto grid = gridFactory.createGrid();
   t.subTest(thoroughGridCheck(*grid));

dune/iga/test/testintersections.cpp

@@ -80,6 +80,7 @@ auto testIntersections(auto& grid, bool trimmed, int refLevel) {
         std::get<1>(resTuple).push_back(intersection.unitOuterNormal({0.5}));
 
         std::cout << "unit outer normal: " << std::setprecision(16) << intersection.unitOuterNormal({0.5}) << std::endl;
+        std::cout << "inside index: " << intersection.indexInInside() << std::endl;
 
         std::get<2>(resTuple).push_back(intersection.centerUnitOuterNormal());
         std::get<3>(resTuple).push_back(intersection.outerNormal({0.5}));

dune/iga/trimmer/defaulttrimmer/entitycontainer.hh

@@ -64,6 +64,32 @@ struct VectorEntityContainer
     __builtin_unreachable();
   }
 
+  size_t subIds(const IdType& elementId, int codim) const {
+    if (codim == 0)
+      return 0;
+    if (codim == 1)
+      return globalEdgesIdOfElementsMap_.at(elementId).size();
+    if (codim == 2)
+      return globalVerticesIdOfElementsMap.at(elementId).size();
+    assert(codim >= 0 and codim <= 2);
+    __builtin_unreachable();
+  }
+
+  int outsideIntersectionIndex(const IdType& insideElementId, const IdType& outsideElementId, int indexInInside) const {
+    const IdType& insideSubId = subId(insideElementId, indexInInside, 1);
+
+    for (const auto i : Dune::range(subIds(outsideElementId, 1))) {
+      const IdType& outsideSubId = subId(outsideElementId, i, 1);
+      if (outsideSubId == insideSubId)
+        return i;
+    }
+    DUNE_THROW(GridError, "outsideIntersectionIndex not successfull");
+  }
+
+  bool isElementTrimmed(const IdType& elementId) const {
+    return infoFromId<0>(elementId).trimmed;
+  }
+
   template <int codim>
   requires(codim >= 0 and codim <= 2)
   const auto& entity(int lvl, int indexInLvlStorage) const {
@@ -98,9 +124,10 @@ struct VectorEntityContainer
     return entity<codim>(infoFromId<codim>(id, lvl));
   }
 
+  // todo specialize for codim 0 and 2 to not need lvl
   template <int codim>
   requires(codim >= 0 and codim <= 2)
-  const auto& infoFromId(const IdType& id, int lvl) const {
+  const auto& infoFromId(const IdType& id, int lvl = std::numeric_limits<int>::max()) const {
     if constexpr (codim == 0)
       return idToElementInfoMap.at(id);
     else if constexpr (codim == 1)
@@ -109,6 +136,7 @@ struct VectorEntityContainer
       return idToVertexInfoMap[lvl].at(id);
   }
 
+
   template <int cc>
   auto subIndexFromId(const IdType& id, int i, int codim, int lvl) const {
     assert(codim > 0);

dune/iga/trimmer/defaulttrimmer/patchgridentityseed.hh

@@ -45,7 +45,9 @@ public:
    */
   explicit PatchGridEntitySeed(const EntityImp& ent)
       : lvl_(ent.getLocalEntity().entityInfo_.lvl),
-        indexInLvlStorage_{ent.getLocalEntity().entityInfo_.indexInLvlStorage} {}
+        indexInLvlStorage_{ent.getLocalEntity().entityInfo_.indexInLvlStorage} {
+    assert(isValid());
+  }
 
   /**
    * @brief Get stored ParameterSpaceGridEntitySeed