ufl.FacetArea for quadrilaterals

General
8

Hi,

I took me while but I understood how to use entity_maps. But now I can compute form_a corresponding to a = ufl.inner(ufl.grad(u) * fa, ufl.grad(v)) * ds in the example below.

However, I don’t know how to compute form_b (see again in the example below) corresponding to b = ufl.inner(ufl.grad(u) * fa, ufl.grad(v)) * ds + ufl.inner(u, v) * dx, the following code raises the error when creating form_b:

Point dim (2) does not match element dim (1).

import dolfinx
from mpi4py import MPI
import numpy as np
from dolfinx.fem.petsc import assemble_matrix
import ufl

msh = dolfinx.mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((0.0, 0.0), (2.0, 1.0)), n=(5, 1), cell_type=dolfinx.mesh.CellType.quadrilateral)
gdim = msh.geometry.dim
dx = ufl.Measure('dx', domain=msh)
ds = ufl.Measure('ds', domain=msh)

V = dolfinx.fem.functionspace(msh, ("DG", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

def compute_facet_area(mesh):
    gdim = mesh.geometry.dim
    bdry_facets = dolfinx.mesh.locate_entities(mesh, gdim-1, lambda x: np.full_like(x[0], 1, dtype=np.int8))
    submesh_facets, entity_map, _, _ = dolfinx.mesh.create_submesh(mesh, gdim-1, bdry_facets)
    dx_e = ufl.Measure("dx", domain=submesh_facets, metadata={"quadrature_degree": 20})

    DG0 = dolfinx.fem.functionspace(submesh_facets, ("DG", 0))
    fa = dolfinx.fem.Function(DG0)
    v = ufl.TestFunction(DG0)

    fa.x.array[:] = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(v * dx_e))
    fa.x.scatter_forward()
    return fa, entity_map
fa, entity_map = compute_facet_area(msh)

a = ufl.inner(ufl.grad(u) * fa, ufl.grad(v)) * ds
form_a = dolfinx.fem.form(a, entity_maps=[entity_map])

b = ufl.inner(ufl.grad(u) * fa, ufl.grad(v)) * ds + ufl.inner(u, v) * dx
# form_b = dolfinx.fem.form(b) # KO
form_b = dolfinx.fem.form(b, entity_maps=[entity_map]) # KO