Efficiently evaluatuating surface integrals on each cell in mesh

General
1

Question

Is there a way to efficiently evaluate line/surface integrals on each cell?

Below is some example code which should demonstrate what I’m trying to do. My issue seems to be that repeated evaluations of either form() or assemble_scalar() is very slow.

For the case of volume/area integrals, I found that I could project the function of interest onto a DG 0 function space, and then use assemble_vector. I’m having trouble adapting that approach to the boundaries, as it requires some additional restriction on the discontinuity.

Code


from mpi4py import MPI     # Parallel communication
from dolfinx import mesh   # Mesh tools

from dolfinx import fem

import numpy as np
import ufl

from petsc4py.PETSc import ScalarType


#%% Domain

nx = 5
ny = 5
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.quadrilateral)


#%% Mesh Characterization
    
tdim = domain.topology.dim
domain.topology.create_connectivity(tdim, 1)     # Compute cells
c_to_e = domain.topology.connectivity(tdim, 1)   # Map cells to edges

num_cells = domain.topology.index_map(tdim).size_local + domain.topology.index_map(tdim).num_ghosts

# Tag each edge
num_edges = domain.topology.index_map(1).size_local + domain.topology.index_map(1).num_ghosts
tag_list = np.arange(num_edges, dtype=np.int32)
edge_tags = mesh.meshtags(domain, 1, tag_list, tag_list)

# Form Components
dE = ufl.Measure("dS", domain, subdomain_data=edge_tags)
cnst = fem.Constant(domain, ScalarType(1))


#%% Integration

L_calc = np.zeros(num_cells)
for itr in range(num_cells):
    
    L_val = 0
    for itrE in c_to_e.links(itr):
        L_val += fem.assemble_scalar( fem.form(cnst*dE(itr)) )

    L_calc[itr] = L_val
2

Could you please add the code where you tried to use a DG-0 test function?

3

For the boundary integrals, I am trying something like:

#%% DG-0 Projection


Vd = fem.FunctionSpace(domain, ("DG", 0))
v = ufl.TestFunction(Vd)

F_form = fem.form( ufl.inner(cnst, v)*dE )
F_asse = fem.assemble_vector(F_form)

Which gives the error:

UFLException: Discontinuous type Argument must be restricted.

My understanding is that I must specify something like the inner/outer edge of the discontinuity, so I try things like v(‘+’), or the following

F_form = fem.form( ufl.inner(cnst, ufl.avg(v))*dE )
F_asse = fem.assemble_vector(F_form)

F_vals = F_asse.array

The second attempt will produce results, but I expect it should be an array with all values equal to 0.8, the perimeter of the quadrilateral mesh elements.
I should note that I am just trying things to see what works at this point. I don’t understand the “restrictions” very well.

4

So I’ve made some progress on this, but I still don’t quite understand how restrictions in DG-0 work.
I’m able to get the answer I expect by using 2*avg(v), but I don’t know why the factor of 2 is necessary. My thinking was that the value on either side of the discontinuity should be the same, as I am projecting a constant c = 1 onto the space, so
\text{avg}(v) = \dfrac{1}{2}(v^+ \ + \ v^-) = v^+ = v^- \ ,
but this does not appear to be the case.

Below is a code snippet which seems to correctly evaluate the boundary integrals on each cell, excluding the exterior facets.

Code

from mpi4py import MPI     # Parallel communication
from dolfinx import mesh   # Mesh tools

from dolfinx import fem

import numpy as np
import ufl

from petsc4py.PETSc import ScalarType


#%% Domain

nx = 100
ny = nx
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.quadrilateral)


#%% DG-0 Projection

Vd = fem.FunctionSpace(domain, ("DG", 0))
v = ufl.TestFunction(Vd)

cnst = fem.Constant(domain, ScalarType(1))
F_form = fem.form( ufl.inner(cnst, 2*ufl.avg(v))*ufl.dS ) # *** Not sure why factor of 2 is necessary
F_asse = fem.assemble_vector(F_form)


#%% Checking

# Count the occurances of each integral value
F_vals = F_asse.array
FU_vals, FU_count = np.unique( np.round(F_vals, 4), return_counts=True)
# Formulas for integrals in a unit square with quadrilateral mesh
# Integrals will not be evaluated on boundary, due to use of dS as measure
#   Types of Elements: [corners, edges, interior]
#                      ie, [2 exterior, 1 exterior, 0 exterior]
ref_vals = np.round( np.array([2, 3, 4])/nx, 4)
ref_count = [4, 4*(nx - 2), nx**2 - 4*(nx-1)]

if (all(FU_vals == ref_vals)):
    print('Values Correct')

if (all(FU_count == ref_count)):
    print('Count Correct')