HowTo: Use FindPointsGSLIB for highorder interpolation
FindPointsGSLIB
provides a wrapper for highorder interpolation
via findpts
, a set of routines that were developed as
a part of the gatherscatter library, gslib.
While findpts
was originally developed for interpolation of grid functions in
H1
for meshes with quadrilateral or hexahedron elements, FindPointsGSLIB
also enables interpolation of functions in L2, H(div), H(curl)
on meshes with
triangle and tetrahedral elements.
The key steps of using FindPointsGSLIB
, as demonstrated in the
gslib miniapps
are:

First, setup the internal data structures required by the
gslib
library for the mesh of interest. This is done by using theFindPointsGSLIB::Setup(mesh)
method with the desiredmfem::Mesh
ormfem::ParMesh
. 
Next, use the
FindPointsGSLIB::FindPoints(xyz)
method with themfem::Vector xyz
of physicalspace coordinates where we seek to interpolate the desired grid function. The ordering of the coordinates inxyz
must be byvdim
, i.e. (x = x_{1},x_{2},x_{3}...,x_{N}, y_{1},y_{2},y_{3}...,y_{N}, z_{1},z_{2},z_{3}...,z_{N}).At this step,
findpts
determines the computational coordinates (q_{j} = {e_{j}, r_{j}, p_{j}}) for each point. These computational coordinates include the element number (e_{j} inmfem::Array<int> gsl_elem
) in which the point is found, the referencespace coordinates (r_{j} inmfem::Vector gsl_ref
) inside e_{j}, and the MPI rank that the element is partitioned on (p_{j} inmfem::Array<int> gsl_proc
).FindPoints
also returns a code (mfem::Array<int> gsl_code
) to indicate weather the point was found inside an element (gsl_code[j] = 0
), on the edge/face of an element (gsl_code[j] = 1
), or not found at all (gsl_code[j] = 2
) for the case when the point is located outside the mesh.Note that if a point (x_{j}) is located outside the mesh within a certain tolerance,
findpts
tries to find the closest location on the mesh surface (i.e.gsl_code[j] = 1
) and returns the distance (mfem::Vector gsl_dist
) between the sought point and the point found on the mesh surface. 
Finally, use
FindPointsGSLIB::Interpolate(u, ui)
to interpolate the desiredmfem::(Par)GridFunction u
at the physicalspace coordinates given byxyz
and return the interpolated values inmfem::Vector ui
.If
u
is inH1
, we usefindpts
for interpolation. Otherwise, we usefindpts
only for communicating computational coordinates of each point across MPI ranks, followed by MFEM's internal methods (mfem::GridFunction::GetValues
) for interpolation. 
Note, the
FindPointsGSLIB::FreeData()
method must be used before the program is terminated to free up the memory setup internally byfindpts
during thesetup
phase.
For convenience, FindPointsGSLIB
class provides methods such as
FindPointsGSLIB::Interpolate(mesh, xyz, u, ui)
which combines the three steps
described above (setup, finding the computational coordinates of the sought points, and
interpolation) into a single method. Please see the class definition
for more details.
Application of FindPointsGSLIB
The gslib
miniapps demonstrate several application of FindPointsGSLIB
:

findpts/pfindpts miniapps demonstrate highorder interpolation of a function in
H1
,L2
,H(div)
, orH(curl)
at an arbitrary set of points in physical space. 
fielddiff miniapp demonstrates comparison of grid functions defined on two different meshes.

fieldinterp miniapp demonstrates transfer of a grid function from one mesh on to another mesh.

schwarz_ex1/ex1p miniapp demonstrates use of overlapping Schwarz method to solve the Poisson problem in overlapping meshes. Here, we use
FindPointsGSLIB
to transfer solution between overlapping meshes to enforce Dirichlet conditions at the interdomain boundaries. 
cht Navier miniapp demonstrates how a conjugate heat transfer problem can be solve with the incompressible NavierStokes equations and the unsteady heat equation solved on different grids. Here,
FindPointsGSLIB
is used to transfer the solution from one mesh to another to couple the two PDEs.