Skip to content

Add parametrized (curved) geometry

Simon Praetorius requested to merge feature/curved-geometries into master

Summary

This is the first version of a parametrized geometry that allows for curved elements. The parametrization is based on a local finite-element map. Note, the local finite-element is just a template parameter and thus not a dependency to dune-localfunctions.

Example 1

Parametrize a geometry by 2nd order Lagrange elements:

// construct a base grid e.g. with piecewise flat elements
YaspGrid<2> grid({1.0,1.0}, {10,10});

// geometry polynomial order
const int order = 2;
// define a local finite-element (from dune-localfunction)
using LocalFiniteElement = LagrangeCubeLocalFiniteElement<double,double,2,order>;
LocalFiniteElement localFE{};

for (const auto& e : elements(grid.leafGridView()))
{
  // projection from local coordinates in the reference element 
  // to 2d coordinates in the base grid to 3d global coordinates
  auto X = [geo=e.geometry()](const auto& local) -> FieldVector<double,3>
  {
    auto x = geo.global(local);
    return {x[0], x[1], std::sin(x[0] * x[1])};
  };

  // construct the parametrized geometry
  ParametrizedGeometry geometry{referenceElement(e), localFE, X};
}

Example 2

First construct a global grid parametrization in form of a vector holding the coordinates of all Lagrange points, then locally build the corresponding geometry. Note, the example uses a global basis from dune-functions, but it can be implemented also low-level

// construct a base grid e.g. with piecewise flat elements
YaspGrid<2> grid({1.0,1.0}, {10,10});

using namespace Dune::Functions::BasisFactory;
auto parametrizationBasis 
  = makeBasis(grid.leafGridView(), power<3>(lagrange<2>(), blockedInterleaved()));

// create avector that holds all the Lagrange points on the curved grid
std::vector<FieldVector<double,3>> points(parametrizationBasis.size());
Functions::interpolate(parametrizationBasis, points, 
  [](auto const& x) -> FieldVector<double,3> {
    return {x[0], x[1], std::sin(x[0] * x[1])};
  });

auto localView = parametrizationBasis.localView();
for (const auto& e : elements(grid.leafGridView()))
{
  localView.bind(e);
  auto const& node = localView.tree().child(0);

  // extract the points on the current elements
  std::vector<FieldVector<double,3>> localPoints(node.size());
  for (std::size_t i = 0; i < node.size(); ++i) {
    std::size_t idx = localView.index(node.localIndex(i))[0];
    localPoints[i] = points[idx];
  }

  // construct the parametrized geometry
  ParametrizedGeometry geometry{referenceElement(e), node.finiteElement(), localPoints};
}
Edited by Simon Praetorius

Merge request reports