work

doc/doxygen/modules.txt

@@ -2,6 +2,6 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later
 /*!
   @defgroup ParameterSpace ParameterSpace
- @details The group of available trimmers, see the Dune::IGA::IdentityTrim::ParameterSpace for the interface definiton
+ @details The group of available trimmers, see the Dune::IGA::IdentityParameterSpace::ParameterSpace for the interface definiton
 */
 

dune/iga/geometrykernel/algorithms.hh

@@ -2,6 +2,7 @@
 // SPDX-License-Identifier: LGPL-2.1-or-later
 #pragma once
 
+#include <dune/common/rangeutilities.hh>
 #include <dune/iga/geometrykernel/closestpointprojection.hh>
 #include <dune/iga/geometrykernel/geohelper.hh>
 
@@ -214,4 +215,48 @@ auto isPointOnLineSegment(const Coordinate& p, const Coordinate& linePoint0, con
   return false;
 }
 
+template <typename T>
+concept GeometryConcept = requires(T t, int i) {
+  T::coorddimension;
+  typename T::ctype;
+  { t.corner(i) } -> std::same_as<FieldVector<typename T::ctype, T::coorddimension>>;
+  { t.corners() } -> std::convertible_to<int>;
+};
+
+/**
+ *  @brief This returns the center of mass for a polygon
+ *
+ *  @details See https://en.wikipedia.org/wiki/Centroid#Of_a_polygon for details
+ *   TODO This actually only works if the corners are ordered consecutively along the polygons perimeter
+ */
+template <GeometryConcept Geometry>
+auto centerOfMass(const Geometry& geometry) -> FieldVector<typename Geometry::ctype, Geometry::coorddimension> {
+  int n = geometry.corners();
+
+  using Point = FieldVector<typename Geometry::ctype, Geometry::coorddimension>;
+  std::vector<Point> corners;
+  corners.reserve(n);
+
+  for (const auto i : Dune::range(n))
+    corners.push_back(geometry.corner(i));
+
+  double a  = 0.0;
+  double cx = 0.0;
+  double cy = 0.0;
+  for (const auto i : Dune::range(n - 1)) {
+    cx += (corners[i][0] + corners[i + 1][0]) * (corners[i][0] * corners[i + 1][1] - corners[i + 1][0] * corners[i][1]);
+    cy += (corners[i][1] + corners[i + 1][1]) * (corners[i][0] * corners[i + 1][1] - corners[i + 1][0] * corners[i][1]);
+    a += corners[i][0] * corners[i + 1][1] - corners[i + 1][0] * corners[i][1];
+  }
+  // Last index
+  cx += (corners[n - 1][0] + corners[0][0]) * (corners[n - 1][0] * corners[0][1] - corners[0][0] * corners[n - 1][1]);
+  cy += (corners[n - 1][1] + corners[0][1]) * (corners[n - 1][0] * corners[0][1] - corners[0][0] * corners[n - 1][1]);
+  a += corners[n - 1][0] * corners[0][1] - corners[0][0] * corners[n - 1][1];
+
+  a *= 0.5;
+  cx *= 1 / (6 * a);
+  cy *= 1 / (6 * a);
+  return {cx, cy};
+}
+
 } // namespace Dune::IGA::GeometryKernel

dune/iga/geometrykernel/nurbspatchgeometrylocalview.hh

