interpolation and projection

[emi0000]

Hi,

I tried to compare interpolation and projection in 1d and would be happy if a few questions could be clarified:

Interpolation is: $f(x)=\sum a_j \phi_j(x)$ with $a_j=\phi_j^*(f)$ (dual basis) and
Projection can be found here : https://scikit-fem.readthedocs.io/en/latest/howto.html

1. Is there a general interpolate function with a function as input like basis.interpolate(lambda x: ......). For Lagrange functions the coefficients are simply the values of the functions, but for more involved function (e.g. Nedelec) thats different, so a general interpolate function could then be useful

2. How can $f(x)$ i.e. the sum be calculated for different values of x which would be useful e.g. for P2 functions.

3. I was wondering why basis.project does not directly give back the $a_j$. I did a project myself and it was different by a factor of number_of_cells (see following code)

import matplotlib.pyplot as plt
import skfem as fem
import numpy as np

def MyProject(func,basis):
  # weak Form
  @fem.BilinearForm
  def K(u, v, w):
    return u*v

  @fem.LinearForm
  def D(v, w):
    x = w.x[0]  
    f = func(x)
    return f*v

  A = fem.asm(K, basis)
  b = fem.asm(D, basis)

  return fem.solve(A,b)  


mesh = fem.MeshLine(np.linspace(0, 1, 6))

e = fem.ElementLineP1()
basis = fem.Basis(mesh, e)

def function(x):
  val = 2*x*np.sin(2*np.pi*x)+3 
  return val

yp=MyProject(function,basis)
fp = basis.project(function)

t=np.linspace(0,1,1000)
plt.plot(t,function(t),label=r'$f(x)=2x\,sin(2\pi x)+3$')
plt.plot(mesh.p[0],function(mesh.p[0]),'-o',label='Interpolation')
plt.plot(mesh.p[0],yp,'-X',label='Projection')
plt.legend()
plt.show()

print(fp/yp)

[kinnala]

basis.interpolate is for interpolating at the integration points which is the most common interpolation operation required in FEM.

Interpolating at other points (defined by a global coordinate x) is cumbersome because it first requires finding in which element x lies in and only then interpolating using its local basis. This has been automated by the function basis.interpolator which returns a function that can be used to get the value at any global coordinate.

Your code looks like it should give the correct result but I cannot test it right now on mobile.

[emi0000]

Do you have a link or a small piece of code how to use basis.interpolator? thanks.

[kinnala]

Here I improvised:

In [1]: from skfem import *

In [2]: m = MeshTri().refined(3)

In [3]: basis = Basis(m, ElementTriP1())

In [4]: f = m.p[0]  # this represents the function f(x) = x in P1 basis

In [5]: interp = basis.interpolator(f)

In [6]: import numpy as np

In [7]: interp(np.array([[0.0], [0.0]]))
Out[7]: array([0.])

In [8]: interp(np.array([[1.0], [0.0]]))
Out[8]: array([1.])

[kinnala]

It can sometimes also help to read the docstring:

In [9]: help(basis.interpolator)

interpolator(y: numpy.ndarray) -> Callable[[numpy.ndarray], numpy.ndarray] method of skfem.assembly.basis.cell_basis.CellBasis instance
    Return a function handle, which can be used for finding
    values of the given solution vector `y` on given points.

Sometimes there are also examples in the docstring, e.g., for basis.get_dofs

[emi0000]

many thanks
I was trying to extend your code to Nedelec N1 functions, but that doesn't work

  import matplotlib.pyplot as plt
  import skfem as fem
  import numpy as np

  poi=np.array([[1.2,1.5],[2.1,1.6],[2.2,2.7]]).T
  tri=np.array([[0,1,2]]).T
  mesh = fem.MeshTri(poi, tri)
  plt.triplot(poi[0],poi[1],tri.T)
  pp=np.array([[1.6],[1.8]])

  e=fem.ElementTriP1()
  e=fem.ElementTriP2()
  e=fem.ElementTriN1()
  Vh = fem.Basis(mesh,e)

  # expansion coefficient
  a_j=np.arange(Vh.Nbfun)+0.11
  interp = Vh.interpolator(a_j)
  val = interp(pp)
  print("Phi(",pp,")=",val)

I looked into the code and figured out that the dimensions of the transformation matrix are not correctly mapped. The following part in cell_basis.py come up with the error

  return coo_matrix(
        (
            phis,
            (
                np.tile(np.arange(x.shape[1]), self.Nbfun),
                self.element_dofs[:, cells].flatten(),
            ),
        ),
        shape=(x.shape[1], self.N),
    )

For Nedelec N1 we have 3 base functions, each having two coordinates. So we have finally 6 values for the phis (len(phis) is indeed 6). The first two values belong to the first function, the next two to the second and so on, which can be verified with

  pts = Vh.mapping.invF(pp[:, :, np.newaxis], tind=[0])
  phis_=np.array([e.gbasis(Vh.mapping, pts, k, tind=[0])[0] for k in range(Vh.Nbfun)])

So each second value of the array phis should be moved to a second row in the coo_matrix and the shape must be adjusted accordingly.

Maybe the method is not intended for vector base functions? many thanks

[kinnala]

You are correct, interpolator does not work for other than scalar valued functions. You can project the components of a vector valued function onto a scalar FEM basis and then use interpolator.