Function at quadrature points of facets

General
1

TLDR: Can I define a function at quadrature points of facets and use it for a boundary integral?

I want to integrate on the boundary of mesh A a function coming from mesh B. From my understanding I should not mix functions on different meshes in a ufl form. Hence I have to interpolate first from mesh B to mesh A.

Since it is for a boundary integral I would like to have a function space on mesh A that has DOFs at quadrature points of facets. However I have trouble finding the right function space. Here is an incomplete MWE:

from dolfinx import fem, mesh, io
from mpi4py import MPI
import numpy as np
import ufl
import basix.ufl

def mesh_rectangle(min_x, min_y, max_x, max_y, elsize, cell_type=mesh.CellType.triangle):
    nx = np.ceil((max_x-min_x)/elsize).astype(np.int32)
    ny = np.ceil((max_y-min_y)/elsize).astype(np.int32)
    return mesh.create_rectangle(MPI.COMM_WORLD,
                                 [np.array([min_x,min_y]),np.array([max_x,max_y])],
                                 [nx,ny],
                                 cell_type,)

class Problem:
    def __init__(self, domain):
        self.domain = domain
        self.dim    = domain.topology.dim
        self.V    = fem.functionspace(self.domain,("CG",1))
        self.uh   = fem.Function(self.V, name="uh")

el_size = 0.25
mesh_left  = mesh_rectangle(0, 0, 0.5, 1, elsize=el_size)
mesh_right = mesh_rectangle(0.5, 0, 1, 1, elsize=el_size, cell_type=mesh.CellType.quadrilateral)

p_left  = Problem(mesh_left)
p_right = Problem(mesh_right)

x = ufl.SpatialCoordinate(mesh_right)
p_right.uh.interpolate(fem.Expression(x[0]**2+x[1]**2,
                                     p_right.uh.function_space.element.interpolation_points()))
dg0_right     = fem.functionspace(mesh_right,("DG",0))
func2interpolate = fem.Function(dg0_right)
func2interpolate.interpolate(fem.Expression(ufl.grad(p_right.uh)[0],
                             func2interpolate.function_space.element.interpolation_points()))

rfacets = mesh.locate_entities_boundary(p_left.domain,p_left.dim-1,lambda x:np.isclose(x[0],0.5))
# I want a function that:
# 1. Has dofs at quadrature points of rfacets on mesh_left
# 2. I can put in a ufl form with other functions on mesh_left

#Qe = basix.ufl.quadrature_element(mesh_left.topology.entity_types[-2][0].name,
                                  #degree=1,)
#Qs = fem.functionspace(mesh_left,Qe)

The function I want to interpolate is DG0 on mesh B. I tried making a quadrature space from a facet element as you can see in the last lines but I get an error.

Non-matching UFL cell and mesh cell shapes.

I am working in a docker container with the latest versions of the Fenicsx suit.

2

A question:
Will the facets of the two meshes align?

A few things you could try:

  • Use the non-matching interpolation (with a tiny padding of the element size if you want to work with DG-0) and specify the subset of cells for those facets close to the boundary (using compute_incident_entities on mesh A.
    or
  • Find the midpoint of each facet on mesh B. Use determine point ownership to find a cell in mesh A that has this midpoint in a cell. Then extract value from cell connected to facet of mesh B and assign value to DG 0 space of mesh A.
1 Like
3

In this case they do align:

image
image520×522 10.1 KB

Thank you for your suggestions. It seems that I should choose the same space as the function I am trying to project. I understand I can’t have a function on quadrature points of facets which I can use for boundary integrals.

4

The problem with choosing DG-0 is that the degrees of freedom are located at the center of the cell, making it more challenging to do with the non-matching interpolation algorithm (you could of course use alternative 2, which works nicely with DG-0).

As far a I am aware there is no support for facet quadrature spaces yet. I guess the work of @jpdean will make this possible, as one would make a mesh-view of the boundary and a quadrature space on this.

1 Like