Merge branch 'feature/curvedgeometry3' into 'master'

Curved Geometry from LocalFunction

See merge request iwr/dune-curvedsurfacegrid!5

dune/curvedsurfacegrid/gridfunctions/gridfunction.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_CURVED_SURFACE_GRID_GRIDFUNCTION_HH
+#define DUNE_CURVED_SURFACE_GRID_GRIDFUNCTION_HH
+
+#include <type_traits>
+#include <dune/curvedsurfacegrid/concepts.hh>
+#include <dune/curvedsurfacegrid/gridfunctions/analyticgridfunction.hh>
+
+namespace Dune
+{
+  namespace Impl
+  {
+    template< class ES, class = void >
+    struct GridOf;
+
+    template< class ES >
+    struct GridOf<ES, std::void_t<typename ES::GridView>>
+    {
+      using type = typename ES::GridView::Grid;
+    };
+
+    template< class ES >
+    struct GridOf<ES, std::void_t<typename ES::Grid>>
+    {
+      using type = typename ES::Grid;
+    };
+
+  } // end namespace Impl
+
+  //! Type-Traits to extract the grid-type from a GridFunction `GF`
+  template< class GF >
+  using GridOf_t = typename Impl::GridOf<std::decay_t<decltype(std::declval<GF>().entitySet())>>::type;
+
+  //! Conditionally define the GridFunction type.
+  /**
+   * If the type `GF` is a GridFunction on the grid `Grid`, use this type as GridFunction type,
+   * Otherwise it is a functor and construct an \ref AnalyticGridFunction on the `Grid`.
+   **/
+  template< class Grid, class GF >
+  using GridFunctionOf_t
+    = std::conditional_t<Dune::Concept::isGridFunction<GF, Grid>(), GF, AnalyticGridFunction<Grid,GF>>;
+
+} // end namespace Dune
+
+#endif // DUNE_CURVED_SURFACE_GRID_GRIDFUNCTION_HH

dune/curvedsurfacegrid/gridfunctions/normalgridviewfunction.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_CURVED_SURFACE_GRID_NORMAL_GRIDVIEWFUNCTION_HH
+#define DUNE_CURVED_SURFACE_GRID_NORMAL_GRIDVIEWFUNCTION_HH
+
+#include <array>
+#include <vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+#include <dune/functions/functionspacebases/basistags.hh>
+#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+#include <dune/grid/utility/hierarchicsearch.hh>
+#include <dune/istl/bvector.hh>
+
+
+namespace Dune 
+{
+  //! Grid-view function representing averaged normal vector
+  /**
+   * \tparam GridView   The grid-view this grid-view-function is defined on
+   * \tparam ORDER      Polynomial order of the lagrange bases used for representing the normals
+   * \tparam T          Value type used for the basis and the coefficients
+   **/
+  template< class GridView, int ORDER = -1, class T = double >
+  class NormalGridViewFunction
+  {
+    static auto makeBasis (const GridView& gridView, int order)
+    {
+      namespace BF = BasisFactory;
+      return BF::makeBasis(gridView, BF::power<GridView::dimensionworld>(BF::lagrange<T>(order), BF::blockedInterleaved()));
+    }
+
+    using Basis = decltype(makeBasis(std::declval<GridView>(), ORDER));
+
+  public:
+    using EntitySet = GridViewEntitySet<GridView,0>;
+
+    using Domain = typename EntitySet::GlobalCoordinate;
+    using Range = FieldVector<T,GridView::dimensionworld>;
+
+    using VectorType = BlockVector<Range>;
+
+  private:
+    class LocalFunction
+    {
+      using LocalView = typename Basis::LocalView;
+      using LocalContext = typename LocalView::Element;
+
+      using Domain = typename EntitySet::LocalCoordinate;
+      using Range = typename NormalGridViewFunction::Range;
+
+    public:
+      LocalFunction (LocalView&& localView, const VectorType& normals)
+        : localView_(std::move(localView))
+        , normals_(normals)
+      {}
+
+      //! Collect the normal vector from all element DOFs into a local
+      //! vector that can be accessed in the operator() for interpolation
+      void bind (const LocalContext& element)
+      {
+        localView_.bind(element);
+
+        const auto& leafNode = localView_.tree().child(0);
+        localNormals_.resize(leafNode.size());
+
+        // collect local normal vectors
+        for (std::size_t i = 0; i < localNormals_.size(); ++i) {
+          auto idx = localView_.index(leafNode.localIndex(i));
+          localNormals_[i] = normals_[idx[0]];
+        }
+
+        bound_ = true;
+      }
+
+      void unbind ()
+      {
+        localView_.unbind();
+        bound_ = false;
+      }
+
+      // evaluate normal vectors in local coordinate
+      // by interpolation of stored local normals.
+      Range operator() (const Domain& local) const 
+      {
+        assert(bound_);
+
+        const auto& leafNode = localView_.tree().child(0);
+        const auto& lfe = leafNode.finiteElement();
+
+        // evaluate basis functions in local coordinate
+        lfe.localBasis().evaluateFunction(local, shapeValues_);
+        assert(localNormals_.size() == shapeValues_.size());
+
+        Range n(0);
+        for (std::size_t i = 0; i < localNormals_.size(); ++i)
+          n.axpy(shapeValues_[i], localNormals_[i]);
+
+        // return normalized vector
+        return n / n.two_norm();
+      }
+
+    private:
+      LocalView localView_;
+      const VectorType& normals_;
+
+      std::vector<Range> localNormals_;
+      mutable std::vector<FieldVector<T,1>> shapeValues_;
+      bool bound_ = false;
+    };
+
+  public:
+    //! Constructor of the grid function. 
+    /** 
+     * Creates a global basis of a power of langrange nodes of given order.
+     * The constructor argument `order` is defaulted to the class template parameter.
+     **/
+    NormalGridViewFunction (const GridView& gridView, int order = ORDER)
+      : entitySet_(gridView)
+      , basis_(makeBasis(gridView, order))
+    {
+      update(gridView);
+    }
+
+    //! Epdate the grid function. 
+    /**
+     * This calculates a mean average of normal vectors in the  DOFs of the basis. 
+     * Those averages are stored normalized in the coefficients vector.
+     **/
+    void update (const GridView& gridView)
+    {
+      entitySet_ = EntitySet{gridView};
+      basis_.update(gridView);
+
+      normals_.resize(basis_.size());
+      normals_ = 0;
+
+      // compute normal vectors by mean averaging
+      auto localView = basis_.localView();
+      for (const auto& e : elements(basis_.gridView()))
+      {
+        localView.bind(e);
+        auto geometry = e.geometry();
+
+        const auto& leafNode = localView.tree().child(0);
+        const auto& lfe = leafNode.finiteElement();
+
+        // interpolate normal of geometry
+        std::vector<Range> localNormals;
+        lfe.localInterpolation().interpolate([&](const auto& local) -> Range {
+          return Dune::normal(geometry, local);
+        }, localNormals);
+
+        // copy to global vector
+        for (std::size_t i = 0; i < localNormals.size(); ++i) {
+          auto idx = localView.index(leafNode.localIndex(i));
+          normals_[idx[0]] += localNormals[i];
+        }
+      }
+
+      // normalize vector
+      for (std::size_t i = 0; i < normals_.size(); ++i)
+        normals_[i] /= normals_[i].two_norm();
+    }
+
+    //! Evaluate normal vectors in global coordinates
+    // NOTE: expensive
+    Range operator() (const Domain& x) const 
+    {
+      using Grid = typename GridView::Grid;
+      using IS = typename GridView::IndexSet;
+
+      const auto& gv = basis_.gridView();
+      HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
+
+      auto element = hsearch.findEntity(x);
+      auto geometry = element.geometry();
+      auto localFct = localFunction(*this);
+      localFct.bind(element);
+      return localFct(geometry.local(x));
+    }
+
+    //! Create a local function of this grifunction
+    friend LocalFunction localFunction (const NormalGridViewFunction& gf)
+    {
+      return LocalFunction{gf.basis_.localView(), gf.normals_};
+    }
+
+    //! obtain the stored \ref GridViewEntitySet
+    const EntitySet& entitySet () const
+    {
+      return entitySet_;
+    }
+
+  private:
+    EntitySet entitySet_;
+    Basis basis_;
+    VectorType normals_;
+  };
+
+} // end namespace Dune 
+
+#endif // DUNE_CURVED_SURFACE_GRID_NORMAL_GRIDVIEWFUNCTION_HH

