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

[Christian Engwer]

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