Make interpolator work for vector valued base functions

[emi0000]

With the interpolator function the value of an expansion at multiple x can be calculated according to

$f(x)=\sum\limits_{n=1}^{dofs} a_n \phi_n(x)$.

  poi=np.array([[0.3,0.4],[1.2,0.4],[1.2,1.6],[0.3,1.6]])
  tri=np.array([[0,1,2],[3,0,2]])
  mesh = fem.MeshTri(poi.T, tri.T)      
  e=fem.ElementTriP1()
  Vh = fem.Basis(mesh,e)
  xx=np.array([[0.8,0.5],[0.7,1.2]])
  aj=[1,0,0,0]
  interp = Vh.interpolator(aj)
  val = interp(xx)
  print(val)

However for vector valued base functions ( e.g. e=fem.ElementTriN1() ) this does not work

  ValueError: row, column, and data array must all be the same length

The error is due to the fact that the rows and cols are not correctly assembled.

This pull request fixes the error. In cell_basis.py the probes() function returns a matrix. The rows refer to the different points in the array xx, the cols refer to the dof-values.
In the fixed version the rows are simply extended by the additional components, so the first k (k number of points the sum has to calculate) rows belong to the first component, the next k rows to the second component and so on.
The following code is an example for Nedelec, first order in a triangle.

  import matplotlib.pyplot as plt
  import skfem as fem
  import numpy as np
  poi=np.array([[0.3,0.4],[1.2,0.4],[1.2,1.6],[0.3,1.6]])
  tri=np.array([[0,1,2],[3,0,2]])

  mesh = fem.MeshTri(poi.T, tri.T)      
  e=fem.ElementTriN1()
  Vh = fem.Basis(mesh,e)

  N=13
  aa=np.linspace(0.3,1.2,N)
  bb=np.linspace(0.4,1.6,N)
  xv, yv = np.meshgrid(aa, bb)
  xx=np.array([xv.flatten(),yv.flatten()])

  plt.figure(figsize=(8, 12))
  aj=Vh.dofs.N*[0]
  # change ii
  ii=1
  aj[ii]=1

  interp = Vh.interpolator(aj)
  val = interp(xx)
  qq=val.reshape(2,-1)

  plt.quiver(xx[0],xx[1],qq[0],qq[1])
  plt.triplot(mesh.p[0],mesh.p[1],mesh.t.T)
  plt.axis('equal')
  x0=Vh.doflocs[0][ii]
  y0=Vh.doflocs[1][ii]
  plt.plot(x0,y0,'or',ms=7,label='dof')
  plt.legend()
  plt.show()

I was not able to find a parameter that directly gives the number of components of a element or base, so it is determined in a first step by calculating the first basefunction at the origin of the local coordinate system and counting the number of values. Then the matrix is restructured using the phis and the number of components as described above.
The values returned after the matrix vector multiplication are therefore ordered in the same way. Currently they are not yet reshaped.

remarks

1. If there is a parameter like number_of_components I did not see
   comp = len(self.elem.lbasis(self.elem.dim*[0],0)[0])
   could be replaced

2. In the example above the returned array val is reshaped externally to separate the components. This could be done also in the interpolator function. A comp parameter would be helpful here too.

[kinnala]

There exists an open issue about this feature, #973.

There were some remaining open questions, e.g.,

• Shouldn't it support also matrix-valued elements (e.g., ElementTriHHJ1)?
• What should the output shape be for vector and matrix valued elements?

I didn't yet pursue analysing these questions in detail because the workaround to the issue was very simple: project the interesting component onto a scalar finite element basis and interpolate as follows:

Vh_scalar = Vh.with_element(ElementDG(ElementTriP1()))
ajx = Vh_scalar.project(Vh.interpolate(aj)[0])
valx = Vh_scalar.interpolator(ajx)
# ... use valx to interpolate the x-component

Note that I don't oppose solving this issue as presented by this PR, I just think those two questions above should be taken into account by the implementation, and omitting the matrix valued interpolator is also fine if there is a clear path towards how to implement it in future. I would also like to see some pros/cons type of analysis on possible alternative ideas concerning the output shape.

Comments on this PR:

• Tests are missing
• Avoid the use of try-except
• Take into account the discussion in Interpolator with ElementVector #973

[emi0000]

