Mfem: Output/Probe matrix

[jandrej]

I am looking for a method on how to construct an output matrix such that

y = C u

where C is the output matrix and u is the state vector (solution of a linear system). C should describe a mapping from geometric points on the domain (not necessarily gridpoints) to function values y. So C is num_probes x num_states. Is there a convenient method on how to construct such an operator?

We don't necessarily need the C matrix, rather something that computes y from u.

In other words we need a measurement of a function at a geometric location. (In our case this is also a vector function, namely deformation in x y z).

My first idea would be to inherit from GridFunction (or use an existing method?) to hand over an array with coordinates and compute y through that. Since i cannot guarantee that the coordinates are located at quadrature points, i would have to interpolate somehow, thats where i am unsure on how to approach that.

Appreciate any help

[jandrej]

What i think i need is a "consistent interpolation".
Pseudocode ->

GridFunction gf
gf.evaluate_at(x, y)
-> find the element in which (x,y) lies
  -> use the polynomial basis to evaluate its value

It looks like GridFunction::GetValue/s implements the second part, but i still need to find the mapping from (x,y) to the cell.

Maybe this makes my problem more clear.

[tzanio]

This is not 100% robust, but you can use the FindPoints method of Mesh or ParMesh:

https://github.com/mfem/mfem/blob/ef6c759a4ef3d235abacc3153d9cfd92ced7114c/mesh/mesh.hpp#L1093-L1110

This is a recent addition via https://github.com/mfem/mfem/pull/240 and https://github.com/mfem/mfem/pull/265.

[jandrej]

Thanks @tzanio !
This works in combination with GridFunction::GetValue. I attached a small code snippet which demonstrates what i did. If this works with a vector valued problem on Monday i will close the issue.

   Vector point(2);
   point(0) = 0.5;
   point(1) = 0.5;
   DenseMatrix points(dim, 1);
   points.SetCol(0, point);

   Array<int> elem_ids;
   Array<IntegrationPoint> ips;

   mesh->FindPoints(points, elem_ids, ips);
   cout << x.GetValue(elem_ids[0], ips[0]) << endl;

[jandrej]

Works fine for our cases. Thanks again!