How to integrate a vector known only at quadrature points on a boundary?

General
1

Hi everyone,

I’m trying to implement a problem where the right-hand side involves a vector that is only known at integration (quadrature) points.

For example, in a Poisson problem on a square domain, I would have a term like:

\int_{\Gamma_r} f \cdot v \ ds

where v is the usual test function, \Gamma_r is the right side of the domain, and f is a force that I only know at the quadrature points on that boundary.

I’m not sure how to proceed. I tried using basix.ufl.quadrature_element to align the nodes with the integration points, but I couldn’t get it to work.

What is the recommended way in FEniCSx to assemble such a term when the coefficient is only available at quadrature points? Should I manually assemble the contribution facet by facet, or is there a way to use a QuadratureElement or another built-in feature for this?

Best Regards

2

Hi @SamiGiroud,
what you would like can be achieved with submeshes, as illustrated below:

from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
import basix.ufl
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

tdim = mesh.topology.dim
mesh.topology.create_connectivity(tdim-1,tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)

submesh, cell_map, _, _ = dolfinx.mesh.create_submesh(mesh, tdim-1, boundary_facets)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
v = ufl.TestFunction(V)
tag = 3
ft = dolfinx.mesh.meshtags(mesh, tdim-1, boundary_facets, np.full_like(boundary_facets, tag, dtype=np.int32))
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=tag)

q_el = basix.ufl.quadrature_element(submesh.basix_cell(), degree=2)

Q = dolfinx.fem.functionspace(submesh, q_el)
f = dolfinx.fem.Function(Q)
f.interpolate(lambda x: x[0])
L = ufl.inner(f, v)*ds

if hasattr(cell_map, "sub_topology_to_topology"):
    entity_maps = [cell_map]
else:
    num_entities_local = mesh.topology.index_map(tdim-1).size_local + mesh.topology.index_map(tdim-1).num_ghosts 
    parent_to_sub = np.full(num_entities_local, -1, dtype=np.int32)
    parent_to_sub[cell_map] = np.arange(len(cell_map), dtype=np.int32)
    entity_maps = {submesh: parent_to_sub}

L_compiled = dolfinx.fem.form(L, entity_maps=entity_maps)


print(Q.tabulate_dof_coordinates())

print(dolfinx.fem.assemble_vector(L_compiled).array)

This code runs on both 0.9 and nightly DOLFINx.
The key is that you should now populate f with the data that you have. Q.tabulate_dof_coordinates() gives you the location of each of the quadrature point on the boundary.
Note that you can take a sub-set of the boundary by restricting the submesh to fewer facets, I just used all exterior facets for simplicity.

3

Hi @dokken,
Thanks for your help, I’ll look into it.
In the meantime, I tried another method. Do you think it could work as well?

def Compute_Vector_known_qpoints(msh,N,degree,Lnl,b):

    #return the UFL expression of the vector Lnl only known on quadrature points
    #msh is a square of lenght b with N * N quadratic cell

    points_ref, weights_ref = basix.make_quadrature(
        cell = basix.CellType(1),
        degree = degree
    )
    
    x = ufl.SpatialCoordinate(msh)
    Nq_point = int(np.ceil((d+1)/2))  #number of quadrature point by element 
    Res = 0
    
    for e in range(N):

        Xmin = e/N
        Xmax = Xmin + b/N
        Jac = (Xmax - Xmin)/2

        points = Xmin + 2*Jac*points_ref
        
        for i in range(Nq_point):
            
            if i == 0:
                aa = Xmin
            else:
                aa = (points[i-1] + points[i])/2
                aa = np.asarray(aa).reshape(-1)[0]
            
            if i == Nq_point - 1:
                bb = Xmax
            else:
                bb = (points[i] + points[i+1])/2
                bb = np.asarray(bb).reshape(-1)[0]
            
            inside = ufl.And(
                ufl.gt(x[1],aa),
                ufl.lt(x[1],bb)
            )
            
            Res += ufl.conditional(
                inside,
                ufl.as_vector((Lnl[i + Nq_point*e,0],Lnl[i + Nq_point*e,1])),
                ufl.as_vector((0.,0.))
            )
               
    return Res
4

The method you are proposing is going to be insanely expensive, as you are making a nested conditional depending on the number of quadrature points, which in turn will generate a massive kernel.

5

Oh, I see. Even if the number of quadrature points doesn’t exceed 3 or 4?

6

You are looping over all the entities along the boundary, making it dependent on the mesh discretization:

Thus you have N * Nq_point conditions within the generated code.

1 Like