dune/curvedsurfacegrid/gridfunctions/spheregridfunction.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH
+#define DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH
+
+#include <type_traits>
+
+#include <dune/curvedsurfacegrid/gridfunctions/analyticgridfunction.hh>
+#include <dune/functions/common/defaultderivativetraits.hh>
+
+namespace Dune
+{
+  // Sphere functor
+  template< class T >
+  struct SphereProjection
+  {
+    T radius_;
+
+    SphereProjection (T radius)
+      : radius_(radius)
+    {}
+
+    //! project the coordinate to the sphere at origin with `radius`
+    template< class Domain >
+    Domain operator() (const Domain& x) const
+    {
+      return x * (radius_ / x.two_norm());
+    }
+
+    //! derivative of the projection
+    friend auto derivative (const SphereProjection& sphere)
+    {
+      return [r=sphere.radius_](auto const& x)
+      {
+        using Domain = std::decay_t<decltype(x)>;
+        using DerivativeTraits = Functions::DefaultDerivativeTraits<Domain(Domain)>;
+        typename DerivativeTraits::Range out;
+
+        auto nrm = x.two_norm();
+        for (int i = 0; i < out.N(); ++i)
+          for (int j = 0; j < out.M(); ++j)
+            out[i][j] = r * ((i == j ? 1 : 0) - (x[i]/nrm) * (x[j]/nrm)) / nrm;
+        return out;
+      };
+    }
+
+    //! normal vector = x/|x|
+    template< class Domain >
+    Domain normal (const Domain& x) const
+    {
+      return x / x.two_norm();
+    }
+
+    //! mean curvature of the sphere = 1/R
+    template< class Domain >
+    T mean_curvature (const Domain& /*x*/) const
+    {
+      return T(1)/radius_;
+    }
+
+    //! surface area of the sphere = 4*pi*r^2
+    T area () const
+    {
+      return 4*M_PI*radius_*radius_;
+    }
+  };
+
+  //! construct a grid function representing a sphere parametrization
+  template< class Grid, class T >
+  auto sphereGridFunction (T radius)
+  {
+    return analyticGridFunction<Grid>(SphereProjection<T>{radius});
+  }
+
+} // end namespace Dune
+
+#endif // DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH

dune/curvedsurfacegrid/surfacedistance/CMakeLists.txt

@@ -7,7 +7,9 @@ dune_add_library("surfacedistance" OBJECT
   vtureader.hh)
 
 #install headers
-install(FILES datastructures.hh
-              surfacedistance.hh
-              vtureader.hh
-        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/curvedsurfacegrid/surfacedistance)
+install(FILES 
+  datastructures.hh
+  surfacedistance.hh
+  vertexmap.hh
+  vtureader.hh
+  DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/curvedsurfacegrid/surfacedistance)