5

Having

mesh = UnitSquareMesh(2, 2)
Q = FunctionSpace(mesh, 'DG', 0)

this works as expected

R = mesh.ufl_cell().circumradius
R = project(R, Q)

but this fails

h = mesh.ufl_cell().max_facet_edge_length
h = project(h, Q)

when compiling generated code with error: 'facet_area' was not declared...

Does max_facet_edge_length return lengths of the longest edge of a cell as I expected or does it do something else? If so, is there some implementation of this? Or should I loop over cells and their edges to compute h manually?

Jan Blechta
  • 927
  • 5
  • 13

1 Answers1

2

Well, lengths of the longest edge of a cell h can be computed pythonically - i.e. not efficiently

h = Function(Q)
for c in cells(mesh):
    h.vector()[c.index()] = max(e.length() for e in edges(c))

For the clarification

mesh.ufl_cell().max_facet_edge_length

defines facet diameter (length of the longest edge) which can be evaluated on facets. For example interior penalty stabilization can be implemented this way

h_F = mesh.ufl_cell().max_facet_edge_length
n = FacetNormal(mesh)
alpha = Constant(0.01)
F += alpha('+')*h_F**2*inner(jump(grad(u), n), jump(grad(v), n))*dS

u being stabilize quantity, v its test function and alpha stabilization parameter here taken as constant. But may be for example scaled with convecting velocity for convection-dominated problems. See http://dx.doi.org/10.1090/S0025-5718-07-01951-5 and http://dx.doi.org/10.1007/s00211-007-0070-5

Jan Blechta
  • 927
  • 5
  • 13