Efficiently Finding Maximum Function Value in a Mesh with FEniCSx

General
1

Dear FEniCSx Community,

I am currently working on a problem where I need to find the maximum value of a function in a triangular cell. My current approach involves sampling points (defined by the coordinates of the three vertices) and evaluating the function at these sampled points to determine the maximum value.

However, I am concerned that this method may miss the true maximum value due to the limited number of sample points, especially in regions where the function might have sharp peaks or variations.

Is there a more robust method within FEniCSx to accurately find the maximum value of a function over each mesh cell?

2

This will depend on the function space you are using, as you can probably exploit the structure of the basis functions of each cell to give you a good guess as to where the maximum of the function will be.

So the question is then, what is the function space you are using?

3

The functionspace I am using is linear lagrange space. According to your presentation, may i guess the max value exists in the vertex?

4

Since it is a linear space, the max value is definitely at a vertex. Thus you can use dofs = V.dofmap.cell_dofs(i) to get the indices of the vertices in a given cell (index i), and then access the appropriate u.x.array[dofs] (remember to unroll for block size if using a vector space).

5

Thanks for your affirmation. I think I can well access the value at vertex. However, what about lagrange 2-oder functionspace? The max value may be not at the vertex. Is there a method to get it accurately? Furthermore, is there a general method to find the max value of any type function in each cell?

My method is as follows:(dolfinx: v0.7.2)

def max_value_in_cell(function: fem.Function, cell_index: int, mesh: mesh.Mesh) -> float:
    cell_vertices = mesh.topology.connectivity(mesh.topology.dim, 0).links(cell_index)
    points = mesh.geometry.x[cell_vertices]
    num_samples = 10000
    sample_points = numpy.random.rand(num_samples, 2)
    u = sample_points[:, 0]
    v = sample_points[:, 1]
    mask = u + v > 1
    u[mask] = 1 - u[mask]
    v[mask] = 1 - v[mask]
    coords = (1 - u - v)[:, numpy.newaxis] * points[0] + \
          u[:, numpy.newaxis] * points[1] + \
          v[:, numpy.newaxis] * points[2]
    values = function.eval(coords, numpy.full(num_samples, cell_index))
    max_value = numpy.max(values)
    return max_value

Although sampling points are garanteed to locate in the cell, the efficiency seems to be poor and some true max value may be missed.

6

You would need to solve the optimization problem

\max_{\hat x} u(\hat x) = \sum_{i=0,..N} (\mathcal{F}^{-1}(\phi_i))(F(\hat x))c_i

where F and \mathcal{F} is described in DefElement: How to define a finite element, c_i are the local coefficients on the given cell, \hat x is the coordinates on the reference cell.

This function is also not correct, as you are mixing topology and geometry indices.

Here is what I would do in the case of P1:

from  mpi4py import MPI
import dolfinx
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dolfinx.cpp.mesh.CellType.triangle)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

u = dolfinx.fem.Function(V)
u.interpolate(lambda x: x[0] + 2*x[1] + np.cos(x[0]))

dofmap = V.dofmap
bs = dofmap.bs
dof_coords = V.tabulate_dof_coordinates()
num_cells_local = mesh.topology.index_map(mesh.topology.dim).size_local
for c in range(num_cells_local):
    cell_dofs = dofmap.cell_dofs(c)
    # Max per component in space
    cell_vals = u.x.array[cell_dofs]
    pos = np.argmax(cell_vals)
    print(f"Max value in cell {c} is {cell_vals[pos]} at {dof_coords[cell_dofs[pos]]} (other options: {cell_vals}, {dof_coords[cell_dofs]})")
7

My idea is to obtain the coordinates of the three vertices of the cell, not the degrees of freedom, because only in this way can we truly locate the range of the triangle. I think this function is correct if the search range is limited to sampling points.

The test is as follows:

from  mpi4py import MPI
import dolfinx
from dolfinx import fem, mesh
import numpy

