Merge branch 'feature/curved-geometries' into 'master'

Add geometry based on a local-finite-element parametrization

See merge request !169

CHANGELOG.md

@@ -21,6 +21,9 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
   default construction results in an "empty"/invalid geometry that can be assigned
   a valid geometry.
 
+- Add a geometry `LocalFiniteElementGeometry` parametrized by local finite-element
+  basis functions.
+
 ## Deprecations and removals
 
 - `Dune::Transitional::ReferenceElement` is deprecated and will be removed after

dune/geometry/CMakeLists.txt

@@ -13,6 +13,7 @@ install(FILES
   generalvertexorder.hh
   mappedgeometry.hh
   multilineargeometry.hh
+  localfiniteelementgeometry.hh
   quadraturerules.hh
   referenceelement.hh
   referenceelementimplementation.hh

dune/geometry/localfiniteelementgeometry.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
+// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
+#ifndef DUNE_GEOMETRY_PARAMETRIZEDGEOMETRY_HH
+#define DUNE_GEOMETRY_PARAMETRIZEDGEOMETRY_HH
+
+#include <cassert>
+#include <functional>
+#include <limits>
+#include <type_traits>
+#include <vector>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/math.hh>
+#include <dune/common/typetraits.hh>
+#include <dune/common/std/type_traits.hh>
+
+#include <dune/geometry/affinegeometry.hh> // for FieldMatrixHelper
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/referenceelements.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/utility/algorithms.hh>
+#include <dune/geometry/utility/convergence.hh>
+
+namespace Dune {
+
+/**
+ * \brief Geometry implementation based on local-basis function parametrization.
+ *
+ * Parametrization of the geometry by any localfunction interpolated into a local
+ * finite-element space.
+ *
+ * \tparam  LFE   Type of a local finite-element.
+ * \tparam  cdim  Coordinate dimension.
+ */
+template <class LFE, int cdim>
+class LocalFiniteElementGeometry
+{
+  using LocalFiniteElement = LFE;
+  using LocalBasis = typename LFE::Traits::LocalBasisType;
+  using LocalBasisTraits = typename LocalBasis::Traits;
+
+public:
+  /// coordinate type
+  using ctype = typename LocalBasisTraits::DomainFieldType;
+
+  /// geometry dimension
+  static const int mydimension = LocalBasisTraits::dimDomain;
+
+  /// coordinate dimension
+  static const int coorddimension = cdim;
+
+  /// type of local coordinates
+  using LocalCoordinate = FieldVector<ctype, mydimension>;
+
+  /// type of global coordinates
+  using GlobalCoordinate = FieldVector<ctype, coorddimension>;
+
+  /// type of volume
+  using Volume = decltype(power(std::declval<ctype>(),mydimension));
+
+  /// type of jacobian
+  using Jacobian = FieldMatrix<ctype, coorddimension, mydimension>;
+
+  /// type of jacobian transposed
+  using JacobianTransposed = FieldMatrix<ctype, mydimension, coorddimension>;
+
+  /// type of jacobian inverse
+  using JacobianInverse = FieldMatrix<ctype, mydimension, coorddimension>;
+
+  /// type of jacobian inverse transposed
+  using JacobianInverseTransposed = FieldMatrix<ctype, coorddimension, mydimension>;
+
+public:
+  /// type of reference element
+  using ReferenceElements = Dune::ReferenceElements<ctype, mydimension>;
+  using ReferenceElement = typename ReferenceElements::ReferenceElement;
+
+protected:
+  using MatrixHelper = Impl::FieldMatrixHelper<ctype>;
+
+public:
+  /// \brief Default constructed geometry results in an empty/invalid representation.
+  LocalFiniteElementGeometry () = default;
+
+  /**
+   * \brief Constructor from a vector of coefficients of the LocalBasis parametrizing
+   *        the geometry.
+   *
+   * \param[in]  refElement  reference element for the geometry
+   * \param[in]  localFE     Local finite-element to use for the parametrization
+   * \param[in]  vertices    Coefficients of the local interpolation into the basis
+   *
+   * The vertices are the coefficients of a local interpolation in the local
+   * finite-element. For Lagrange local bases, these correspond to vertices on
+   * the curved geometry in the local Lagrange nodes.
+   *
+   * \note The vertices are stored internally, so if possible move an external
+   *       vertex storage to this constructor
+   **/
+  LocalFiniteElementGeometry (const ReferenceElement& refElement,
+                              const LocalFiniteElement& localFE,
+                              std::vector<GlobalCoordinate> vertices)
+    : refElement_(refElement)
+    , localFE_(localFE)
+    , vertices_(std::move(vertices))
+  {
+    assert(localFE_.size() == vertices_.size());
+  }
+
+  /**
+   * \brief Constructor from a local parametrization function, mapping local to
+   *        (curved) global coordinates.
+   *
+   * \param[in]  refElement      reference element for the geometry
+   * \param[in]  localFE         Local finite-element to use for the parametrization
+   * \param[in]  parametrization parametrization function with signature
+   *                             `GlobalCoordinate(LocalCoordinate)`
+   *
+   * The parametrization function is not stored in the class, but interpolated into
+   * the local finite-element basis and the computed interpolation coefficients
+   * are stored.
+   **/
+  template <class Param,
+    std::enable_if_t<std::is_invocable_r_v<GlobalCoordinate,Param,LocalCoordinate>, int> = 0>
+  LocalFiniteElementGeometry (const ReferenceElement& refElement,
+                              const LocalFiniteElement& localFE,
+                              Param&& parametrization)
+    : refElement_(refElement)
+    , localFE_(localFE)
+  {
+    localFE_.localInterpolation().interpolate(parametrization, vertices_);
+  }
+
+  /**
+   * \brief Constructor, forwarding to the other constructors that take a reference-element.
+   *
+   * \param[in]  gt       geometry type
+   * \param[in]  args...  arguments passed to the other constructors
+   **/
+  template <class... Args>
+  explicit LocalFiniteElementGeometry (GeometryType gt, Args&&... args)
+    : LocalFiniteElementGeometry(ReferenceElements::general(gt), std::forward<Args>(args)...)
+  {}
+
+  /// \brief Obtain the polynomial order of the parametrization.
+  int order () const
+  {
+    return localBasis().order();
+  }
+
+  /// \brief Is this mapping affine? This is only true for flat affine geometries.
+  bool affine () const
+  {
+    if (!affine_)
+      affine_.emplace(affineImpl());
+    return *affine_;
+  }
+
+  /// \brief Obtain the name of the reference element.
+  GeometryType type () const
+  {
+    return refElement_.type();
+  }
+
+  /// \brief Obtain number of corners of the corresponding reference element.
+  int corners () const
+  {
+    return refElement_.size(mydimension);
+  }
+
+  /// \brief Obtain coordinates of the i-th corner.
+  GlobalCoordinate corner (int i) const
+  {
+    assert( (i >= 0) && (i < corners()) );
+    return global(refElement_.position(i, mydimension));
+  }
+
+  /// \brief Obtain the centroid of the mapping's image.
+  GlobalCoordinate center () const
+  {
+    return global(refElement_.position(0, 0));
+  }
+
+  /**
+   * \brief Evaluate the coordinate mapping.
+   *
+   * Implements a linear combination of local basis functions scaled by
+   * the vertices as coefficients.
+   *
+   * \f[ global = \sum_i v_i \psi_i(local) \f]
+   *
+   * \param[in] local  local coordinate to map
+   * \returns          corresponding global coordinate
+   **/
+  GlobalCoordinate global (const LocalCoordinate& local) const
+  {
+    thread_local std::vector<typename LocalBasisTraits::RangeType> shapeValues;
+    localBasis().evaluateFunction(local, shapeValues);
+    assert(shapeValues.size() == vertices_.size());
+
+    GlobalCoordinate out(0);
+    for (std::size_t i = 0; i < shapeValues.size(); ++i)
+      out.axpy(shapeValues[i], vertices_[i]);
+
+    return out;
+  }
+
+  /**
+   * \brief Evaluate the inverse coordinate mapping.
+   *
+   * \param[in] y  Global coordinate to map
+   * \param[in] opts  Parameters to control the behavior of the Gauss-Newton
+   *                  algorithm.
+   *
+   * \return For given global coordinate `y` the local coordinate `x` that minimizes
+   *         the function `(global(x) - y).two_norm2()` over the local coordinate
+   *         space spanned by the reference element.
+   *
+   *  \throws In case of an error indicating that the local coordinate can not be
+   *          obtained, an exception is thrown, with an error code from
+   *          \ref `GaussNewtonErrorCode`.
+   *  \note It is not guaranteed that the resulting local coordinate is inside the
+   *        reference element domain.
+   **/
+  LocalCoordinate local (const GlobalCoordinate& y, Impl::GaussNewtonOptions<ctype> opts = {}) const
+  {
+    LocalCoordinate x = refElement_.position(0,0);
+    Impl::GaussNewtonErrorCode err = Impl::gaussNewton(
+      [&](const LocalCoordinate& local) { return this->global(local); },
+      [&](const LocalCoordinate& local) { return this->jacobianTransposed(local); },
+      y, x, opts
+    );
+
+    if (err != Impl::GaussNewtonErrorCode::OK)
+      DUNE_THROW(Dune::Exception,
+        "Local coordinate can not be recovered from global coordinate, error code = " << int(err) << "\n"
+        << "  (global(x) - y).two_norm() = " << (global(x) - y).two_norm()
+        << " > tol = " << opts.absTol);
+
+    return x;
+  }
+
+  /**
+   * \brief Obtain the integration element.
+   *
+   * If the Jacobian of the geometry is denoted by $J(x)$, the integration element
+   * \f$\mu(x)\f$ is given by \f[ \mu(x) = \sqrt{|\det (J^T(x) J(x))|}.\f]
+   *
+   * \param[in]  local  local coordinate to evaluate the integration element in.
+   *
+   * \returns the integration element \f$\mu(x)\f$.
+   **/
+  ctype integrationElement (const LocalCoordinate& local) const
+  {
+    return MatrixHelper::sqrtDetAAT(jacobianTransposed(local));
+  }
+
+  /**
+   * \brief Obtain the volume of the mapping's image.
+   *
+   * Calculates the volume of the entity by numerical integration. Since the
+   * polynomial order of the volume element is not known, iteratively compute
+   * numerical integrals with increasing order of the quadrature rules, until
+   * tolerance is reached.
+   **/
+  Volume volume (Impl::ConvergenceOptions<ctype> opts = {}) const
+  {
+    Volume vol0 = volume(QuadratureRules<ctype, mydimension>::rule(type(), 1));
+    if (affine())
+      return vol0;
+
+    using std::abs;
+    for (int p = 2; p < opts.maxIt; ++p) {
+      Volume vol1 = volume(QuadratureRules<ctype, mydimension>::rule(type(), p));
+      if (abs(vol1 - vol0) < opts.absTol)
+        return vol1;
+
+      vol0 = vol1;
+    }
+    return vol0;
+  }
+
+  /// \brief Obtain the volume of the mapping's image by given quadrature rules.
+  Volume volume (const QuadratureRule<ctype, mydimension>& quadRule) const
+  {
+    Volume vol(0);
+    for (const auto& qp : quadRule)
+      vol += integrationElement(qp.position()) * qp.weight();
+    return vol;
+  }
+
+  /**
+   * \brief Obtain the Jacobian.
+   * \param[in]  local  local coordinate to evaluate Jacobian in
+   * \returns the matrix corresponding to the Jacobian
+   **/
+  Jacobian jacobian (const LocalCoordinate& local) const
+  {
+    thread_local std::vector<typename LocalBasisTraits::JacobianType> shapeJacobians;
+    localBasis().evaluateJacobian(local, shapeJacobians);
+    assert(shapeJacobians.size() == vertices_.size());
+
+    Jacobian out(0);
+    for (std::size_t i = 0; i < shapeJacobians.size(); ++i) {
+      for (int j = 0; j < Jacobian::rows; ++j) {
+        shapeJacobians[i].umtv(vertices_[i][j], out[j]);
+      }
+    }
+    return out;
+  }
+
+  /**
+   * \brief Obtain the transposed of the Jacobian.
+   * \param[in]  local  local coordinate to evaluate Jacobian in
+   * \returns the matrix corresponding to the transposed of the Jacobian
+   **/
+  JacobianTransposed jacobianTransposed (const LocalCoordinate& local) const
+  {
+    return jacobian(local).transposed();
+  }
+
+  /**
+   * \brief Obtain the Jacobian's inverse.
+   *
+   * The Jacobian's inverse is defined as a pseudo-inverse. If we denote the
+   * Jacobian by \f$J(x)\f$, the following condition holds:
+   *   \f[ J^{-1}(x) J(x) = I. \f]
+   **/
+  JacobianInverse jacobianInverse (const LocalCoordinate& local) const
+  {
+    JacobianInverse out;
+    MatrixHelper::leftInvA(jacobian(local), out);
+    return out;
+  }
+
+  /**
+   * \brief Obtain the transposed of the Jacobian's inverse.
+   *
+   * The Jacobian's inverse is defined as a pseudo-inverse. If we denote the
+   * Jacobian by \f$J(x)\f$, the following condition holds:
+   *   \f[ J^{-1}(x) J(x) = I. \f]
+   **/
+  JacobianInverseTransposed jacobianInverseTransposed (const LocalCoordinate& local) const
+  {
+    return jacobianInverse(local).transposed();
+  }
+
+  /// \brief Obtain the reference-element related to this geometry.
+  friend ReferenceElement referenceElement (const LocalFiniteElementGeometry& geometry)
+  {
+    return geometry.refElement_;
+  }
+
+  /// \brief Obtain the local finite-element.
+  const LocalFiniteElement& finiteElement () const
+  {
+    return localFE_;
+  }
+
+  /// \brief Obtain the coefficients of the parametrization.
+  const std::vector<GlobalCoordinate>& coefficients () const
+  {
+    return vertices_;
+  }
+
+  /// \brief The local basis of the stored local finite-element.
+  const LocalBasis& localBasis () const
+  {
+    return localFE_.localBasis();
+  }
+
+private:
+
+  bool affineImpl () const
+  {
+    if constexpr(mydimension == 0)
+      // point geometries are always affine mappings
+      return true;
+    else {
+      if (order() > 1)
+        // higher-order parametrizations are by definition not affine
+        return false;
+      if constexpr(mydimension == 1)
+        // linear line geometries are affine mappings
+        return true;
+      else {
+        if (type().isSimplex())
+          // linear simplex geometries are affine mappings
+          return true;
+
+        // multi-linear mappings on non-simplex geometries might be affine
+        // as well. This is tested explicitly for all vertices by constructing
+        // an affine mapping from dim+1 affine-independent corners and evaluating
+        // at the other corners.
+        auto refSimplex = referenceElement<ctype,mydimension>(GeometryTypes::simplex(mydimension));
+        auto simplexIndices = Dune::range(refSimplex.size(mydimension));
+        auto simplexCorners = Dune::transformedRangeView(simplexIndices,
+          [&](int i) { return this->global(refSimplex.position(i,mydimension)); });
+        AffineGeometry<ctype,mydimension,coorddimension> affineGeo(refSimplex,simplexCorners);
+        using std::sqrt;
+        const ctype tol = sqrt(std::numeric_limits<ctype>::epsilon());
+        for (int i = 0; i < corners(); ++i) {
+          const auto corner = refElement_.position(i,mydimension);
+          if ((affineGeo.global(corner) - global(corner)).two_norm() > tol)
+            return false;
+        }
+        return true;
+      }
+    }
+  }
+
+private:
+  /// Reference of the geometry
+  ReferenceElement refElement_{};
+
+  /// A local finite-element
+  LocalFiniteElement localFE_{};
+
+  /// The (Lagrange) coefficients of the interpolating geometry
+  std::vector<GlobalCoordinate> vertices_{};
+
+  mutable std::optional<bool> affine_ = std::nullopt;
+};
+
+namespace Impl {
+
+// extract the LocalCoordinate type from a LocalFiniteElement
+template <class LFE>
+using LocalCoordinate_t
+  = FieldVector<typename LFE::Traits::LocalBasisType::Traits::DomainFieldType,
+                LFE::Traits::LocalBasisType::Traits::dimDomain>;
+
+} // end namespace Impl
+
+
+// deduction guides
+template <class I, class LFE, class GlobalCoordinate>
+LocalFiniteElementGeometry (Geo::ReferenceElement<I>, const LFE&, std::vector<GlobalCoordinate>)
+  -> LocalFiniteElementGeometry<LFE, GlobalCoordinate::dimension>;
+
+template <class I, class LFE, class F,
+          class Range = std::invoke_result_t<F,Impl::LocalCoordinate_t<LFE>>>
+LocalFiniteElementGeometry (Geo::ReferenceElement<I>, const LFE&, const F&)
+  -> LocalFiniteElementGeometry<LFE, Range::dimension>;
+
+template <class LFE, class GlobalCoordinate>
+LocalFiniteElementGeometry (GeometryType, const LFE& localFE, std::vector<GlobalCoordinate>)
+  -> LocalFiniteElementGeometry<LFE, GlobalCoordinate::dimension>;
+
+template <class LFE, class F,
+          class Range = std::invoke_result_t<F,Impl::LocalCoordinate_t<LFE>>>
+LocalFiniteElementGeometry (GeometryType, const LFE&, const F&)
+  -> LocalFiniteElementGeometry<LFE, Range::dimension>;
+
+} // namespace Dune
+
+#endif // DUNE_GEOMETRY_PARAMETRIZEDGEOMETRY_HH

dune/geometry/test/CMakeLists.txt

 # SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
 # SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 
+dune_add_test(SOURCES benchmark-geometries.cc
+              LINK_LIBRARIES dunegeometry)
+
 dune_add_test(SOURCES test-affinegeometry.cc
               LINK_LIBRARIES dunegeometry)
 
@@ -27,6 +30,9 @@ dune_add_test(SOURCES test-multilineargeometry.cc
 dune_add_test(SOURCES test-nonetype.cc
               LINK_LIBRARIES dunegeometry)
 
+dune_add_test(SOURCES test-localfiniteelementgeometry.cc
+              LINK_LIBRARIES dunegeometry)
+
 dune_add_test(SOURCES test-refinement.cc
               LINK_LIBRARIES dunegeometry)
 

dune/geometry/test/benchmark-geometries.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
+// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
+#include <config.h>
+
+#include <limits>
+#include <cmath>
+#include <type_traits>
+#include <vector>
+
+#include <dune/common/timer.hh>
+#include <dune/geometry/localfiniteelementgeometry.hh>
+#include <dune/geometry/mappedgeometry.hh>
+#include <dune/geometry/multilineargeometry.hh>
+#include <dune/geometry/referenceelements.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/test/checkgeometry.hh>
+#include <dune/geometry/test/comparegeometries.hh>
+#include <dune/geometry/test/localfiniteelements.hh>
+#include <dune/geometry/test/referenceelementgeometry.hh>
+
+
+template <class Geometry>
+bool benchmarkGeometry (const Geometry& geo)
+{
+  bool pass = true;
+  bool isAffine = geo.affine();
+  Dune::GeometryType t = geo.type();
+
+  const auto& quadrature = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geo.type(), 8);
+  for (auto&& [pos,weight] : quadrature)
+  {
+    pass &= (geo.affine() == isAffine);
+    pass &= (geo.type() == t);
+    pass &= (geo.corners() > 0);
+
+    for (int i = 0; i < geo.corners(); ++i) {
+      pass &= ((geo.corner(i) - geo.corner((i+1)%geo.corners())).two_norm() > 0);
+      pass &= ((geo.corner(i) - geo.center()).two_norm() > 0);
+    }
+
+    pass &= (geo.volume() > 0);
+    pass &= (geo.global(pos).size() == Geometry::coorddimension);
+    pass &= (geo.jacobian(pos).M() == Geometry::mydimension);
+    pass &= (geo.jacobianTransposed(pos).N() == Geometry::mydimension);
+    pass &= (geo.jacobianInverse(pos).M() == Geometry::coorddimension);
+    pass &= (geo.jacobianInverseTransposed(pos).N() == Geometry::coorddimension);
+    pass &= (geo.integrationElement(pos) > 0);
+  }
+
+  return pass;
+}
+
+
+template <class ctype, int cdim, Dune::GeometryType::Id id>
+bool benchmarkGeometries (int nIter = 100)
+{
+  bool pass = true;
+
+  constexpr auto gt = Dune::GeometryType{id};
+  auto refElem = Dune::referenceElement<ctype,gt.dim()>(gt);
+  std::cout << "Time(" << refElem.type() << "):" << std::endl;
+
+  using LFE = std::conditional_t<
+    gt.isSimplex(), Dune::Impl::P1LocalFiniteElement<ctype,ctype,gt.dim()>, std::conditional_t<
+    gt.isCube(),    Dune::Impl::Q1LocalFiniteElement<ctype,ctype,gt.dim()>, void>
+    >;
+  auto lfe = LFE{};
+
+  Dune::FieldMatrix<ctype,cdim,gt.dim()> A;
+  for (int i = 0; i < cdim; ++i)
+    for (int j = 0; j < int(gt.dim()); ++j)
+      A[i][j] = ctype(3-std::abs(i-j));
+
+  Dune::FieldVector<ctype,cdim> b;
+  for (std::size_t i = 0; i < cdim; ++i)
+    b[i] = ctype(i+1);
+
+  // mapping to generate coordinates from reference-element corners
+  auto f = [&](Dune::FieldVector<ctype,gt.dim()> const& x) {
+    Dune::FieldVector<ctype,cdim> y;
+    A.mv(x,y);
+    y += b;
+    return y;
+  };
+
+  auto corners = std::vector<Dune::FieldVector<ctype,cdim>>(refElem.size(gt.dim()));
+  for (int i = 0; i < refElem.size(gt.dim()); ++i)
+    corners[i] = f(refElem.position(i, gt.dim()));
+
+  Dune::Timer t;
+
+  // construct a geometry using given parametrization coefficients
+  using Geometry = Dune::LocalFiniteElementGeometry<LFE,cdim>;
+  auto geometry = Geometry{refElem, lfe, corners};
+  t.reset();
+  for (int i = 0; i < nIter; ++i)
+    pass &= benchmarkGeometry(geometry);
+  std::cout << "  LocalFiniteElementGeometry = " << t.elapsed() << "sec" << std::endl;
+
+  // compare against MultiLinearGeometry
+  using MLGeometry = Dune::MultiLinearGeometry<ctype,gt.dim(),cdim>;
+  auto mlgeometry = MLGeometry{refElem, corners};
+  t.reset();
+  for (int i = 0; i < nIter; ++i)
+    pass &= benchmarkGeometry(mlgeometry);
+  std::cout << "  MultiLinearGeometry = " << t.elapsed() << "sec" << std::endl;
+
+  // compare against a MappedGeometry
+  using Mapping = Dune::Impl::LocalFiniteElementFunction<LFE,cdim,ctype>;
+  using RefGeo = Dune::Impl::ReferenceElementGeometry<decltype(refElem)>;
+  using MappedGeometry = Dune::MappedGeometry<Mapping,RefGeo>;
+  auto mapping = Mapping{lfe,corners};
+  auto refGeo = RefGeo{refElem};
+  bool affine = lfe.localBasis().order() == 1 && refElem.template geometry<0>(0).affine();
+  auto mappedgeometry = MappedGeometry{mapping, refGeo, affine};
+  t.reset();
+  for (int i = 0; i < nIter; ++i)
+    pass &= benchmarkGeometry(mappedgeometry);
+  std::cout << "  MappedGeometry = " << t.elapsed() << "sec" << std::endl;
+
+  return pass;
+}
+
+
+template <class ctype>
+static bool benchmarkGeometries (int nIter)
+{
+  bool pass = true;
+
+  // pass &= benchmarkGeometries<ctype, 0, 0, Dune::GeometryTypes::simplex(0)>();
+
+  pass &= benchmarkGeometries<ctype, 1, Dune::GeometryTypes::simplex(1)>(nIter);
+  pass &= benchmarkGeometries<ctype, 2, Dune::GeometryTypes::simplex(1)>(nIter);
+  pass &= benchmarkGeometries<ctype, 3, Dune::GeometryTypes::simplex(1)>(nIter);
+  pass &= benchmarkGeometries<ctype, 4, Dune::GeometryTypes::simplex(1)>(nIter);
+
+  pass &= benchmarkGeometries<ctype, 1, Dune::GeometryTypes::cube(1)>(nIter);
+  pass &= benchmarkGeometries<ctype, 2, Dune::GeometryTypes::cube(1)>(nIter);
+  pass &= benchmarkGeometries<ctype, 3, Dune::GeometryTypes::cube(1)>(nIter);
+  pass &= benchmarkGeometries<ctype, 4, Dune::GeometryTypes::cube(1)>(nIter);
+
+  pass &= benchmarkGeometries<ctype, 2, Dune::GeometryTypes::simplex(2)>(nIter);
+  pass &= benchmarkGeometries<ctype, 3, Dune::GeometryTypes::simplex(2)>(nIter);
+  pass &= benchmarkGeometries<ctype, 4, Dune::GeometryTypes::simplex(2)>(nIter);
+
+  pass &= benchmarkGeometries<ctype, 2, Dune::GeometryTypes::cube(2)>(nIter);
+  pass &= benchmarkGeometries<ctype, 3, Dune::GeometryTypes::cube(2)>(nIter);
+  pass &= benchmarkGeometries<ctype, 4, Dune::GeometryTypes::cube(2)>(nIter);
+
+  pass &= benchmarkGeometries<ctype, 3, Dune::GeometryTypes::simplex(3)>(nIter);
+  pass &= benchmarkGeometries<ctype, 4, Dune::GeometryTypes::simplex(3)>(nIter);
+
+  pass &= benchmarkGeometries<ctype, 3, Dune::GeometryTypes::cube(3)>(nIter);
+  pass &= benchmarkGeometries<ctype, 4, Dune::GeometryTypes::cube(3)>(nIter);
+
+  // pass &= benchmarkGeometries<ctype, 3, 3, Dune::GeometryTypes::pyramid>(nIter);
+  // pass &= benchmarkGeometries<ctype, 3, 4, Dune::GeometryTypes::pyramid>(nIter);
+
+  // pass &= benchmarkGeometries<ctype, 3, 3, Dune::GeometryTypes::prism>(nIter);
+  // pass &= benchmarkGeometries<ctype, 3, 4, Dune::GeometryTypes::prism>(nIter);
+
+  pass &= benchmarkGeometries<ctype, 4, Dune::GeometryTypes::simplex(4)>(nIter);
+  pass &= benchmarkGeometries<ctype, 5, Dune::GeometryTypes::simplex(4)>(nIter);
+
+  pass &= benchmarkGeometries<ctype, 4, Dune::GeometryTypes::cube(4)>(nIter);
+  pass &= benchmarkGeometries<ctype, 5, Dune::GeometryTypes::cube(4)>(nIter);
+
+  return pass;
+}
+
+int main ( int argc, char **argv )
+{
+  bool pass = true;
+
+  int nIter = 100;
+  if (argc > 1)
+    nIter = std::atoi(argv[1]);
+
+  std::cout << ">>> Checking ctype = double" << std::endl;
+  pass &= benchmarkGeometries< double >(nIter);
+  std::cout << ">>> Checking ctype = float" << std::endl;
+  pass &= benchmarkGeometries< float >(nIter);
+
+  if (!pass)
+    std::cerr << "test failed!" << std::endl;
+
+  return (pass ? 0 : 1);
+}

dune/geometry/test/comparegeometries.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
+// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
+#ifndef DUNE_GEOMETRY_TEST_COMPAREGEOMETRIES_HH
+#define DUNE_GEOMETRY_TEST_COMPAREGEOMETRIES_HH
+
+#include <cmath>
+#include <iostream>
+#include <limits>
+#include <type_traits>
+
+namespace Dune {
+
+template <class R>
+R defaultTolerance ()
+{
+  using std::sqrt;
+  return sqrt(std::numeric_limits<R>::epsilon());
+}
+
+template <class Geometry1, class Geometry2,
+  class R = std::common_type_t<typename Geometry1::ctype, typename Geometry2::ctype>>
+bool compareGeometries (const Geometry1& geo1, const Geometry2& geo2,
+                        R tolerance = defaultTolerance<R>())
+{
+  if constexpr(Geometry1::mydimension != Geometry2::mydimension) {
+    std::cerr << "Error: Dimensions do not match." << std::endl;
+    return false;
+  }
+  else if constexpr(Geometry1::coorddimension != Geometry2::coorddimension) {
+    std::cerr << "Error: Coord-dimensions do not match." << std::endl;
+    return false;
+  }
+  else {
+    using std::abs;
+
+    if (geo1.affine() != geo2.affine()) {
+      std::cerr << "Error: Affine-property does not match." << std::endl;
+      return false;
+    }
+
+    if (geo1.type() != geo2.type()) {
+      std::cerr << "Error: GeometryType does not match." << std::endl;
+      return false;
+    }
+
+    if (geo1.corners() != geo2.corners()) {
+      std::cerr << "Error: Number of corners does not match." << std::endl;
+      return false;
+    }
+
+    for (int i = 0; i < geo1.corners(); ++i) {
+      if ((geo1.corner(i) - geo2.corner(i)).two_norm() > tolerance) {
+        std::cerr << "Error: Corner " << i << " does not match." << std::endl;
+        return false;
+      }
+    }
+
+    if ((geo1.center() - geo2.center()).two_norm() > tolerance) {
+      std::cerr << "Error: Center does not match." << std::endl;
+      return false;
+    }
+
+    if (abs(geo1.volume() - geo2.volume()) > tolerance) {
+      std::cerr << "Error: Volume does not match." << std::endl;
+      return false;
+    }
+
+    const auto& quadrature = Dune::QuadratureRules<R, Geometry1::mydimension>::rule(geo1.type(), 4);
+    for (auto&& [pos,weight] : quadrature)
+    {
+      if ((geo1.global(pos) - geo2.global(pos)).two_norm() > tolerance) {
+        std::cerr << "Error: global(" << pos << ") does not match." << std::endl;
+        return false;
+      }
+
+      if ((geo1.jacobian(pos) - geo2.jacobian(pos)).frobenius_norm() > tolerance) {
+        std::cerr << "Error: jacobian(" << pos << ") does not match." << std::endl;
+        return false;
+      }
+
+      if ((geo1.jacobianTransposed(pos) - geo2.jacobianTransposed(pos)).frobenius_norm() > tolerance) {
+        std::cerr << "Error: jacobianTransposed(" << pos << ") does not match." << std::endl;
+        return false;
+      }
+
+      if ((geo1.jacobianInverse(pos) - geo2.jacobianInverse(pos)).frobenius_norm() > tolerance) {
+        std::cerr << "Error: jacobianInverse(" << pos << ") does not match." << std::endl;
+        return false;
+      }
+
+      if ((geo1.jacobianInverseTransposed(pos) - geo2.jacobianInverseTransposed(pos)).frobenius_norm() > tolerance) {
+        std::cerr << "Error: jacobianInverse(" << pos << ") does not match." << std::endl;
+        return false;
+      }
+
+      if (abs(geo1.integrationElement(pos) - geo2.integrationElement(pos)) > tolerance) {
+        std::cerr << "Error: integrationElement(" << pos << ") does not match." << std::endl;
+        return false;
+      }
+    }
+
+    return true;
+  }
+}
+
+} // end namespace Dune
+
+#endif // DUNE_GEOMETRY_TEST_COMPAREGEOMETRIES_HH

dune/geometry/test/localfiniteelements.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
+// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
+#ifndef DUNE_GEOMETRY_TEST_LOCALFINITEELEMENT_HH
+#define DUNE_GEOMETRY_TEST_LOCALFINITEELEMENT_HH
+
+#include <algorithm>
+#include <vector>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/math.hh>
+
+namespace Dune::Impl {
+
+template <class D, class R, unsigned int dim>
+struct ScalarLocalBasisTraits
+{
+  using DomainFieldType = D;
+  using RangeFieldType = R;
+  using DomainType = FieldVector<D,dim>;
+  using RangeType = FieldVector<R,1>;
+  using JacobianType = FieldMatrix<R,1,dim>;
+
+  enum {
+    dimDomain = dim,
+    dimRange = 1
+  };
+};
+
+/// Lagrange shape functions of order 1 on the reference simplex
+template <class D, class R, unsigned int dim>
+class P1LocalBasis
+{
+public:
+  using Traits = ScalarLocalBasisTraits<D,R,dim>;
+
+  /// Number of shape functions
+  static constexpr unsigned int size () { return dim+1; }
+
+  /// Evaluate all shape functions
+  void evaluateFunction (const typename Traits::DomainType& x,
+                          std::vector<typename Traits::RangeType>& out) const
+  {
+    out.resize(size());
+    out[0] = 1.0;
+    for (unsigned int i=0; i<dim; i++)
+    {
+      out[0]  -= x[i];
+      out[i+1] = x[i];
+    }
+  }
+
+  /// Evaluate Jacobian of all shape functions
+  void evaluateJacobian (const typename Traits::DomainType& x,
+                          std::vector<typename Traits::JacobianType>& out) const
+  {
+    out.resize(size());
+    std::fill(out[0][0].begin(), out[0][0].end(), -1);
+
+    for (unsigned int i=0; i<dim; i++)
+      for (unsigned int j=0; j<dim; j++)
+        out[i+1][0][j] = (i==j);
+  }
+
+  /// Polynomial order of the shape functions
+  static constexpr unsigned int order () {  return 1; }
+};
+
+/// Lagrange shape functions of order 1 on the reference cube
+template <class D, class R, unsigned int dim>
+class Q1LocalBasis
+{
+public:
+  using Traits = ScalarLocalBasisTraits<D,R,dim>;
+
+  /// Number of shape functions
+  static constexpr unsigned int size () { return power(2, dim); }
+
+  /// Evaluate all shape functions
+  void evaluateFunction (const typename Traits::DomainType& x,
+                          std::vector<typename Traits::RangeType>& out) const
+  {
+    out.resize(size());
+    for (size_t i=0; i<size(); i++)
+    {
+      out[i] = 1;
+
+      for (unsigned int j=0; j<dim; j++)
+        // if j-th bit of i is set multiply with x[j], else with 1-x[j]
+        out[i] *= (i & (1<<j)) ? x[j] :  1-x[j];
+    }
+  }
+
+  /// Evaluate Jacobian of all shape functions
+  void evaluateJacobian (const typename Traits::DomainType& x,
+                          std::vector<typename Traits::JacobianType>& out) const
+  {
+    out.resize(size());
+
+    // Loop over all shape functions
+    for (unsigned int i=0; i<size(); i++)
+    {
+      // Loop over all coordinate directions
+      for (unsigned int j=0; j<dim; j++)
+      {
+        // Initialize: the overall expression is a product
+        // if j-th bit of i is set to 1, else -11
+        out[i][0][j] = (i & (1<<j)) ? 1 : -1;
+
+        for (unsigned int l=0; l<dim; l++)
+        {
+          if (j!=l)
+            // if l-th bit of i is set multiply with x[l], else with 1-x[l]
+            out[i][0][j] *= (i & (1<<l)) ? x[l] :  1-x[l];
+        }
+      }
+    }
+  }
+
+  /// Polynomial order of the shape functions
+  static constexpr unsigned int order () { return 1; }
+};
+
+
+template <class LB>
+class P1LocalInterpolation
+{
+public:
+  /// Evaluate a given function at the Lagrange nodes
+  template <class F, class C>
+  void interpolate (F f, std::vector<C>& out) const
+  {
+    constexpr auto dim = LB::Traits::dimDomain;
+    out.resize(LB::size());
+
+    // vertex 0
+    typename LB::Traits::DomainType x;
+    std::fill(x.begin(), x.end(), 0);
+    out[0] = f(x);
+
+    // remaining vertices
+    for (int i=0; i<dim; i++) {
+      for (int j=0; j<dim; j++)
+        x[j] = (i==j);
+
+      out[i+1] = f(x);
+    }
+  }
+};
+
+template <class LB>
+class Q1LocalInterpolation
+{
+public:
+  /// Evaluate a given function at the Lagrange nodes
+  template <class F, class C>
+  void interpolate (F f, std::vector<C>& out) const
+  {
+    constexpr auto dim = LB::Traits::dimDomain;
+    out.resize(LB::size());
+
+    typename LB::Traits::DomainType x;
+    for (unsigned int i=0; i<LB::size(); i++)
+    {
+      // Generate coordinate of the i-th corner of the reference cube
+      for (int j=0; j<dim; j++)
+        x[j] = (i & (1<<j)) ? 1.0 : 0.0;
+
+      out[i] = f(x);
+    }
+  }
+};
+
+
+/// Wrapper for local basis and local interpolation
+template <class LB, template <class> class LI>
+class LocalFiniteElement
+{
+public:
+  struct Traits
+  {
+    using LocalBasisType = LB;
+    using LocalInterpolationType = LI<LB>;
+  };
+
+  const auto& localBasis () const { return basis_; }
+  const auto& localInterpolation () const { return interpolation_; }
+
+  /// The number of shape functions
+  static constexpr std::size_t size () { return LB::size(); }
+
+private:
+  LB basis_;
+  LI<LB> interpolation_;
+};
+
+template <class D, class R, int d>
+using P1LocalFiniteElement = LocalFiniteElement<P1LocalBasis<D,R,d>, P1LocalInterpolation>;
+
+template <class D, class R, int d>
+using Q1LocalFiniteElement = LocalFiniteElement<Q1LocalBasis<D,R,d>, Q1LocalInterpolation>;
+
+
+
+template <class LFE, int cdim,
+          class R = typename LFE::Traits::LocalBasisType::Traits::RangeFieldType>
+class LocalFiniteElementFunction
+{
+  using LocalFiniteElement = LFE;
+  using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+  using LocalBasisRange = typename LocalBasis::Traits::RangeType;
+  using LocalBasisJacobian = typename LocalBasis::Traits::JacobianType;
+  using Domain = typename LocalBasis::Traits::DomainType;
+  using Range = FieldVector<R,cdim>;
+  using Jacobian = FieldMatrix<R,cdim,LocalBasis::Traits::dimDomain>;
+
+  static_assert(LocalBasis::Traits::dimRange == 1);
+
+public:
+  LocalFiniteElementFunction () = default;
+  LocalFiniteElementFunction (const LocalFiniteElement& lfe, std::vector<Range> coefficients)
+    : lfe_(lfe)
+    , coefficients_(std::move(coefficients))
+  {}
+
+  Range operator() (const Domain& local) const
+  {
+    thread_local std::vector<LocalBasisRange> shapeValues;
+    lfe_.localBasis().evaluateFunction(local, shapeValues);
+    assert(shapeValues.size() == coefficients_.size());
+    Range range(0);
+    for (std::size_t i = 0; i < shapeValues.size(); ++i)
+      range.axpy(shapeValues[i], coefficients_[i]);
+    return range;
+  }
+
+  friend auto derivative (const LocalFiniteElementFunction& f)
+  {
+    return [&lfe=f.lfe_,coefficients=f.coefficients_](const Domain& local) -> Jacobian
+    {
+      thread_local std::vector<LocalBasisJacobian> shapeJacobians;
+      lfe.localBasis().evaluateJacobian(local, shapeJacobians);
+      assert(shapeJacobians.size() == coefficients.size());
+      Jacobian jacobian(0);
+      for (std::size_t i = 0; i < shapeJacobians.size(); ++i) {
+        for (int j = 0; j < Jacobian::rows; ++j) {
+          shapeJacobians[i].umtv(coefficients[i][j], jacobian[j]);
+        }
+      }
+      return jacobian;
+    };
+  }
+
+private:
+  LocalFiniteElement lfe_{};
+  std::vector<Range> coefficients_{};
+};
+
+} // end namespace Dune::Impl
+
+#endif // DUNE_GEOMETRY_TEST_LOCALFINITEELEMENT_HH

dune/geometry/test/test-localfiniteelementgeometry.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
+// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
+#include <config.h>
+
+#include <limits>
+#include <cmath>
+#include <type_traits>
+#include <vector>
+
+#include <dune/geometry/localfiniteelementgeometry.hh>
+#include <dune/geometry/mappedgeometry.hh>
+#include <dune/geometry/multilineargeometry.hh>
+#include <dune/geometry/referenceelements.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/test/checkgeometry.hh>
+#include <dune/geometry/test/comparegeometries.hh>
+#include <dune/geometry/test/localfiniteelements.hh>
+#include <dune/geometry/test/referenceelementgeometry.hh>
+
+template <class ctype, int cdim, Dune::GeometryType::Id id>
+bool checkLocalFiniteElementGeometry ()
+{
+  bool pass = true;
+
+  constexpr auto gt = Dune::GeometryType{id};
+  auto refElem = Dune::referenceElement<ctype,gt.dim()>(gt);
+
+  using LFE = std::conditional_t<
+    gt.isSimplex(), Dune::Impl::P1LocalFiniteElement<ctype,ctype,gt.dim()>, std::conditional_t<
+    gt.isCube(),    Dune::Impl::Q1LocalFiniteElement<ctype,ctype,gt.dim()>, void>
+    >;
+  auto lfe = LFE{};
+
+  Dune::FieldMatrix<ctype,cdim,gt.dim()> A;
+  for (int i = 0; i < cdim; ++i)
+    for (int j = 0; j < int(gt.dim()); ++j)
+      A[i][j] = ctype(3-std::abs(i-j));
+
+  Dune::FieldVector<ctype,cdim> b;
+  for (std::size_t i = 0; i < cdim; ++i)
+    b[i] = ctype(i+1);
+
+  // mapping to generate coordinates from reference-element corners
+  auto f = [&](Dune::FieldVector<ctype,gt.dim()> const& x) {
+    Dune::FieldVector<ctype,cdim> y;
+    A.mv(x,y);
+    y += b;
+    return y;
+  };
+
+  auto corners = std::vector<Dune::FieldVector<ctype,cdim>>(refElem.size(gt.dim()));
+  for (int i = 0; i < refElem.size(gt.dim()); ++i)
+    corners[i] = f(refElem.position(i, gt.dim()));
+
+  using Geometry = Dune::LocalFiniteElementGeometry<LFE,cdim>;
+
+  // construct a geometry using given parametrization coefficients
+  auto geometry = Geometry{refElem, lfe, corners};
+  pass &= checkGeometry(geometry);
+
+  // construct a geometry using local interpolation
+  auto geometry2 = Geometry{refElem, lfe, f};
+  pass &= checkGeometry(geometry2);
+
+  // compare against MultiLinearGeometry
+  using MLGeometry = Dune::MultiLinearGeometry<ctype,gt.dim(),cdim>;
+  auto mlgeometry = MLGeometry{refElem, corners};
+  pass &= Dune::compareGeometries(geometry, mlgeometry);
+
+  // compare against a MappedGeometry
+  using Mapping = Dune::Impl::LocalFiniteElementFunction<LFE,cdim,ctype>;
+  using RefGeo = Dune::Impl::ReferenceElementGeometry<decltype(refElem)>;
+  using MappedGeometry = Dune::MappedGeometry<Mapping,RefGeo>;
+  auto mapping = Mapping{lfe,corners};
+  auto refGeo = RefGeo{refElem};
+  auto mappedgeometry = MappedGeometry{mapping, refGeo, true};
+  pass &= Dune::compareGeometries(geometry, mappedgeometry);
+
+  return pass;
+}
+
+
+template <class ctype>
+static bool checkLocalFiniteElementGeometry ()
+{
+  bool pass = true;
+
+  // pass &= checkLocalFiniteElementGeometry<ctype, 0, 0, Dune::GeometryTypes::simplex(0)>();
+
+  pass &= checkLocalFiniteElementGeometry<ctype, 1, Dune::GeometryTypes::simplex(1)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 2, Dune::GeometryTypes::simplex(1)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 3, Dune::GeometryTypes::simplex(1)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 4, Dune::GeometryTypes::simplex(1)>();
+
+  pass &= checkLocalFiniteElementGeometry<ctype, 1, Dune::GeometryTypes::cube(1)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 2, Dune::GeometryTypes::cube(1)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 3, Dune::GeometryTypes::cube(1)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 4, Dune::GeometryTypes::cube(1)>();
+
+  pass &= checkLocalFiniteElementGeometry<ctype, 2, Dune::GeometryTypes::simplex(2)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 3, Dune::GeometryTypes::simplex(2)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 4, Dune::GeometryTypes::simplex(2)>();
+
+  pass &= checkLocalFiniteElementGeometry<ctype, 2, Dune::GeometryTypes::cube(2)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 3, Dune::GeometryTypes::cube(2)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 4, Dune::GeometryTypes::cube(2)>();
+
+  pass &= checkLocalFiniteElementGeometry<ctype, 3, Dune::GeometryTypes::simplex(3)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 4, Dune::GeometryTypes::simplex(3)>();
+
+  pass &= checkLocalFiniteElementGeometry<ctype, 3, Dune::GeometryTypes::cube(3)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 4, Dune::GeometryTypes::cube(3)>();
+
+  // pass &= checkLocalFiniteElementGeometry<ctype, 3, 3, Dune::GeometryTypes::pyramid>();
+  // pass &= checkLocalFiniteElementGeometry<ctype, 3, 4, Dune::GeometryTypes::pyramid>();
+
+  // pass &= checkLocalFiniteElementGeometry<ctype, 3, 3, Dune::GeometryTypes::prism>();
+  // pass &= checkLocalFiniteElementGeometry<ctype, 3, 4, Dune::GeometryTypes::prism>();
+
+  pass &= checkLocalFiniteElementGeometry<ctype, 4, Dune::GeometryTypes::simplex(4)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 5, Dune::GeometryTypes::simplex(4)>();
+
+  pass &= checkLocalFiniteElementGeometry<ctype, 4, Dune::GeometryTypes::cube(4)>();
+  pass &= checkLocalFiniteElementGeometry<ctype, 5, Dune::GeometryTypes::cube(4)>();
+
+  return pass;
+}
+
+int main ( int argc, char **argv )
+{
+  bool pass = true;
+
+  std::cout << ">>> Checking ctype = double" << std::endl;
+  pass &= checkLocalFiniteElementGeometry< double >();
+  std::cout << ">>> Checking ctype = float" << std::endl;
+  pass &= checkLocalFiniteElementGeometry< float >();
+
+  if (!pass)
+    std::cerr << "test failed!" << std::endl;
+
+  return (pass ? 0 : 1);
+}