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

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
9

What version of DOLFINx are you running? I cannot reproduce the error on my system. This error reminds me of before the support of codim-1 and codim-2 support, ref:

10

Sorry, I did something wrong in my script. I confirm that the script in the above post works.