[dune-pdelab] Evaluate solution through a native DOF vector
Christian Engwer
christian.engwer at uni-muenster.de
Thu Nov 24 12:42:01 CET 2016
Hi Lukas,
> thanks a lot for your reply! I see that your suggestion is what I eventually want to have, because I can use measurements that are not point-like.
> Since I’m not that much into numerics, I have a couple of questions:
as usual the answer starts with "it depends"...
> - 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.
> - 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.
> - Is it sufficient to place the local operator into a grid operator and call ‘residual()’? Does this give me the vector I’m looking for?
yes
> 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.
Ciao
Christian
> Here’s the code (not containing all of the typedefs):
>
> std::vector<RF> observation_operator (const typename Traits::Domain pos)
> {
> SolutionContainer u_raw = Dune::PDELab::Backend::native(u);
> std::vector<RF> res(u_raw.size());
>
> typedef LocalFunctionSpace<GFS> LFS;
> typedef LFSIndexCache<LFS> LFSCache;
> typedef typename U::template ConstLocalView<LFSCache> XView;
> typedef FiniteElementInterfaceSwitch<
> typename Dune::PDELab::LocalFunctionSpace<GFS>::Traits::FiniteElementType
> > FESwitch;
> LFS lfs(gfs);
> LFSCache lfs_cache(lfs);
> XView x_view(u);
> std::vector<typename Traits::Scalar> fem_factors(gfs.maxLocalSize()); // FEM coefficient container
>
> // find entity containing the coordinate
> auto it = gv.template begin<0>();
> auto pos_local = pos;
> for( ; it != gv.template end<0>(); ++it){
> const auto geo_type = it->type();
> const auto& ref = Dune::ReferenceElements<typename Traits::DF,Traits::dim>::general(geo_type);
> pos_local = it->geometry().local(pos);
> if(ref.checkInside(pos_local))
> break;
> }
> if(it == gv.template end<0>())
> DUNE_THROW(Dune::Exception,"Measurement position is not inside grid domain!");
>
> // evaluate local basis and find appropriate container index
> lfs.bind(*it);
> lfs_cache.update();
> x_view.bind(lfs_cache);
> FESwitch::basis(lfs.finiteElement()).evaluateFunction(pos_local,fem_factors);
> for(std::size_t i = 0; i<fem_factors.size(); ++i){
> const auto index = x_view.cache().containerIndex(i)[0];
> res[index] = fem_factors[i][0];
> }
>
> return res;
> }
>
>
> > On 22 Nov 2016, at 23:17, Christian Engwer <christian.engwer at uni-muenster.de> wrote:
> >
> > Dear Lukas,
> >
> > actually there is a much simpler (and judging from your code snippet
> > also more efficient) approach. Although I have to admit that it is not
> > obvious ;-)
> >
> > If I understand you correctly, what you want to do is to evaluate for
> > a given measurement function m (e.g. point evaluation, or
> > mean-value).
> >
> > Basically what you want to do is (given your solution u):
> > * compute v = \int_\Omega m(x) u(x) dx
> > * for a point evaluation m(x) is a dirac function, for the mean it is
> > 1/\int_\Omega 1 dx
> >
> > You can now the following. You write a new local operator which
> > contains the element wise weights, needed to compute the integral:
> > v = \sum_{E\in\Omega} v_E
> > with
> > v_E = \int_E m(x) u(x) dx
> > = \sum_i \sum_j w_i m(x_i) u_j \phi_j(x_i) detJ(x_i)
> > = \sum_j [\sum_j w_i m(x_i) \phi_j(x_i) detJ(x_i)] u_j
> > = \sum_j w_j u_j
> > => v = w * u
> >
> > You input data is the space of function you want to evaluate and you
> > map to the same space, as you want to store one coefficient for each
> > of the discrete coefficients of u. This local operator can now be
> > written without thinking about local or global numbering and without
> > special knowledge of the shape functions, as this is (as always)
> > handled by PDELab.
> >
> > Then you compute the weights w and the actual evaluation is computed
> > as the scalar product v = w*j.
> >
> > Hope this helps
> >
> > Ciao
> > Christian
> >
> >
> > On Tue, Nov 22, 2016 at 09:03:15PM +0100, Lukas Riedel wrote:
> >> Hi guys,
> >>
> >> for my current project I need to create a ‘measurement operator’ for a DOF vector. I.e., I want to multiply the native version (PDELab::Backend::native(u)) of this vector with another vector that contains the factors of the Finite Elements at the appropriate indices. Basically, I want to perform the same task as the evaluate() function of a DiscreteGridFunction outside the PDELab environment for a fixed position.
> >>
> >> To do this, I need to know which LocalFunctionSpace index of the given entity maps to which index of the native global DOF vector. I see that I can access the correct element of the global DOF vector directly via the ConstUncachedVectorView, and that I can get the DOFIndex from LFSIndexCache::containerIndex(i), but I cannot access the native container with this data type. Is there a mapping from a local index to the index of the native DOF vector?
> >>
> >> Here’s an excerpt from my code:
> >>
> >> auto u_native = Dune::PDELab::Backend::native(u); // native DOF vector
> >> std::vector<typename Traits::RangeType> evaluator(u_native.size()); // evaluation vector
> >>
> >> typedef LocalFunctionSpace<GFS> LFS;
> >> typedef LFSIndexCache<LFS> LFSCache;
> >> typedef typename U::template ConstLocalView<LFSCache> XView;
> >> LFS lfs(gfs);
> >> LFSCache lfs_cache(lfs);
> >> XView x_view(u);
> >> std::vector<typename Traits::RangeFieldType> xl(gfs.maxLocalSize());
> >> std::vector<typename Traits::RangeType> yb(gfs.maxLocalSize());
> >> typedef FiniteElementInterfaceSwitch<
> >> typename Dune::PDELab::LocalFunctionSpace<GFS>::Traits::FiniteElementType
> >> > FESwitch;
> >>
> >> auto it = gv.template begin<0>(); // take the first entity for now
> >> auto pos_global = it->geometry().center(); // take center position for now
> >> auto pos_local = it->geometry().local(pos_global);
> >> lfs.bind(*it);
> >> lfs_cache.update();
> >> x_view.bind(lfs_cache);
> >> x_view.read(xl);
> >> FESwitch::basis(lfs.finiteElement()).evaluateFunction(pos_local,yb);
> >>
> >> for(std::size_t i = 0; i<xl.size(); ++i){
> >> if(xl[i] == x_view[i]){ // accessing the correct element of the solution vector
> >> std::cout << “It’s a match!” << std::endl;
> >> }
> >> std::cout << “This is a strange object, but its printout contains the index I’m looking for: " << x_view.cache().containerIndex(i) << std::endl;
> >> // u_native[x_view.cache().containerIndex(i)]; // this doesn’t work
> >>
> >> /* what I want to have:
> >> auto global_index = x_view.cache().containerIndex(i) // ???
> >> evaluator[global_index] = yb[i];
> >> */
> >> }
> >>
> >>
> >> Thank you in advance
> >> Cheers,
> >> Lukas
> >> _______________________________________________
> >> dune-pdelab mailing list
> >> dune-pdelab at dune-project.org
> >> http://lists.dune-project.org/mailman/listinfo/dune-pdelab
More information about the dune-pdelab
mailing list