fenics-adapter
fenics-adapter copied to clipboard
Compute CouplingExpression from point data via Minimization on P0DG function space
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
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.
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!
Note: to create a mesh of disconnected vertices (i.e. a point cloud) you need tdim=0
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
.
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.
Relevant PR https://github.com/FEniCS/ufl/pull/30
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 aself.topology.cell_name()
-
self.topology.cell_name()
provides a string based onself.cell_type()
-
cell_type()
is provided byTopology
-
Topology
has amesh::CellType
-
But
mesh::CellType
is apoint
, not avertex
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.
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.
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!