[dune-pdelab] Evaluate solution through a native DOF vector

[Lukas Riedel]

Hi Christian,

>> - I have a DG space. Which pattern assembly and residual assembly
>> flags do I need? Is an ‘alpha_volume()’ implementation sufficient?
> 
> a) you will not need any jacobian terms.
> b) in most cases the lambda_volume term should be sufficient.
> 
> regarding DG, you might encounter situations, where you don't want to
> evaluate the int_\Omega m u immediatelly, but if your m contains some
> gradient operator, you might want to use partial integration. In this
> case you will need additional skeleton terms. It would help, if you
> told us, what you measurment function looks like.

I’m thinking very simple here: My solution is a scalar quantity, and I want to measure it synthetically at a certain point in the domain. First point-like, then with a normalised Gaussian emulating an extended measurement volume.
Naively, I assume I can use the quadrature points in this case and evaluate the Gaussian, the shape functions and the integration factor like

template<typename EG, typename LFSV, typename R>
void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
    {
      const int order = lfsv.finiteElement().localBasis().order();
      const int intorder = intorderadd + quadrature_factor * order;
      auto gt = eg.geometry();
      for (const auto& it : quadratureRule(gt,intorder))
      {
        std::vector<Scalar> phi(lfsv.size());
        lfsv.finiteElement().localBasis().evaluateFunction(it.position(),phi);
        const RF factor = it.weight() * eg.geometry().integrationElement(it.position());
        const RF m = gaussian(eg.entity(),it.position());
        for (unsigned int i = 0; i<lfsv.size(); i++)
          r.accumulate(lfsv,i, m * phi[i] * factor);
      }
    }

>> - How would I implement a ‘narrow’ measurement function (or even a dirac function) in a local operator? How can I circumvent that it is zero at every quadrature point queried?
> 
> you can simply evaluate the shapefunctions at the point of your delta
> functions. There is no need to use numerical integration.

Perfect. I’ll try that.

>> On a side note, I achieved what I originally wanted by taking the first element of the MultiIndex I receive from the containerIndex() function. This works because my GFS is scalar.
> 
> a) it will not work anymore, when you use other measurment functions
> b) searching the pos is quite expensive. you might use the hierarchic
>   search, but still it might be expensive.

I have few measurement positions and I only need to assemble this ‘measurement operator’ once for every position. So the expensive search is not a large concern.

Cheers,
Lukas