How to apply edge size and normal gradient

General
1

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
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!

2

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.

1 Like
3

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?

4

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

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

Similarly you can use the ufl.jump 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 which may suit your needs.

As an example you could examine the biharmonic demo 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.

5

That is exactly what I need, thanks very much!

6

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

But the FacetArea cause some error:

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:

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.

7

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.

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)))
2 Likes
8

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.

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

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!

9

You are setting the order to 2

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