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?

1 Like
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?

1 Like
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

15

Dear @dokken, I have a possibly naive follow-up: is there a more efficient or vectorized way to compute the maximum value of the DOFs per cell across the entire mesh?

As far as I can tell, dofmap.cell_dofs() only accepts a single cell index at a time, which means I have to loop over all cells in Python. This works, but becomes a bottleneck for large meshes. Just wondering if thereā€™s a faster alternative ā€” maybe a way to get all cell DOFs in one shot, or use lower-level access?

16

dofmap.list will give you an adjacency list that you can access the raw data from (through .array and .offsets:

If you check that the offsets for each cell is constant, which is a 99% certainty, you can reshape dofmap.list.array as dofmap.list.array.reshape(num_cells_local+num_ghost_cells, -1) to get them all at once.

1 Like
17

Would it make sense to implement such a functionality in Fenicsx or do you think itā€™s not worth it because this can be done with e.g. Paraview?

18

In general, finding the maximum value within a cell is a quite expensive operation (when the basis function is not linear). However, it can be done, I just havenā€™t had the bandwidth to sit down and implement it (as it involves a Newton solver for each individual cell).

1 Like
19

Based on the suggestion from leader @dokken, I implemented a function that computes the maximum value of all DoFs in each cell. I tried both a vectorized version and a loop-based version of the function. To my surprise, I found that vectorization had almost no impact on the actual computation speed. Iā€™m sharing both versions of the code here, in case others are interested in exploring this further or spotting potential improvements.(FEniCSx version: 0.7.2)

def max_dofs_in_cells(fem_function: fem.Function, cell_indices, mesh):
    num_cells = mesh.topology.index_map(mesh.topology.dim).size_local + \
                mesh.topology.index_map(mesh.topology.dim).num_ghosts
    cell_dofs_array = fem_function.function_space.dofmap.list.reshape(num_cells, -1)
    target_cell_dofs = cell_dofs_array[cell_indices, :]
    max_dofs_in_cells = numpy.max(fem_function.x.array[target_cell_dofs], axis=1)
    return max_dofs_in_cells

def max_dofs_in_cells(fem_function: fem.Function, cell_indices, mesh):
    max_dofs_in_cells = numpy.zeros(len(cell_indices))
    for i, c in enumerate(cell_indices):
        cell_dofs = fem_function.function_space.dofmap.cell_dofs(c)
        max_dofs_in_cells[i] = numpy.max(fem_function.x.array[cell_dofs])
    return max_dofs_in_cells
1 Like
20

Hi @dokken, I just tested this on DOLFINx 0.9.0, and it seems that dofmap.list itself is already a 2D NumPy array with shape (num_cells_local + num_ghost_cells, ...), and it doesnā€™t have an .array attribute. So it looks like thereā€™s no need to access .array or .offsets, and also no need to reshape ā€” itā€™s already structured for per-cell DoF access.

I wonder if there might have been some mix-up between the Python and C++ layer APIs, or perhaps a change from earlier versions?

Could you confirm if thatā€™s the case? Just wanted to make sure Iā€™m not missing anything.