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

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. — MiroK, May 16, 2014
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. — jenom, May 16, 2014

MiroK · Answer 1 · May 16, 2014

Best answer

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()

answered May 16, 2014 by MiroK FEniCS Expert (80,920 points)
selected May 16, 2014 by jenom

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

commented May 16, 2014 by jenom FEniCS Novice (690 points)

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

commented May 18, 2014 by Charles FEniCS User (4,220 points)

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

commented Jun 13, 2015 by yytamu FEniCS Novice (130 points)

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.

commented Jun 16, 2015 by jenom FEniCS Novice (690 points)

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

1 Answer