Fix facet and cell averages

[jorgensd]

Currently we do not support FacetAvg or CellAvg in ufl forms, as the following mwe illustrates:

import ufl

cell = ufl.triangle
c_el = ufl.VectorElement("Lagrange", cell, 1)
mesh = ufl.Mesh(c_el)
el = ufl.FiniteElement("Lagrange", cell, 1)
V = ufl.FunctionSpace(mesh, el)
u = ufl.Coefficient(V)

a = ufl.cell_avg(u) * ufl.dx
L = ufl.facet_avg(u) * ufl.ds
forms = [a, L]

Traceback (most recent call last):
  File "/usr/lib/python3.10/runpy.py", line 196, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/usr/lib/python3.10/runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "/root/shared/ffcx/ffcx/__main__.py", line 12, in <module>
    sys.exit(main())
  File "/root/shared/ffcx/ffcx/main.py", line 68, in main
    code_h, code_c = compiler.compile_ufl_objects(
  File "/root/shared/ffcx/ffcx/compiler.py", line 102, in compile_ufl_objects
    ir = compute_ir(analysis, object_names, prefix, parameters, visualise)
  File "/root/shared/ffcx/ffcx/ir/representation.py", line 197, in compute_ir
    irs = [
  File "/root/shared/ffcx/ffcx/ir/representation.py", line 198, in <listcomp>
    _compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
  File "/root/shared/ffcx/ffcx/ir/representation.py", line 487, in _compute_integral_ir
    integral_ir = compute_integral_ir(itg_data.domain.ufl_cell(), itg_data.integral_type,
  File "/root/shared/ffcx/ffcx/ir/integral.py", line 71, in compute_integral_ir
    S = build_scalar_graph(expression)
  File "/root/shared/ffcx/ffcx/ir/analysis/graph.py", line 79, in build_scalar_graph
    scalar_expressions = rebuild_with_scalar_subexpressions(G)
  File "/root/shared/ffcx/ffcx/ir/analysis/graph.py", line 179, in rebuild_with_scalar_subexpressions
    ws = reconstruct(expr, wops)
  File "/root/shared/ffcx/ffcx/ir/analysis/reconstruct.py", line 168, in reconstruct
    raise RuntimeError("Not expecting expression of type %s in here." % type(o))
RuntimeError: Not expecting expression of type <class 'ufl.averaging.CellAvg'> in here.

This is due to: https://bitbucket.org/fenics-project/ufl/pull-requests/89 which never got fixed in old DOLFIN.

[michalhabera]

Substitution CellAvg(u) --> u / CellVolume does not look correct to me, these things have also different units (u-units vs. u-units/volume). And it is also wrong for constant functions, since CellAvg(u) = u for u in DG0.

[jorgensd]

Substitution CellAvg(u) --> u / CellVolume does not look correct to me, these things have also different units (u-units vs. u-units/volume). And it is also wrong for constant functions, since CellAvg(u) = u for u in DG0.

Im then not sure what a cell_avg should be. As far as I can tell this is what they do in tsfc (https://github.com/firedrakeproject/tsfc/pull/158)

I see that this is wrong for DG 0. I guess what we actually want is cell_avg(var)*dx->(var/CellVolume*dx)*dx

@michalhabera do we support such operations in ffcx?

[michalhabera]

Yes, the correct implementation of entity averages requires an integral within integral. Unfortunately, FFCx is not flexible enough for this at the moment.

This should become easier once we have separate codepaths for handling of expressions (integrands) and their evaluation at quadrature points to represent integrals. This is a longer standing refactoring issue.

On 25. Jun 2022, at 09:47, Jørgen Schartum Dokken @.***> wrote:

 Substitution CellAvg(u) --> u / CellVolume does not look correct to me, these things have also different units (u-units vs. u-units/volume). And it is also wrong for constant functions, since CellAvg(u) = u for u in DG0.

Im then not sure what a cell_avg should be. As far as I can tell this is what they do in tsfc (firedrakeproject/tsfc#158)

I see that this is wrong for DG 0. I guess what we actually want is cell_avg(var)dx->(var/CellVolumedx)*dx

@michalhabera do we support such operations in ffcx?

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.

[wence-]

As far as I can tell this is what they do in tsfc (firedrakeproject/tsfc#158)

TSFC implements CellAvg(expr)*dx (and FacetAvg) as cellwise_integral(expr/CellVolume *dx) * dx as @michalhabera is advocating for