Evaluation of UFL expression at mesh points

wvanharten · December 23, 2021, 4:09pm

Hi all,

Background
I want to perform a coordinate transform of my PDE, the Helmholtz equation. This coordinate transform is prescribed by another PDE which I solve in disjoint regions. Because this is solved in disjoint regions it is not continuously differentiable. To finally transform the Helmholtz equation I introduce, among other terms, a discontinuous Jacobian.

To be able to use the calculated transformation in the weak formulation of the Helmholtz equation, I have to map the values of the computed coordinate transformation onto a function that is defined on the mesh of the complete space.

To obtain the values at the vertices of the mesh of the solution x_transform I can do

values_x_transform = x_transform.x.array

However, I cannot do the same to get the values at the vertices of the derivative of x_transform, Dx(x_transform, 0). I would like to do something like

x_transform_derivative = Dx(x_transform, 0)
values_x_transform_derivative = x_transform_derivative.x.array

But this is not possible. I currently do this by solving an auxillary equation:

u = TrialFunction(V)
v = TestFunction(V)

a = inner(u, v) * dx
L = inner(x_transform_derivative, v) * dx

uh = Function(V)

problem = fem.LinearProblem(a, L, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.solve()

values = uh.x.array  # Extract sought vector

but this sounds rather wasteful on the computer resource front because a complete system has to be solved for this.

Question

Can I evaluate an UFL expression g at the mesh points of the underlying mesh? I have tried g(mesh.geometry.x) but this results in an infinite recursion loop.

A minimal working example (MWE) would be

import dolfinx
from dolfinx import FunctionSpace, Function, fem
from mpi4py import MPI
from dolfinx.generation import UnitSquareMesh
from ufl import TrialFunction, TestFunction, inner, dx, Dx

mesh = UnitSquareMesh(MPI.COMM_WORLD, 8, 8, dolfinx.cpp.mesh.CellType.quadrilateral)
V = FunctionSpace(mesh, ("CG", 1))

f = Function(V)
f.interpolate(lambda x: (x[0] ** 2) * (1 + 1j))  # An interesting nonzero example

g = Dx(f, 0)  # Some ufl expression we want to evaluate

g(mesh.geometry.x)  # Segfault due to infinite recursion:

'''
Infinite recursion encountered:
(...)
  File "/../ufl/core/expr.py", line 318 in __float__
  File "/../ufl/core/terminal.py", line 50 in evaluate
  File "/../ufl/exproperators.py", line 320 in _eval
  File "/../ufl/exproperators.py", line 330 in _call
  File "/../ufl/core/expr.py", line 313 in _ufl_evaluate_scalar_
  File "/../ufl/core/expr.py", line 318 in __float__
  File "/../ufl/core/terminal.py", line 50 in evaluate
  File "/../ufl/exproperators.py", line 320 in _eval
  File "/../ufl/exproperators.py", line 330 in _call
  File "/../ufl/core/expr.py", line 313 in _ufl_evaluate_scalar_
(...)
'''

# ---- Known solution, presumably rather slow because a full linear system needs to be solved ---- #
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u, v) * dx
L = inner(g, v) * dx
uh = Function(V)
problem = fem.LinearProblem(a, L, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.solve()
values = uh.x.array  # Extract sought vector

Can I extract the values of the UFL expression at the mesh points without solving an auxillary problem as I have done above?

More information:
I use the Docker image build from the 0.3.0 release with complex mode enabled.

Many thanks in advance!

-Wouter

dokken · December 23, 2021, 5:10pm

You can use dolfinx.Expression as shown in

github.com

FEniCS/dolfinx/blob/114c268aab752b6669abc2718837116cec90dba6/python/test/unit/fem/test_expression.py#L70

    
      
          cell_geometry = np.zeros((num_dofs_x, gdim), dtype=np.float64)
          points_ref = np.zeros((len(local_map), tdim))
          
          
# Map cells on process back to reference element
          for i, cell in enumerate(closest_cell):
              cell_geometry[:] = x[x_dofs[cell], :gdim]
              point_ref = mesh.geometry.cmap.pull_back(points[local_map[i]][:gdim].reshape(1, -1), cell_geometry)
              points_ref[i] = point_ref
          
          
# Eval using Expression
          expr = Expression(grad_u, points_ref)
          grad_u_at_x = expr.eval(closest_cell).reshape(len(closest_cell), points_ref.shape[0], gdim)
          
          
# Compare solutions
          for i, cell in enumerate(closest_cell):
              point = points[local_map[i]]
              assert np.allclose(grad_u_at_x[i, i], gradf(point))

Some nicer front-end UI will hopefully be merged soon with: Interpolation of ufl expressions by jorgensd · Pull Request #1875 · FEniCS/dolfinx · GitHub

wvanharten · December 24, 2021, 9:36am

Dear Dokken,

Thank you for your quick reply. If I’m not mistaken the Expression.eval(...) will produce a full matrix solution; the full Expression evaluated at all the points points. Finally, in the example you provided, the diagonal is extracted.

However, for large number of points to evaluate (in my case, all the vertices), the memory requirements scale with the square of the number of vertices. Therefore, using the proposed method is not really feasible.

Is this solved with the (soon to be merged) interface?

Thank you!

And if you celebrate christmas, happy holidays

dokken · December 24, 2021, 2:26pm

I do not follow your logic here. Say we have a triangle, and you want to evalutate at every vertex of every cell, you can call expression with the corners of the reference element, i.e, (0,0),(1,0),(0,1).
This would scale as 3*num_cells.

Expression.Eval produces a matrix of size (num_cells, num_vertices*size_of_expression), where size_of_expression is 1 if the expression is a scalar, gdim if the expression is a vector, gdim^2 if the expression is a tensor.

wvanharten · January 3, 2022, 1:30pm

Hi,

First of all, all the best in '22!

Let me try again to formulate my question.

Suppose I have a mesh with 121 vertices and 200 faces:

with some function f defined on the mesh.

In order to find the value of g=Dx(f,1) at all the 121 different vertices I would construct an expression with the size size_of_expression=121 (or some different, more optimal, subset of the faces covering all vertices). This is also what happens in the test I was pointed to: test_expression.py
Now, Expression.Eval will evaluate to num_cells * size_of_expression \propto num_cels^2.

Should I maybe construct an different Expression instance per face I want to evaluate? Or shouldn’t I construct an Expression of size 121?

Many Thanks in advance!

Wouter

dokken · January 3, 2022, 2:16pm

Please note that Dx(f, 1) is not uniquely defined at vertices (when f is a CG 1 function). Thus there is ambiguity in what the value should be at a vertex.

What you can do if you insist of evaluating the expression at the vertices, is that you:

Interpolate the ufl expression into an appropriate space, as for instance shown here: dolfinx/test_interpolation.py at c9f3c459a75f02abcf255a9dd05185350eb1edde · FEniCS/dolfinx · GitHub
Use the point evaluation function for the vertices, as shown here: Implementation — FEniCSx tutorial

However, note that you can send in the cells you want to use the dolfinx.Expression at, which you can reduce to a subset of cells spanning all the unique vertices, reducing the size of the output array from cells * num_vertices_per_cell=3*num_cells in your example above, to 3*subset_of_cells, which will depend on which cells you choose to cover all vertices.

Topic		Replies	Views
Shape Derivatives in UFL dolfin-adjoint	21	2724	June 10, 2020
Evaluating Finite Element Function General	10	572	February 1, 2024
Evaluate function on whole mesh General	2	282	October 11, 2023
Updating dof of dolfinx.fem.function.Function using Midpoint time stepping General	2	343	May 31, 2023
Solving ufl variable by equation solve dolfinx	0	130	February 25, 2024

Evaluation of UFL expression at mesh points

Related topics