5
$\begingroup$

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?

| cite | improve this question | | | | |
$\endgroup$
2
$\begingroup$

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

| cite | improve this answer | | | | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.