Min/Max on Cells and Edges of a Mesh for an Expression

Question

Hello,

I try to find the Min/Max value of some expression on a given mesh. I managed to do that for the min/max on cells as follows. First I interpolated the expression with some high degree k onto the mesh with a DG_k space. Then I ran over all the cells and took the max/min of the local dofs as the (approximate) max/min of the expression:

    for cell in cells( mesh ):
        i = cell.index()
        dofInt = dofmapInt.cell_dofs( i )
        dofMinMax = dofmapMinMax.cell_dofs( i )
        fMinVec[dofMinMax] = fIntVec[dofInt].min()
        fMaxVec[dofMinMax] = fIntVec[dofInt].max()

Here fMinVec is the vector of some DG0 function living on the mesh. Now I have the question:

How can I find the max/min over the edges?

I tried tabulate_facet_dofs in combination with MeshFunction but could not succeed. Any help would be very much appreciated,
Jo

Comments

Hi, I am curious if you insist on DG_k space? As you noticed, DG_k has all dofs bound to cell and none to lower order mesh entities. With CG_k the solution would be straightforward.

---

well actually that is not a constraint. in fact the expression is continuous so it makes more sense to use CG_k for the cells as well.

Answer 1

Hi, based on your comment, here's a solution that uses CG space for interpolation

from dolfin import *

mesh = UnitSquareMesh(20, 20)
V = FunctionSpace(mesh, 'CG', 4)
dm = V.dofmap()

f = Expression('sin(pi*x[0])*cos(pi*x[1])')
f = interpolate(f, V)
F = f.vector()

# Use edge function to store values
F_min = EdgeFunction('double', mesh, F.max())
F_max = EdgeFunction('double', mesh, F.min())

for cell in cells(mesh):
    cell_index = cell.index()
    cell_dofs = dm.cell_dofs(cell_index)

    for i, edge in enumerate(edges(cell)):
        # These are dofs internal to edge
        edge_dofs = list(cell_dofs[dm.tabulate_facet_dofs(i)])

        # Add vertex dofs
        vertex_dofs = [cell_dofs[dm.tabulate_facet_dofs(cell.index(v))][0]
                       for v in vertices(edge)]

        # Combine
        dofs = edge_dofs + vertex_dofs

        # Find min/max
        F_min[edge] = min(F[dofs].min(), F_min[edge])
        F_max[edge] = max(F[dofs].max(), F_max[edge])

plot(f)
plot(F_min, title='Min')
plot(F_max, title='Max')
interactive()

Comments

wow, thank you very much for your effort. That helped a lot.

---

Thanks for this, I had been doing this in post processing (paraview) - which annoyingly can't look at all of the DOFs.

---

Does this work when the code is running parallel? Thanks.

---

I'm far from being an expert on this topic but as far as I know the parallel code can still be quite brittle and parallelism is only applied in certain areas where it is mainly supported through the back end. So I would say in this particular situation it will just do nothing. I assumed that you are referring to the parallel options in FEniCS. If anybody knows better, please correct me.