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?
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?
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).
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?
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]})")
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}")
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.
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:
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 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.