1. Structure of matrix for different base functions.
   -- One could have a scalar like probes() function which is called multiple times i.e. for each component of the base. I do not oversee all the implications of such a method. For the probes function one would probably need the component index as additional parameter. However when calling gbasis all components are given back, so why not using them directly?
   -- Therefore the other way would be to encode the different components into the rows, as I did in this PR. If the base function is matrix/tensor like the components are numbered in C- or fortran order. I would go for this method, as at least for the interpolator function the probes function is hidden, so the way the data are orderd for the calculation is not important for the user and invisible. That might be different if probes is used at other places too. There might be other reasons this is not a good idea I do not know, but as the call of probes so far does not work for vector like base functions I do not see any problems.

To be concrete, for three points and multiple components (vector like or matrix like)

point 1, first component --> $a_{11}$ $a_{12}$ $a_{13}$ $.....$ $a_{1N}$
point 2, first component --> $a_{21}$ $a_{22}$ $a_{23}$ $.....$ $a_{2N}$
point 3, first component --> $a_{31}$ $a_{32}$ $a_{33}$ $.....$ $a_{3N}$
$...............................................................................$
$...............................................................................$
point 1, last component --> $a_{K-2,1}$ $a_{K-2,2}$ $a_{K-2,3}$ $.....$ $a_{K-2N}$
point 2, last component --> $a_{K-1,1}$ $a_{K-1,2}$ $a_{K-1,3}$ $.....$ $a_{K-1,N}$
point 3, last component --> $a_{K,1}$ $a_{K,2}$ $a_{K,3}$ $.....$ $a_{K,N}$

2. output of the matrix multiplication
   For scalar base functions the output is as far as I see out.shape = (number of points,).
   -- For vector like base functions I would propose to return an array ( number of components, number of points). I would take this ordering, as in general points in skfem (e.g. mesh.p) are ordered that way. First all x values, then all y-values and so on. For example for ElementTriN1() the output would be (2, number of points)
   -- For matrix like base functions the ordering of the vectors above could be generalized. The general rule would be to use the structure of the base function, and each componental position contain the results for all points. For example for ElementTriHHJ1 the base function is a 2x2 matrix. So the output would be a (2,2,number of points) array. Another posibility would be to give back something like (number of components, number of points) in C/fortran order for matrix like base functions too.
   One could also think about letting the user decide about the output structure when calling the interpolator function by giving a tuple on input.

I think it would be necessary to define is a parameter/function in element or basis indicating the structure of the base function, e.g.
1 for scalar
(number of components) for vectors
(N,M) for matrix like base functions

In principle I did that by

  self.elem.lbasis(self.elem.dim*[0],0)

but maby that could be done routinely when calling e.g. fem.ElementTriN1() or fem.Basis(mesh,e)

Your are mentioning Basis._detect_tensor_order. Is there already such a parameter? I did not found it, to be honest.

BTW. What happens by doing ElementVector(ElementTriP1()) ?

[kinnala]

Basis._detect_tensor_order <- this does not yet exist but could be implemented. There are two options:

1. as you have done with the try-except check
2. add Element.order or similar property to each element one by one as you suggest

ElementVector(ElementTriP1()) <- this creates a vector-valued element with both x- and y-components represented by ElementTriP1(). It is a common element in classical mechanics, e.g, for linear elasticity and flow problems

[emi0000]

I introduced Basis._base_tensor_order
Its a tuple for tensor like base function and 1 for a scalar base function
e.g.
For ElementTriP1() Basis._base_tensor_order=1
For ElementTriN1() Basis._base_tensor_order=(2,)
For ElementTriHHJ1() Basis._base_tensor_order=(2,2)

This parameter is used to compute the number of components (1, 2, 4 for the examples above) and the matrix is assembled as described above.

on output of out = probes(x) @ y I reshaped the output
e.g.
For ElementTriN1()
out[0] : all x components of expansion for all points
out[1] : all y components of expansion for all points
For ElementHHJ1()
out[0,0] : component 1,1 of expansion for all points
out[0,1] : component 1,2 of expansion for all points
out[1,0] : component 2,1 of expansion for all points
out[1,1] : component 2,2 of expansion for all points

I also checked if ElementVector(ElementTriN1()) works
I did not yet write any testing, this is just to see if it makes sense to proceed.

[kinnala]

Sorry for the delay in responding - I was busy with teaching tasks this autumn. I think it makes sense. Are you still willing to push this forward and proceed with tests and such? I can also continue from your prototype at some point if you don't have the time.