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

[Lukas Riedel]

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