# How to apply edge size and normal gradient

**URL:** <https://fenicsproject.discourse.group/t/how-to-apply-edge-size-and-normal-gradient/13420>\
**Category:** General\
**Tags:** dolfinx\
**Created:** [January 26, 2024, 2:21pm UTC](https://fenicsproject.discourse.group/t/how-to-apply-edge-size-and-normal-gradient/13420 "2024-01-26T14:21:37Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![WjjSteve](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wjjsteve/32/4389_2.png) [@WjjSteve](https://fenicsproject.discourse.group/u/WjjSteve)\
**Post date:** [January 26, 2024, 2:21pm UTC](https://fenicsproject.discourse.group/t/how-to-apply-edge-size-and-normal-gradient/13420/1 "2024-01-26T14:21:37Z")

</div>

Hello,

I want to ask some quesions about basix and ufl.

If I have the following term in my formulation:

![Screenshot from 2024-01-26 15-08-02](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/ffca34aa460d9928a4a6be3373674409847649e6.png)  
here I want to compute the intergral over **edges** , in the above formulation, h is the **edge length** , and jump(\partial\_n u) means the **jump of the normal gradient across the edge**. Since the above terms are not commonly used, I would like the ask:

1: Is it possible to apply the edge size, I know that ufl has some function like CellDiameter, but that is for the cell in stead of edge. (My mesh is not in 1D)  
2: How can I compute the normal gradient across the edge

Thanks very much!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 26, 2024, 6:33pm UTC](https://fenicsproject.discourse.group/t/how-to-apply-edge-size-and-normal-gradient/13420/2 "2024-01-26T18:33:35Z")

</div>

Is your mesh 3D, and you want to compute a 1D integral?  
Or is your mesh 2D, and you want to compute a jump over an edge (which is equivalent to a facet)?

If you do want 3D - 1D coupling, how is the jump defined? A normal jump contains to cells, but an edge can be connected to an arbitrary set of cells.

3D - 1D is something @cdaversin and I hope to look into later this spring.

---

<div class="post-metadata">

**Author:** ![WjjSteve](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wjjsteve/32/4389_2.png) [@WjjSteve](https://fenicsproject.discourse.group/u/WjjSteve)\
**Post date:** [January 26, 2024, 10:45pm UTC](https://fenicsproject.discourse.group/t/how-to-apply-edge-size-and-normal-gradient/13420/3 "2024-01-26T22:45:35Z")

</div>

Thanks for quick reply.

Currently I only need work with the 2D case, as you said I need the jump over the facet. But I don’t know how to calculate the normal gradient, it seems like we can only use the partial gradient like w.dx

Also, when I calculate the intergral over edge, I need a constant which is the edge size, can I get property of the facet instead of the cell?

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [January 27, 2024, 5:07am UTC](https://fenicsproject.discourse.group/t/how-to-apply-edge-size-and-normal-gradient/13420/4 "2024-01-27T05:07:59Z")

</div>

It’s difficult to understand exactly what you require since you may be confusing terminology of face, facet and edge. See for example [https://arxiv.org/pdf/1205.3081.pdf](https://arxiv.org/pdf/1205.3081.pdf) Table 1.

Regarding normal fluxes, this is straightforward to achieve. Consider a cell \kappa in a finite element mesh where we define

\partial\_{n\_\kappa} u = \nabla u \cdot n\_\kappa

which is represented in UFL via something along the lines of

```python
mesh = ...
u = ...
n = ufl.FacetNormal(mesh)
flux = ufl.dot(ufl.grad(u), n)

```

Similarly you can use the [`ufl.jump`](https://fenics.readthedocs.io/projects/ufl/en/latest/manual/form_language.html#jump-jump-v) function as required for the jump over facets.

In interior penalty methods, typically the measure h is not strictly a facet area (cf. “edge length” in the 2D case); rather a measure of cell size. Nevertheless UFL offers [`ufl.FacetArea`](https://github.com/FEniCS/ufl/blob/main/ufl/geometry.py#L764) which may suit your needs.

As an example you could examine the [biharmonic demo](https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_biharmonic.py#L204-L208) which employs a C^0-SIPG formulation employing these operators in 2D.

However, if you require edge (topology dimension D-2) measures in 3D computations, things get a lot more complicated as @dokken has already mentioned.

---

<div class="post-metadata">

**Author:** ![WjjSteve](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wjjsteve/32/4389_2.png) [@WjjSteve](https://fenicsproject.discourse.group/u/WjjSteve)\
**Post date:** [January 27, 2024, 12:12pm UTC](https://fenicsproject.discourse.group/t/how-to-apply-edge-size-and-normal-gradient/13420/5 "2024-01-27T12:12:15Z")

</div>

That is exactly what I need, thanks very much!

---

<div class="post-metadata">

**Author:** ![WjjSteve](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/wjjsteve/32/4389_2.png) [@WjjSteve](https://fenicsproject.discourse.group/u/WjjSteve)\
**Post date:** [January 27, 2024, 12:58pm UTC](https://fenicsproject.discourse.group/t/how-to-apply-edge-size-and-normal-gradient/13420/6 "2024-01-27T12:58:59Z")

</div>

Hello, I just checked and I think the jump operator works well, thanks.

But the **FacetArea** cause some error:

```auto
UserWarning: Only know how to compute the facet area of an affine cell.
  warnings.warn("Only know how to compute the facet area of an affine cell.")

Traceback (most recent call last):
  File "/home/junjie/anaconda3/envs/fenicsx-env/bin/ffcx", line 10, in <module>
    sys.exit(main())
  File "/home/junjie/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/main.py", line 67, in main
    code_h, code_c = compiler.compile_ufl_objects(
  File "/home/junjie/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/compiler.py", line 107, in compile_ufl_objects
    code = generate_code(ir, options)
  File "/home/junjie/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/codegeneration.py", line 50, in generate_code
    code_integrals = [integral_generator(integral_ir, options) for integral_ir in ir.integrals]
  File "/home/junjie/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/codegeneration.py", line 50, in <listcomp>
    code_integrals = [integral_generator(integral_ir, options) for integral_ir in ir.integrals]
  File "/home/junjie/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/integrals.py", line 43, in generator
    parts = ig.generate()
  File "/home/junjie/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/integrals.py", line 211, in generate
    all_preparts += self.generate_piecewise_partition(rule)
  File "/home/junjie/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/integrals.py", line 353, in generate_piecewise_partition
    pre_definitions, parts = self.generate_partition(arraysymbol, F, "piecewise", None)
  File "/home/junjie/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/integrals.py", line 399, in generate_partition
    vaccess = self.backend.access.get(mt.terminal, mt, tabledata, quadrature_rule)
  File "/home/junjie/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/access.py", line 64, in get
    raise RuntimeError(f"Not handled: {type(e)}")
RuntimeError: Not handled: <class 'ufl.geometry.FacetArea'>

```

I write my py code like this:

```auto
import basix
from basix.ufl_wrapper import create_element
from ufl import Coefficient, TrialFunction, TestFunction, Constant, FacetArea, FacetNormal, jump, dx, dS, quadrilateral, inner, grad

cell = "quadrilateral"
element = create_element("Lagrange", cell, 1, basix.LagrangeVariant.equispaced)

w = TestFunction(element)

f0 = Coefficient(element)
vloc = Constant(cell)
h = FacetArea(cell)
n = FacetNormal(cell)

L = vloc * h * h * inner(jump(grad(f0), n),jump(grad(w), n)) * dS

```

is there anything wrong inside it? Thanks.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 28, 2024, 10:24am UTC](https://fenicsproject.discourse.group/t/how-to-apply-edge-size-and-normal-gradient/13420/7 "2024-01-28T10:24:07Z")

</div>

The issue is that quadrilateral meshes can be non-affine, which means that the area cannot be computed directly as:  
`ufl.geometry.FacetJacobianDeterminant(mesh)*ufl.geometry.ReferenceFacetVolume(mesh)`.  
as this would vary per quadrature point (this would mean that to get the exact area of a facet, one would have to first compute it as an integral over the facet).

Nevertheless, as your grid is a structured unit square, the FacetJacobianDeterminant will be constant, and thus you can use this product directly (see below for an example).

But do note that if your cell is non-affine (i.e. the quadrilaterals are not parallelograms, this would give an unexact result, as the Jacobian will vary over the facet.

```python
from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, 2, 3, cell_type=dolfinx.mesh.CellType.quadrilateral)

area = ufl.geometry.FacetJacobianDeterminant(mesh)
ref_volume = ufl.geometry.ReferenceFacetVolume(mesh)

print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(area*ref_volume*ufl.ds)))

```

---

<div class="post-metadata">

**Author:** ![simons](https://avatars.discourse-cdn.com/v4/letter/s/ea5d25/32.png) [@simons](https://fenicsproject.discourse.group/u/simons)\
**Post date:** [February 15, 2024, 5:36pm UTC](https://fenicsproject.discourse.group/t/how-to-apply-edge-size-and-normal-gradient/13420/8 "2024-02-15T17:36:48Z")

</div>

I encounter the same situation/error message, even though I am using a triangular mesh (generated with gmsh). Aren’t all triangles affine transformations of the base triangle? Or am I missing something obvious ?  
Below is a MWE of my situation.

```auto
import gmsh
from mpi4py import MPI
from dolfinx.fem import ( FunctionSpace, form )
from dolfinx.io import gmshio
from ufl import ( FacetNormal, FiniteElement, FacetArea, TestFunction, TrialFunction, dS, grad, jump )

gmsh.initialize()
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0

rectangle1 = gmsh.model.occ.addRectangle(0, 0, 0, 1.5, 0.41, tag=1)
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=gdim)
domain_marker = 1
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], domain_marker)
gmsh.model.setPhysicalName(volumes[0][0], domain_marker, "Domain")

gmsh.option.setNumber("Mesh.MeshSizeFactor", .2)
gmsh.option.setNumber("Mesh.Algorithm", 5) # Delaunay triangulation
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)

mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)

s_cg1 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
V = FunctionSpace(mesh, s_cg1)

u = TrialFunction(V)
v = TestFunction(V)

ha = FacetArea(mesh)
n = FacetNormal(mesh)

print(mesh.basix_cell())

bil_form = ha("+") * ha("-") * jump(grad(u), n) * jump(grad(v) , n) * dS
bil_form = form(bil_form)

```

And the error message is

```auto
Traceback (most recent call last):
  File "/home/simon/fail.py", line 38, in <module>
    bil_form = form(bil_form)
               ^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 188, in form
    return _create_form(form)
           ^^^^^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 183, in _create_form
    return _form(form)
           ^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 141, in _form
    ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
                              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/dolfinx/jit.py", line 56, in mpi_jit
    return local_jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/dolfinx/jit.py", line 204, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 199, in compile_forms
    raise e
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 190, in compile_forms
    impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 260, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/ffcx/compiler.py", line 107, in compile_ufl_objects
    code = generate_code(ir, options)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/ffcx/codegeneration/codegeneration.py", line 54, in generate_code
    code_integrals = [integral_generator(integral_ir, options) for integral_ir in ir.integrals]
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/ffcx/codegeneration/C/integrals.py", line 36, in generator
    parts = ig.generate()
            ^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/ffcx/codegeneration/integral_generator.py", line 132, in generate
    all_preparts += self.generate_piecewise_partition(rule)
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/ffcx/codegeneration/integral_generator.py", line 261, in generate_piecewise_partition
    pre_definitions, parts = self.generate_partition(arraysymbol, F, "piecewise", None)
                             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/ffcx/codegeneration/integral_generator.py", line 305, in generate_partition
    vaccess = self.backend.access.get(mt.terminal, mt, tabledata, quadrature_rule)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/simon/miniconda3/envs/fenicsx/lib/python3.12/site-packages/ffcx/codegeneration/access.py", line 64, in get
    raise RuntimeError(f"Not handled: {type(e)}")
RuntimeError: Not handled: <class 'ufl.geometry.FacetArea'>

```

Thank you in advance!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [March 20, 2024, 8:03pm UTC](https://fenicsproject.discourse.group/t/how-to-apply-edge-size-and-normal-gradient/13420/9 "2024-03-20T20:03:08Z")

</div>

> [@simons](#):
>
> Aren’t all triangles affine transformations of the base triangle?

You are setting the order to 2

> [@simons](#):
>
> ```auto
> gmsh.model.mesh.generate(gdim)
> gmsh.model.mesh.setOrder(2)
> 
> ```

A second order triangle has curved edges, and is non-affine.