[Dune] Position of dofs for finite elements

[Timo Koch]

Dear Dune users & developers,

we are using dune-functions and the core modules in a simple FEM code. We want to strongly enforce Dirichlet boundary conditions.
The boundary conditions are specified locally for each element as a function of the element and the local position.

For globally defined analytical functions, we could use a global interpolate to get the coefficient / Dirichlet values. But how do we do this for “piece-wise”-defined boundary conditions?
We came up with something like this (restriction to simple Lagrange basis, nothing composite and assuming all dofs are Dirichlet for now):

Dune::BlockVector<Dune::FieldVector<double, numEq>> dirichletValues;
auto getDirichletValues = [&] (auto localIndex, const auto& localView, const auto& intersection)
{
    const auto& element = localView.element();

    // get the local coordinates of the dofs of the element
    std::vector<double> coords;
    const auto& interp = localView.tree().finiteElement().localInterpolation();
    interp.interpolate([&] (const LocalCoordinate& x) { return x[0]; }, coords);
    localDofCoords.resize(coords.size());
    for (unsigned int i = 0; i < coords.size(); ++i)
        localDofCoords[i][0] = coords[i];

    interp.interpolate([&] (const LocalCoordinate& x) { return x[1]; }, coords);
    for (unsigned int i = 0; i < coords.size(); ++i)
        localDofCoords[i][1] = coords[i];

    const auto index = localView.index(localIndex);
    for (unsigned int eqIdx = 0; eqIdx < numEq; ++eqIdx)
        dirichletValues[index] = problem.dirichlet(element, localDofCoords[localIndex]);
};

Dune::Functions::forEachBoundaryDOF(lagrangeBasis, getDirichletValues);

which is far from nice. In this way the coordinates are recomputed every time although it would be enough to compute them once per element type. Also we could figure out how to get the full position instead of only the x,y,z-components separately.
There seems to be no obvious way to get the local position of a DOF within an element, although this sounds like a very basic thing to do.
With the `localKey` it is possible to obtain the subEntity and its center position via the reference element. But this doesn’t (necessarily) correspond to the DOF position of course.

Is there some obvious way that we are missing on how to solve this problem without writing an excessive amount of code?
Do you have any suggestion or comments on how “piece-wise” defined strongly enforced Dirichlet boundary conditions could be realized?

Thanks a lot for your help
Timo

-- 
_________________________________________________

Timo Koch                                      phone: +49 711 685 64676
IWS, Universität Stuttgart                  fax:   +49 711 685 60430
Pfaffenwaldring 61         email: timo.koch at iws.uni-stuttgart.de
D-70569 Stuttgart             url: www.iws.uni-stuttgart.de/en/lh2/
_________________________________________________

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20200619/750d8a3d/attachment.htm>