def max_value_in_cell(function: fem.Function, cell_index: int, mesh) -> float:
    cell_vertices = mesh.topology.connectivity(mesh.topology.dim, 0).links(cell_index)
    points = mesh.geometry.x[cell_vertices]
    num_samples = 10000
    sample_points = numpy.random.rand(num_samples, 2)
    u = sample_points[:, 0]
    v = sample_points[:, 1]
    mask = u + v > 1
    u[mask] = 1 - u[mask]
    v[mask] = 1 - v[mask]
    coords = (1 - u - v)[:, numpy.newaxis] * points[0] + \
          u[:, numpy.newaxis] * points[1] + \
          v[:, numpy.newaxis] * points[2]
    values = function.eval(coords, numpy.full(num_samples, cell_index))
    max_pos = numpy.argmax(values)
    return values[max_pos], coords[max_pos], points

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dolfinx.cpp.mesh.CellType.triangle)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

u = dolfinx.fem.Function(V)
u.interpolate(lambda x: x[0] + 2*x[1])

num_cells_local = mesh.topology.index_map(mesh.topology.dim).size_local
for c in range(num_cells_local):
    value, coord, vertices_coords = max_value_in_cell(u, c, mesh)
    print(f"Max value in cell {c} is {value} at {coord}\n", 
          f"the vertices coordinates in cell {c} are \n{vertices_coords}")
8

What do you mean by this? If you use a P1 function space, the dofs are located at the vertices, whatever those vertices are are not affected by V.tabulate_dof_coordinates() which I just used to illustrate to you want the coordinates are.
The true value at the vertex is the value that is in u.x.array, u.eval will be affected by floating point accuracy in the push forward.

Secondly:

is not correct, as I mentioned earlier. See for instance: What is the difference between compute_incident_entities and entities_to_geometry? - #2 by dokken
and Getting coordinates of facets in dolfinx - #2 by dokken

9

Okay. What I mean is, if the function space is not “Lagrange 1” and the degrees of freedom are not located at the vertices of each cell, is there aother method to obtain the coordinates of the cell vertices?

Thank you for reminding me again. Is it correct to use the compute_incident_entities like this:

cell_vertices = compute_incident_entities(mesh.topology, cell_index, 2, 0)
points = mesh.geometry.x[cell_vertices]

Is the entities_to_geometry also possible like this:(dolfinx:v0.7.2)

cell_vertices = dolfinx.cpp.mesh.entities_to_geometry(mesh._cpp_object, 2, cell_index, False)
points = mesh.geometry.x[cell_vertices[0]]
10

There is a vertex to dof map:

which you can use to map vertices to dogs.

Combine that with entities_to_geoemtry, and you can get the values at vertices for any space that has dofs associated with vertices.

Note that for higher order spaces you are then excluding all other dofs, and you will likely not get the correct maximum value.

11

Hi dokken,

Apologies for getting back to this so late! I’ve been busy with other commitments recently, but now I’m finally revisiting this approach and trying to implement it properly. However, I’ve hit a bit of a roadblock and was hoping to ask for some additional guidance.

If I understand correctly, the idea is to first obtain the expression of the function on each cell in FEniCSx and then use a numerical optimization tool in Python to find the maximum value. Regarding obtaining the function expression, you previously shared a link about defining finite elements, which I believe is closely related to this process. However, the link now seems to be broken. Could you kindly share an updated link or point me to any other resources that might help?

Thank you very much for your assistance!

12

The new link is:

13

Thank you for sharing the updated link! It’s really helpful for understanding the finite element definitions in more detail.

I’ve been thinking about this and trying to start the implementation, but I find it really difficult to even begin the process—I feel stuck at the very first step. From my understanding, the first step involves combining the coefficients with the basis functions on the reference element to form the target function. Is that correct? I also tried searching the community for similar issues and came across some potentially helpful posts, like Accessing the basis function of a Finite Element, but I’m not sure if they are relevant or necessary for this optimization problem.

Could I ask for your further guidance on this? Thanks again for your patience and support—I really appreciate it.

14

Hello @dokken,

I am looking for the function values on the individual cells, in the similar manner as you shown above but on the vectorial function space
V = dolfinx.fem.functionspace(mesh, ("CG", 1, (mesh.geometry.dim,)))

I can manage to do this , if I need the individual components of the vectorial space seperately

      u0, V0_to_V = V.sub(0).collapse()
      v0, V1_to_V = V.sub(1).collapse()
      ux, uy = Function(u0), Function(v0)
      ux.x.array[:] = uh.x.array[V0_to_V]
      uy.x.array[:] = uh.x.array[V1_to_V]

But that’s not what I need. I want the vectorial values on individual cells.

Looking forward for your help. Thanks