fenics-adapter icon indicating copy to clipboard operation
fenics-adapter copied to clipboard

Compute CouplingExpression from point data via Minimization on P0DG function space

Open BenjaminRodenberg opened this issue 3 years ago • 9 comments

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:

https://github.com/precice/fenics-adapter/blob/be1b49de1adb3fef61b496ee84cd0183fed6a335/fenicsprecice/expression_core.py#L101-L114

BenjaminRodenberg avatar Mar 23 '21 14:03 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 avatar Apr 21 '21 11:04 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 https://github.com/precice/fenics-adapter/issues/119#issuecomment-824079207, creating the FunctionSpace still fails, but we get a few steps further. Thanks!

BenjaminRodenberg avatar Apr 21 '21 13:04 BenjaminRodenberg

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

ReubenHill avatar Apr 21 '21 13:04 ReubenHill

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.

BenjaminRodenberg avatar Apr 21 '21 14:04 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 avatar Apr 21 '21 14:04 ReubenHill

Relevant PR https://github.com/FEniCS/ufl/pull/30

ReubenHill avatar Apr 21 '21 14:04 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 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:

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.

BenjaminRodenberg avatar Apr 21 '21 15:04 BenjaminRodenberg

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.

ReubenHill avatar Apr 21 '21 16:04 ReubenHill

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!

BenjaminRodenberg avatar Apr 24 '21 11:04 BenjaminRodenberg