construct evaluation matrix

[nschloe]

Given a domain and a number of points within this domain, I would like to construct the small and fat evaluation matrix E which, given coefficient vector x returns E @ x, the values of the corresponding finite element function at the points. In FEniCS, I could do

  def _build_eval_matrix(V, points):
      """Build the sparse m-by-n matrix that maps a coefficient set for a
      function in V to the values of that function at m given points.
      """
      mesh = V.mesh()

      bbt = BoundingBoxTree()
      bbt.build(mesh)
      dofmap = V.dofmap()
      el = V.element()
      sdim = el.space_dimension()

      rows = []
      cols = []
      data = []
      for i, x in enumerate(points):
          cell_id = bbt.compute_first_entity_collision(Point(*x))
          cell = Cell(mesh, cell_id)
          coordinate_dofs = cell.get_vertex_coordinates()

          rows.append(numpy.full(sdim, i))
          cols.append(dofmap.cell_dofs(cell_id))

          v = el.evaluate_basis_all(x, coordinate_dofs, cell_id)
          data.append(v)

      rows = numpy.concatenate(rows)
      cols = numpy.concatenate(cols)
      data = numpy.concatenate(data)

      m = len(points)
      n = V.dim()
      matrix = sparse.csr_matrix((data, (rows, cols)), shape=(m, n))
      return matrix

The explicit loop doesn't hurt much since the number of points is usually small.

Any hints?

[gdmcbain]

This is very closely related to the recent feature ‘flip interpolator as probes’ #573. That defined a probes method of an InteriorBasis that took an ndarray of points and returned a function from V to ℝ^m. Internally, the returned function does

scikit-fem/skfem/assembly/basis/interior_basis.py

Line 147 in b29f231

return np.sum(phis * y[indices], 0)

where phis and indices are ndarrays known to the function (it's a closure).

Is this good enough? It'd be easy enough to rejig InteriorBasis.probes to return E if required, I think.

[kinnala]

It'd be easy enough to rejig InteriorBasis.probes to return E if required, I think.

Isn't E essentially phis inside probes with different shape and some additional zeros?

Should we have a separate method for returning E or is this sufficient?

[gdmcbain]

Isn't E essentially phis inside probes with different shape and some additional zeros?

Yes.

Alternatively one could define E as a LinearOperator to get the @ syntax, if that were desired; I'm O. K. with the function-call syntax, myself.

[nschloe]

I would need the linear operator and its transpose, so I'd be happy about a matrix of sorts. (Two operators would also work.)

[gdmcbain]

Interesting. I can't think what the transpose might be for, but it won't be hard to implement this so I'll have a go. Feel free to share the motivation if it's not too secret. I plan to soon be using probes or this or something like it heavily when I look into implementing the ‘method of characteristics’ #600.

[gdmcbain]

How's #601?

[nschloe]

Feel free to share the motivation if it's not too secret.

For sure! This is for smoothfit. The number 1 complaint is that people can't install it because of dolfin (understandable), so I thought I try out scikit-fem. In smoothfit, I have to solve a linear least-squares problem of the form

where E is the evaluation matrix and Delta is the Laplace with appropriate boundary conditions. In words: Get x to minimize the "bend" and minimize the deviation from prescribed values y0.

The can be solved with LSQR for which E.T is needed.

[nschloe]

How's #601?

LGTM.

[gdmcbain]

Well, smoothfit is very interesting indeed; thank you very much.

So in N-dimensions, eh? skfem.Mesh doesn't go beyond Mesh3D today, nor skfem.element beyond hexahedra and tetrahedra, but perhaps that has been merely for lack of a motivating application.

[nschloe]

I still get

AttributeError: 'InteriorBasis' object has no attribute 'probes_matrix'

when using skfem from master. Just probes works.

[nschloe]

Ah, I get it: probes is the probes_matrix.

[kinnala]

Yeah, probes was unreleased so we decided to use this approach.

[nschloe]

I've ported smoothfit to skfem and I'm pretty happy with the result. Finally I can actually run tests on gh-actions without having pull some crazy PPA tricks! Tests are still failing because probes hasn't been released yet, so I'm eagerly waiting for the next release. Thanks a lot so far!

[kinnala]

Alright, next release will be 3.0.0. I was expecting to further test and document some new features but we can always go for 3.0.x if any issues arise. I need to go through the issues and the changelog one more time and see if we are able to tag 3.0.0 in two weeks or so.