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.
2 Likes
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
5

I tried to follow your second suggestion and got this solution:

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

comm = MPI.COMM_WORLD
rank = comm.rank

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

def interpolate_dg0_at_facets(sending_f,
                             receiving_f,
                             facets,bb_tree_ext,):
    domain           = receiving_f.function_space.mesh
    ext_domain       = sending_f.function_space.mesh
    cdim = domain.topology.dim
    function_dim = 1 if (len(receiving_f.ufl_shape) == 0) else receiving_f.ufl_shape[0]
    # Build Gamma midpoints array
    local_interface_midpoints = np.zeros((len(facets), 3), np.double)
    for i, ifacet in enumerate(facets):
        local_interface_midpoints[i,:] = mesh.compute_midpoints(domain,cdim-1,np.array([ifacet],dtype=np.int32))

    facet_counts  = np.zeros(comm.size, dtype=np.int32)
    facets_offsets = np.zeros(comm.size, dtype=np.int32)
    comm.Allgather(np.array([len(facets)], np.int32), facet_counts)
    facets_offsets[1:] = np.cumsum(facet_counts[:-1])
    total_facet_count = np.sum(facet_counts, dtype=int)

    global_interface_midpoints = np.zeros((total_facet_count,3), dtype=np.double, order='C')
    comm.Allgatherv(local_interface_midpoints,[global_interface_midpoints,facet_counts*local_interface_midpoints.shape[1],facets_offsets*local_interface_midpoints.shape[1],MPI.DOUBLE])

    # Collect values at midpoints
    local_vals  = np.zeros((total_facet_count,function_dim),dtype=np.double,order='C')
    global_vals = np.zeros((total_facet_count,function_dim),dtype=np.double,order='C')
    found_local  = np.zeros(total_facet_count,dtype=np.double,order='C')
    found_global = np.zeros(total_facet_count,dtype=np.double,order='C')
    for idx in range(total_facet_count):
        candidate_parents_ext = geometry.compute_collisions_points(bb_tree_ext,global_interface_midpoints[idx,:])
        potential_parent_els_ext = geometry.compute_colliding_cells(ext_domain, candidate_parents_ext, global_interface_midpoints[idx,:])
        #potential_parent_els_ext = potential_parent_els_ext.array[np.flatnonzero( ext_activation_tag.values[ potential_parent_els_ext.array] ) ]
        if len(potential_parent_els_ext.array)>0:
            idx_owner_el = potential_parent_els_ext.array[0]
            if idx_owner_el < ext_domain.topology.index_map(cdim).size_local:
                local_vals[idx,:]  = sending_f.eval(global_interface_midpoints[idx,:], idx_owner_el)
                found_local[idx] = 1
    comm.Allreduce([local_vals, MPI.DOUBLE], [global_vals, MPI.DOUBLE])
    comm.Allreduce([found_local, MPI.DOUBLE], [found_global, MPI.DOUBLE])

    f_to_c_left = domain.topology.connectivity(1,2)

    # build global parent el array for facets
    global_parent_els_proc = np.zeros(len(facets), np.int32)
    for idx, ifacet in enumerate(facets):
        parent_els  = f_to_c_left.links(ifacet)
        #parent_els  = parent_els[np.flatnonzero(activation_tag.values[parent_els])]
        assert (len(parent_els)) >= 1
        parent_el_glob  = domain.topology.index_map(domain.geometry.dim).local_to_global(parent_els)
        global_parent_els_proc[idx] = parent_el_glob[0]

    global_parent_els = np.zeros(total_facet_count, np.int32)
    comm.Allgatherv(global_parent_els_proc,[global_parent_els,facet_counts,facets_offsets,MPI.INT])
    local_parent_els  = domain.topology.index_map(domain.geometry.dim).global_to_local(global_parent_els)

    for idx, el in enumerate(local_parent_els):
        if el < 0:
            continue
        flat_idx    = el*function_dim
        receiving_f.x.array[flat_idx:flat_idx+2] = global_vals[idx]

    receiving_f.vector.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)


class Problem:
    def __init__(self, domain, name="case"):
        self.domain = domain
        self.name = name
        self.V    = fem.functionspace(self.domain,("CG",1))
        self.dg0  = fem.functionspace(self.domain, ("DG", 0))
        dg0_dim2_element = basix.ufl.element("DG", self.domain.basix_cell(), 0, shape=(2,))
        self.dg0_dim2 = fem.functionspace(self.domain, dg0_dim2_element)
        self.uh   = fem.Function(self.V, name="uh")
        self.dim  = domain.topology.dim

        self.domain.topology.create_entities(self.dim-1)
        self.cell_map  = self.domain.topology.index_map(self.dim)
        self.facet_map = self.domain.topology.index_map(self.dim-1)

        self.bb_tree = geometry.bb_tree(domain,self.dim)
        self.domain.topology.create_connectivity(1,2)


    def writepos(self, funcs=[]):
        with io.VTKFile(self.domain.comm, f"out/{self.name}.pvd", "w") as ofile:
            ofile.write_function(funcs)

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

p_right = Problem(mesh_right, name="right")
p_left  = Problem(mesh_left, name="left")

d_f   = fem.Function(p_right.dg0,name="d_f")
der   = fem.Function(p_right.dg0_dim2,name="der")
c_f   = fem.Function(p_right.V,name="c_f")
x = ufl.SpatialCoordinate(mesh_right)
c_f.interpolate( fem.Expression(x[0]**2 + x[1]**2,c_f.function_space.element.interpolation_points()))
d_f.interpolate( fem.Expression(x[0],d_f.function_space.element.interpolation_points()))
der.interpolate( fem.Expression(ufl.grad(c_f),der.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))

der_python = fem.Function(p_left.dg0_dim2,name="der")
interpolate_dg0_at_facets(der,
                          der_python,
                          rfacets,
                          p_right.bb_tree)

p_left.writepos(funcs=[der_python])
p_right.writepos(funcs=[c_f,d_f,der])

image
image750×511 57.8 KB

I am satisfied with the interpolation and am now concerned with performance.

I call here sending mesh the mesh of the function to interpolate and receiving mesh the mesh of the interpolate. So what I am doing currently in the attached script is:

  1. Compute midpoints of local interface facets of receiving mesh.
  2. Global communication so that midpoints end up on all processors.
  3. Each processor checks if it owns the point on the sending mesh. If it does, evaluate the function to interpolate.
  4. Global communication so that values end up on all processors.
  5. Processors set the solution on the cells they own on the receiving end.

I would like to reduce the global communication. I think I should build a bounding box tree of the midpoints, create a global tree and go from there but I don’t quite understand what is a global tree. I am asking for help in sketching out what I should do.

6

I would suggest you look into how determine_point_ownership works, as you could use it to find a unique owner for each midpoint in the other mesh, which in turn can be used with MPI neighbourhood communicators to make to communication scalable.

1 Like