Compute CouplingExpression from point data via Minimization on P0DG function space

[BenjaminRodenberg]

At FEniCS2021 @ReubenHill gave a talk on using point data in Firedrake. This might be an alternative approach to our CustomExpression and interpolation based approach, where the CustomExpression just wraps an interpolating function:

fenics-adapter/fenicsprecice/expression_core.py

Lines 101 to 114 in be1b49d

	def eval(self, value, x):
	"""
	Evaluates expression at x using self.interpolate(x) and stores result to value.
	Parameters
	----------
	value : double
	Buffer where result has to be returned to.
	x : double
	Coordinate where expression has to be evaluated.
	"""
	return_value = self.interpolate(x)
	for i in range(self._vals.ndim):
	value[i] = return_value[i]

[BenjaminRodenberg]

I investigated this issue a bit. Main problem with FEniCS is that it does not provide the VertexOnlyMesh needed for this purpose. We could try adapting the firedrake version of VertexOnlyMesh to also make it available in FEniCS. This is a bigger issue, though.

[BenjaminRodenberg]

I tried to create a mesh myself. I will provide the example code below. Note that the example is not working. I guess it is just not that simple to have a mesh without connectivity in FEniCS.

import fenics
import numpy as np
from fenics import Mesh, MeshEditor, FunctionSpace


# using https://bitbucket.org/fenics-project/dolfin/src/b55804ecca7d010ff976967af869571b56364975/dolfin/generation/IntervalMesh.cpp#lines-76:98 as template

N = 5  # we want to work with 5 vertices on the mesh
gdim = 2  # geometric dimension
tdim = 0  # topological dimension
vertices = np.random.rand(N, gdim)
mesh = Mesh()  # empty mesh
editor = MeshEditor()
editor.open(mesh, type="interval", tdim=tdim, gdim=gdim)  # here type="point" or type="vertex" would be a reasonable choice, but also leads to an error.
editor.init_vertices_global(N,N)

for i in range(N):
    editor.add_vertex(i, vertices[i,:])

editor.close()

V = FunctionSpace(mesh, "DG", 0)

Edit: Updated following advice in #119 (comment), creating the FunctionSpace still fails, but we get a few steps further. Thanks!

[ReubenHill]

Note: to create a mesh of disconnected vertices (i.e. a point cloud) you need tdim=0

[BenjaminRodenberg]

Neither

editor.open(mesh, type="point", tdim=tdim, gdim=gdim)

nor

editor.open(mesh, type="vertex", tdim=tdim, gdim=gdim)

work. Looking at the ufl package it looks like vertex is the correct type. However, the editor only accepts point.

[ReubenHill]

UFL should support vertex cells - at least I added such support to the github hosted FEniCS/ufl repository which I believe is the current UFL core and appears to be newer than the bitbucket one you pointed to. Not sure which you get when you import fenics as you do in your example mind you.

[ReubenHill]

Relevant PR FEniCS/ufl#30

[BenjaminRodenberg]

UFL should support vertex cells - at least I added such support to the github hosted FEniCS/ufl repository which I believe is the current UFL core and appears to be newer than the bitbucket one you pointed to. Not sure which you get when you import fenics as you do in your example mind you.

@ReubenHill On the ufl side vertex makes totally sense. I'm currently confused about the point which the editor expects. Looking at the dolfinx implementation I still see point, but not a vertex. I have the impression that dolfinx provides a point to ufl, but ufl expects a vertex:

• ufl.Cell gets a self.topology.cell_name()
• self.topology.cell_name() provides a string based on self.cell_type()
• cell_type() is provided by Topology
• Topology has a mesh::CellType
• But mesh::CellType is a point, not a vertex

Since I'm lacking a bit background in DolfinX and UFL, I'm not sure whether this is an issue or I'm just missing something.

[ReubenHill]

point and vertex are often used synonymously, but it seems that, in UFL at least, vertex is the term used to refer to a 0D simplex (see https://github.com/FEniCS/ufl/blob/master/ufl/cell.py).

It looks like the issue you are running into is trying to get the translation from the symbolic representation of a mesh of disconnected vertices (which you should be able to create in pure UFL as a ufl.mesh with a UFL ufl.FiniteElement(the_ufl_mesh, "DG", 0) defined on it) to the numerical implementation. The numerical implementation I have done is in firedrake and you're, I believe, trying to do the same thing in FEniCSx.

Assuming that FEniCSx makes use of the UFL domain and mesh abstractions for symbolically describing meshes then the symbolic bit is done and shouldn't need any UFL changes.

Looking at what I had to do in firedrake might be instructive to see what you need to change in FEniCSx to get this numerically represented but I honestly can't say for sure since I'm not a FEniCSx developer. Hopefully the above is helpful nonetheless.

[BenjaminRodenberg]

I asked a question about this issue in the FEniCS discourse forum (https://fenicsproject.discourse.group/t/unclear-how-to-use-mesheditor-with-pointcloud/5660). @ReubenHill Thanks for the explanation!