ufl.FacetArea for quadrilaterals

General
1

Hello,

Below is a minimum working example showing that ufl.FacetArea can not be used with quadrilateral mesh.

from mpi4py import MPI
import numpy as np
import dolfinx
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.triangle) # OK
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)
V = dolfinx.fem.functionspace(msh, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

fa = ufl.FacetArea(msh)

a = ufl.dot(ufl.grad(u) * fa, ufl.grad(v)) * ufl.ds
assemble_matrix(dolfinx.fem.form(a), bcs=[])

I know that there is a workaround for ufl.CellVolume (see here). Is there a similar workaround for ufl.FacetArea?

Best,

Lucas

2

You can create a sub-mesh of the facets of the mesh, and use the h function on the submesh: dolfinx/python/dolfinx/mesh.py at b694be6e66f9043de039ebc18c166733d044409f · FEniCS/dolfinx · GitHub

3

Thank you for the answer, I understand the idea but I don’t know how to implement it. I tried:

facets = dolfinx.mesh.locate_entities(msh, msh.geometry.dim - 1, lambda x: np.full_like(x[0], 1, dtype=np.int8))
submesh_facets, _, _, _ = dolfinx.mesh.create_submesh(msh, msh.geometry.dim-1, facets)
Ve = dolfinx.fem.functionspace(submesh_facets, ("DG", 0))
facet_area = dolfinx.fem.Function(Ve)
facet_area.x.array[:] = submesh_facets.h(self._mesh.geometry.dim - 1, facets)
facet_area.x.scatter_forward()

But I get the error:

ValueError: Multiple domains found, making the choice of integration domain ambiguous.

So I guess that Ve should rather be defined on msh than submesh_facets but then I do not understand what to do.

4

Something like the following may help. Here I’ve used both the Mesh::h function as suggested by @dokken and an assembly of \int_\kappa \mathrm{d}\vec{x}, \kappa \in \mathcal{T}^h(\Omega), exploiting a degree 0 DG space. Both in this instance yield the same result; however, if you have a high order mesh, the latter will be correct (provided appropriate quadrature is employed).

from mpi4py import MPI
import numpy as np
import dolfinx
import dolfinx.fem
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
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)

facets = dolfinx.mesh.locate_entities(msh, msh.geometry.dim - 1, lambda x: np.full_like(x[0], 1, dtype=np.int8))
submesh_facets, _, _, _ = dolfinx.mesh.create_submesh(msh, msh.geometry.dim-1, facets)

Ve = dolfinx.fem.functionspace(submesh_facets, ("DG", 0))
child_facets = np.arange(submesh_facets.topology.index_map(1).size_local, dtype=np.int32)

facet_area_h = dolfinx.fem.Function(Ve)
facet_area_h.x.array[:] = submesh_facets.h(submesh_facets.geometry.dim - 1, child_facets)
facet_area_h.x.scatter_forward()

facet_area_q = dolfinx.fem.Function(Ve)
v = ufl.TestFunction(Ve)
# Make sure you use appropriate number of quadrature points
dx_e = ufl.Measure("dx", domain=submesh_facets)
facet_area_q.x.array[:] = assemble_vector(dolfinx.fem.form(v * dx_e))
facet_area_q.x.scatter_forward()

which yields in both cases:

image
image450×311 3.23 KB

5

Thank you Nate for your help. However I still have the error

ValueError: Multiple domains found, making the choice of integration domain ambiguous.

The issue is not really how to compute the facet areas but rather how to use in a form like ufl.dot(ufl.grad(u) * fa, ufl.grad(v)) * ufl.ds

Here is the full script:

from mpi4py import MPI
import numpy as np
import dolfinx
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.triangle)
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)
V = dolfinx.fem.functionspace(msh, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# fa = ufl.FacetArea(msh)
facets = dolfinx.mesh.locate_entities(msh, msh.geometry.dim - 1, lambda x: np.full_like(x[0], 1, dtype=np.int8))
submesh_facets, _, _, _ = dolfinx.mesh.create_submesh(msh, msh.geometry.dim-1, facets)

Ve = dolfinx.fem.functionspace(submesh_facets, ("DG", 0))
child_facets = np.arange(submesh_facets.topology.index_map(1).size_local, dtype=np.int32)

fa = dolfinx.fem.Function(Ve)
fa.x.array[:] = submesh_facets.h(submesh_facets.geometry.dim - 1, child_facets)
fa.x.scatter_forward()

a = ufl.dot(ufl.grad(u) * fa, ufl.grad(v)) * ufl.ds
assemble_matrix(dolfinx.fem.form(a), bcs=[])

As I explained in my previous post, I think that Ve should rather be defined on msh than submesh_facets.

6

Hello,

I still did not find any solution to my problem (basically computing something like ufl.dot(ufl.grad(u) * facet_area, ufl.grad(v)) * ufl.dS'). In the examples above, facet_area is defined as a function on submesh_facets then assembling the above expression fails because u and v are defined on msh, not submesh_facets. Should not facet_area be defined on msh?

7

I think you should use entity_maps, see for instance