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

[Lukas Riedel]

Hi Christian,

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:
- 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?
- I have a DG space. Which pattern assembly and residual assembly flags do I need? Is an ‘alpha_volume()’ implementation sufficient?
- 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?

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.
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