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};
}