@@ -139,11 +139,10 @@ namespace GeometryKernel {
      * @return Global coordinates of the center.
      */
     [[nodiscard]] GlobalCoordinate center() const {
-      // TODO Test this? Also not really good like that
       if constexpr (mydimension == gridDimension and not ParameterSpace::isAlwaysTrivial)
-        return global(parameterSpaceGeometry->local(parameterSpaceGeometry->center()));
+        return centerOfMass(*parameterSpaceGeometry);
       else
-        return global(LocalCoordinate(0.5));
+        return LocalCoordinate(0.5);
     }
 
     /**
@@ -208,7 +207,7 @@ namespace GeometryKernel {
      * @brief Compute the volume of the patch.
      * @return Volume of the patch.
      */
-    [[nodiscard]] double volume() const { // TODO: Test this
+    [[nodiscard]] double volume() const {
       int order = mydimension * *std::ranges::max_element(patchGeometry_->patchData_.degree);
 
       auto rule = [&]() {
@@ -221,8 +220,8 @@ namespace GeometryKernel {
 
       Volume vol = 0.0;
       for (auto& gp : rule)
-        vol +=
-            integrationElement(globalInParameterSpace(gp.position())) * gp.weight() * parameterSpaceGeometry->volume();
+        vol += integrationElement(gp.position()) * gp.weight();
+        // vol += integrationElement(globalInParameterSpace(gp.position())) * gp.weight() * parameterSpaceGeometry->volume();
       return vol;
     }
 

dune/iga/hierarchicpatch/patchgrid.hh

@@ -32,7 +32,7 @@ class NurbsPreBasis;
 }
 namespace Dune::IGA {
 
-namespace IdentityTrim {
+namespace IdentityParameterSpace {
   template <int dim, int dimworld, typename ScalarType = double>
   struct PatchGridFamily;
 }
@@ -79,7 +79,7 @@ struct HostGridAccess;
  *
  * @endcode
  */
-template <int dim, int dimworld, template <int, int, typename> typename GridFamily_ = IdentityTrim::PatchGridFamily,
+template <int dim, int dimworld, template <int, int, typename> typename GridFamily_ = IdentityParameterSpace::PatchGridFamily,
           typename ScalarType = double>
 class PatchGrid : public GridDefaultImplementation<dim, dimworld, ScalarType, GridFamily_<dim, dimworld, ScalarType>>
 {
@@ -161,15 +161,15 @@ public:
                      const std::optional<PatchTrimData>& patchTrimData            = std::nullopt,
                      const typename GridFamily::ParameterSpaceTraits::ParameterType& par = {})
       : patchGeometries_(1, GeometryKernel::NURBSPatch<dim, dimworld, ctype>(patchData)),
-        trimmer_(std::make_unique<ParameterSpace>(*this, patchTrimData, par)) {
+        parameterSpace_(std::make_unique<ParameterSpace>(*this, patchTrimData, par)) {
     patchGeometriesUnElevated = patchGeometries_;
   }
 
   PatchGrid& operator=(PatchGrid&& other) noexcept {
     patchGeometries_          = std::move(other.patchGeometries_);
     patchGeometriesUnElevated = std::move(other.patchGeometriesUnElevated);
-    trimmer_                  = std::move(other.trimmer_);
-    trimmer_->update(this);
+    parameterSpace_                  = std::move(other.parameterSpace);
+    parameterSpace_->update(this);
     return *this;
   }
 
@@ -178,7 +178,7 @@ public:
    * Levels are numbered 0 ... maxlevel with 0 the coarsest level.
    */
   [[nodiscard]] int maxLevel() const {
-    return trimmer_->maxLevel();
+    return parameterSpace_->maxLevel();
   }
 
   // Iterator to first entity of given codim on level
@@ -238,7 +238,7 @@ public:
   /** @brief returns the number of boundary segments within the macro grid
    */
   [[nodiscard]] size_t numBoundarySegments() const {
-    return trimmer_->numBoundarySegments();
+    return parameterSpace_->numBoundarySegments();
   }
 
   // number of leaf entities per codim in this process
@@ -258,31 +258,31 @@ public:
 
   /** @brief Access to the GlobalIdSet */
   const typename Traits::GlobalIdSet& globalIdSet() const {
-    return *trimmer_->globalIdSet_;
+    return *parameterSpace_->globalIdSet_;
   }
 
   /** @brief Access to the LocalIdSet */
   const typename Traits::LocalIdSet& localIdSet() const {
-    return *trimmer_->localIdSet_;
+    return *parameterSpace_->localIdSet_;
   }
 
   /** @brief Access to the LevelIndexSets */
   const typename Traits::LevelIndexSet& levelIndexSet(int level) const {
     if (level < 0 || level > maxLevel()) {
       DUNE_THROW(GridError, "levelIndexSet of nonexisting level " << level << " requested!");
     }
-    return *trimmer_->levelIndexSets_[level];
+    return *parameterSpace_->levelIndexSets_[level];
   }
 
   /** @brief Access to the LeafIndexSet */
   const typename Traits::LeafIndexSet& leafIndexSet() const {
-    return *trimmer_->leafIndexSet_;
+    return *parameterSpace_->leafIndexSet_;
   }
 
   /** @brief Create Entity from EntitySeed */
   template <class EntitySeed>
   typename Traits::template Codim<EntitySeed::codimension>::Entity entity(const EntitySeed& seed) const {
-    return trimmer_->entity(seed);
+    return parameterSpace_->entity(seed);
   }
 
   /** @name Grid Refinement Methods */
@@ -310,7 +310,7 @@ public:
       patchGeometries_.emplace_back(newfinestPatchData, newUniqueKnotVecs);
       patchGeometriesUnElevated.emplace_back(patchGeometries_.back());
 
-      trimmer_->globalRefine(1);
+      parameterSpace_->globalRefine(1);
     }
   }
 
@@ -388,52 +388,52 @@ public:
    */
   bool mark(int refCount, const typename Traits::template Codim<0>::Entity& e) {
     // @todo trim this does not do the right thing! the knotspans should also be aware of this change
-    return false; // trimmer_->parameterSpaceGrid().mark(refCount, getHostEntity<0>(e));
+    return false; // parameterSpace->parameterSpaceGrid().mark(refCount, getHostEntity<0>(e));
   }
 
   /** @brief Return refinement mark for entity
    *
    * \return refinement mark (1,0,-1)
    */
   int getMark(const typename Traits::template Codim<0>::Entity& e) const {
-    return 0; // trimmer_->parameterSpaceGrid().getMark(getHostEntity<0>(e));
+    return 0; // parameterSpace->parameterSpaceGrid().getMark(getHostEntity<0>(e));
   }
 
   /** @brief returns true, if at least one entity is marked for adaption */
   bool preAdapt() {
-    return trimmer_->parameterSpaceGrid().preAdapt();
+    return parameterSpace_->parameterSpaceGrid().preAdapt();
   }
 
   // Triggers the grid refinement process
   bool adapt() {
-    return trimmer_->parameterSpaceGrid().adapt();
+    return parameterSpace_->parameterSpaceGrid().adapt();
   }
 
   /** @brief Clean up refinement markers */
   void postAdapt() {
-    return trimmer_->parameterSpaceGrid().postAdapt();
+    return parameterSpace_->parameterSpaceGrid().postAdapt();
   }
 
   /*@}*/
 
   /** @brief Size of the overlap on the leaf level */
   unsigned int overlapSize(int codim) const {
-    return trimmer_->parameterSpaceGrid().leafGridView().overlapSize(codim);
+    return parameterSpace_->parameterSpaceGrid().leafGridView().overlapSize(codim);
   }
 
   /** @brief Size of the ghost cell layer on the leaf level */
   unsigned int ghostSize(int codim) const {
-    return trimmer_->parameterSpaceGrid().leafGridView().ghostSize(codim);
+    return parameterSpace_->parameterSpaceGrid().leafGridView().ghostSize(codim);
   }
 
   /** @brief Size of the overlap on a given level */
   unsigned int overlapSize(int level, int codim) const {
-    return trimmer_->parameterSpaceGrid().levelGridView(level).overlapSize(codim);
+    return parameterSpace_->parameterSpaceGrid().levelGridView(level).overlapSize(codim);
   }
 
   /** @brief Size of the ghost cell layer on a given level */
   unsigned int ghostSize(int level, int codim) const {
-    return trimmer_->parameterSpaceGrid().levelGridView(level).ghostSize(codim);
+    return parameterSpace_->parameterSpaceGrid().levelGridView(level).ghostSize(codim);
   }
 
 #if 0
@@ -455,13 +455,13 @@ public:
   /** @brief Communicate data of level gridView */
   template <class DataHandle>
   void communicate(DataHandle& handle, InterfaceType iftype, CommunicationDirection dir, int level) const {
-    // trimmer_->parameterSpaceGrid().levelGridView(level).communicate(handle, iftype, dir);
+    // parameterSpace->parameterSpaceGrid().levelGridView(level).communicate(handle, iftype, dir);
   }
 
   /** @brief Communicate data of leaf gridView */
   template <class DataHandle>
   void communicate(DataHandle& handle, InterfaceType iftype, CommunicationDirection dir) const {
-    // trimmer_->parameterSpaceGrid().leafGridView().communicate(handle, iftype, dir);
+    // parameterSpace->parameterSpaceGrid().leafGridView().communicate(handle, iftype, dir);
   }
 
   // **********************************************************
@@ -470,10 +470,10 @@ public:
 
   // Returns the hostgrid this PatchGrid lives in
   const ParameterSpaceGrid& parameterSpaceGrid() const {
-    return trimmer_->parameterSpaceGrid();
+    return parameterSpace_->parameterSpaceGrid();
   }
   ParameterSpaceGrid& parameterSpaceGrid() {
-    return trimmer_->parameterSpaceGrid();
+    return parameterSpace_->parameterSpaceGrid();
   }
 
   // Returns the hostgrid entity encapsulated in given PatchGrid entity
@@ -488,8 +488,8 @@ public:
     return patchGeometries_[lvl].numberOfSpans();
   }
 
-  const auto& trimmer() const {
-    return *trimmer_;
+  const auto& parameterSpace() const {
+    return *parameterSpace_;
   }
   const auto& patchGeometry(int i) const {
     return patchGeometries_.at(i);
@@ -498,11 +498,13 @@ public:
     return patchGeometries_.back();
   }
 
+  /** Obtain Integration Rule Generator (as a std::function) */
   auto integrationRule() const
   requires(not ParameterSpace::isAlwaysTrivial and dim == 2)
   {
     return integrationRuleHolder_.integrationRule();
   }
+  /** Dynamically set integration Rule (as a std::function, see IntegrationRuleHolder for more details) */
   void integrationRule(typename IntegrationRuleHolder<PatchGrid>::FunctionType integrationRule)
   requires(not ParameterSpace::isAlwaysTrivial and dim == 2)
   {
@@ -514,7 +516,7 @@ private:
   std::vector<GeometryKernel::NURBSPatch<dim, dimworld, ctype>> patchGeometries_;
   std::vector<GeometryKernel::NURBSPatch<dim, dimworld, ctype>> patchGeometriesUnElevated;
 
-  std::unique_ptr<ParameterSpace> trimmer_;
+  std::unique_ptr<ParameterSpace> parameterSpace_;
 
   IntegrationRuleHolder<PatchGrid> integrationRuleHolder_{};
 

dune/iga/hierarchicpatch/patchgridentity.hh

@@ -108,7 +108,7 @@ public:
 
   /** @brief returns  EntitySeed */
   EntitySeed seed() const {
-    return patchGrid_->trimmer_->seed(*this);
+    return patchGrid_->parameterSpace().seed(*this);
   }
 
   /** @brief  level of this element */
@@ -240,7 +240,7 @@ public:
 
   /** @brief Return EntitySeed */
   [[nodiscard]] EntitySeed seed() const {
-    return patchGrid_->trimmer_->seed(*this);
+    return patchGrid_->parameterSpace().seed(*this);
   }
 
   /** @brief  Level of this element */
@@ -278,22 +278,22 @@ public:
 
   // First level intersection
   [[nodiscard]] LevelIntersectionIterator ilevelbegin() const {
-    return patchGrid_->trimmer_->ilevelbegin(*this);
+    return patchGrid_->parameterSpace().ilevelbegin(*this);
   }
 
   // Reference to one past the last neighbor
   LevelIntersectionIterator ilevelend() const {
-    return patchGrid_->trimmer_->ilevelend(*this);
+    return patchGrid_->parameterSpace().ilevelend(*this);
   }
 
   // First leaf intersection
   LeafIntersectionIterator ileafbegin() const {
-    return patchGrid_->trimmer_->ileafbegin(*this);
+    return patchGrid_->parameterSpace().ileafbegin(*this);
   }
 
   // Reference to one past the last leaf intersection
   LeafIntersectionIterator ileafend() const {
-    return patchGrid_->trimmer_->ileafend(*this);
+    return patchGrid_->parameterSpace().ileafend(*this);
   }
 
   // returns true if Entity has NO children

dune/iga/hierarchicpatch/patchgridfactory.hh

@@ -86,7 +86,7 @@ public:
 
 template <int dim, int dimworld, typename ScalarType = double>
 auto makePatchGridFactory() {
-  using PatchGrid = IGA::PatchGrid<dim, dimworld, IGA::IdentityTrim::PatchGridFamily, ScalarType>;
+  using PatchGrid = IGA::PatchGrid<dim, dimworld, IGA::IdentityParameterSpace::PatchGridFamily, ScalarType>;
   return GridFactory<PatchGrid>{};
 }
 
@@ -102,7 +102,7 @@ constexpr auto withTrimmingCapabilities() {
 
 template <int dim, int dimworld, typename ScalarType = double>
 auto makePatchGridFactory(Impl::UseTrimmingCapabilities) {
-  using PatchGrid = IGA::PatchGrid<dim, dimworld, IGA::DefaultTrim::PatchGridFamily, ScalarType>;
+  using PatchGrid = IGA::PatchGrid<dim, dimworld, IGA::DefaultParameterSpace::PatchGridFamily, ScalarType>;
   return GridFactory<PatchGrid>{};
 }
 

dune/iga/hierarchicpatch/patchgridgeometry.hh

@@ -75,6 +75,10 @@ public:
     return geometryLocalView_.center();
   }
 
+  [[nodiscard]] Volume volume() const {
+    return geometryLocalView_.volume();
+  }
+
   /** @brief Maps a local coordinate within reference element to
    * global coordinate in element  */
   [[nodiscard]] GlobalCoordinate global(const LocalCoordinate& local) const {