# Element-wise post-processing

Hello everyone,

I’m working with a hybrid method that provides me a solution t_h for the trace of the pressure field on every edge of my mesh.

To recover the pressure field in the whole domain, an element-wise post-processing needs to be used. Basically, for every element K we want to find p_K \in P2(K) such that

a(p_K, v) = f(v), \quad \forall v \in P2(K)
p_K = t_h \textrm{ on } \partial K,

where P2(K) denotes the space of quadratic polynomials over K and a(\cdot, \cdot) is my bilinear form. The final pressure approximation p_h belongs to the discontinuous global space DG2 and is given by

p_h|_K = p_K, \quad \forall K.

My question is how can I implement this kind of post-processing in FEniCSx. I know from here that I can generate local assemble matrices using the function assemble_local.

Can I then use such matrices to solve the local problems? If so, how do I cast the solutions of each local problem into the global DG2 space?

Giovanni

I’ve made a local assembler here:

You can then invert each local matrix to Get the solution per cell, and cast it to the u.x.array by accessing the cell dofs in the dofmap.

Thanks for your reply Dokken. I think I got the general idea but I still have to figure out some details.

Is the code you provided also able to assemble the linear form f(v) locally? In my case, the right-hand side is given by

f(v) = \int \textrm{div}(u_h) v \textrm{ dx}

Where u_h is a vector function on an approximation space V1 that was previously computed.

Also, provided a finite element space V on a mesh M, is there a simple way to locate all the dofs of V on the edges of M? Perhaps using the locate_dofs_topological function? This is a bit unrelated to the previous questions but is important for my code too.

Just by a slight modification of

I’ve made a local assembler (and compared it to a global assembly)

# Copyright (C) 2023 Jørgen S. Dokken
#
# Local assembly in DOLFINx using numpy
#

import time

import cffi
import dolfinx.fem.petsc as petsc
import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import div, ds, dx, grad, inner

class LocalAssembler():
def __init__(self, form):
self.form = fem.form(form)
self.update_coefficients()
self.update_constants()
subdomain_ids = self.form.integral_ids(fem.IntegralType.cell)
assert(len(subdomain_ids) == 1)
assert(subdomain_ids[0] == -1)
is_complex = np.issubdtype(ScalarType, np.complexfloating)
nptype = "complex128" if is_complex else "float64"
self.kernel = getattr(self.form.ufcx_form.integrals(fem.IntegralType.cell)[0], f"tabulate_tensor_{nptype}")
self.active_cells = self.form.domains(fem.IntegralType.cell, -1)
assert len(self.form.function_spaces) == self.form.rank
self.local_shape = [0 for _ in range(self.form.rank)]
for i, V in enumerate(self.form.function_spaces):
self.local_shape[i] = V.dofmap.dof_layout.block_size * V.dofmap.dof_layout.num_dofs

e0 = self.form.function_spaces[0].element
if self.form.rank == 1:
e1 = self.form.function_spaces[0].element
elif self.form.rank == 2:
e1 = self.form.function_spaces[1].element
else:
raise RuntimeError(f"Assembly for tensor of tensor of {rank=} is not supported")

needs_transformation_data = e0.needs_dof_transformations or e1.needs_dof_transformations or \
self.form.needs_facet_permutations
if needs_transformation_data:
raise NotImplementedError("Dof transformations not implemented")

self.ffi = cffi.FFI()
V = self.form.function_spaces[0]
self.x_dofs = V.mesh.geometry.dofmap
self.local_tensor = np.zeros(self.local_shape, dtype=ScalarType)

def update_coefficients(self):
self.coeffs = fem.assemble.pack_coefficients(self.form)[(fem.IntegralType.cell, -1)]

def update_constants(self):
self.consts = fem.assemble.pack_constants(self.form)

def update(self):
self.update_coefficients()
self.update_constants()

def assemble(self, i:int):
x = self.form.function_spaces[0].mesh.geometry.x
geometry = np.zeros((len(x_dofs), 3), dtype=np.float64)
geometry[:, :] = x[x_dofs]

facet_index = np.zeros(0, dtype=np.intc)
facet_perm = np.zeros(0, dtype=np.uint8)
if self.coeffs.shape == (0, 0):
coeffs = np.zeros(0, dtype=ScalarType)
else:
coeffs = self.coeffs[i,:]
ffi_fb = self.ffi.from_buffer
self.local_tensor.fill(0)
self.kernel(ffi_fb(self.local_tensor), ffi_fb(coeffs), ffi_fb(self.consts), ffi_fb(geometry),
ffi_fb(facet_index), ffi_fb(facet_perm))
return self.local_tensor

N = 20
msh = mesh.create_box(MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]),
np.array([2.0, 1.0, 1.0])], [N, N, N],
mesh.CellType.tetrahedron)

Q = fem.FunctionSpace(msh, ("DG", 2))

p = ufl.TrialFunction(Q)
q = ufl.TestFunction(Q)
a = inner(p, q) *dx
lhs_assembler = LocalAssembler(a)

W = fem.VectorFunctionSpace(msh,  ("Lagrange", 1))
uh = fem.Function(W)
uh.interpolate(lambda x: (x[0], 2*x[1], x[2]**2))

l = inner(div(uh), q)*dx
rhs_assembler = LocalAssembler(l)

w1 = fem.Function(Q)
bs = Q.dofmap.bs
start = time.perf_counter()
for cell in range(msh.topology.index_map(msh.topology.dim).size_local):
A = lhs_assembler.assemble(cell)
b = rhs_assembler.assemble(cell)
cell_dofs = Q.dofmap.cell_dofs(cell)
unrolled_dofs = np.zeros(len(cell_dofs)*bs, dtype=np.int32)
for i, dof in enumerate(cell_dofs):
for j in range(bs):
unrolled_dofs[i*bs+j] = dof*bs+j
w1.x.array[unrolled_dofs] = np.linalg.solve(A, b)
w1.x.scatter_forward()
end = time.perf_counter()
print(f"Local assembly {end-start:.2e}")

start = time.perf_counter()
w2 = fem.Function(Q)
problem = petsc.LinearProblem(a, l, u=w2)
problem.solve()
w2.x.scatter_forward()
end = time.perf_counter()
print(f"Global assembly {end-start:.2e}")

from dolfinx.io import VTXWriter

with VTXWriter(msh.comm, "w1.bp", [w1]) as vtx:
vtx.write(0)

with VTXWriter(msh.comm, "w2.bp", [w2]) as vtx:
vtx.write(0)

assert np.allclose(w1.x.array, w2.x.array)


Hello Dokken,

Thanks to the code you provided, I managed to implement the original local post-processing I was interested in.

Now I want to study a different local post-processing defined as: For each element K, find (z_h, p_h) \in RT_1(K) \times P_1(K) such that

\int_K z_h \cdot v_h \textrm{ dx} \,- \int_K p_h \, \textrm{div} v_h \textrm{ dx} = - \int_{\partial K} s_e \, (v_h \cdot n) \textrm{ ds} \quad \forall v_h \in RT_1(K)
\int_K \textrm{div} z_h \, q_h \textrm{dx} = \int_K u_e \, q_h \textrm{ dx} \quad \forall q_h \in P_1(K),

where RT_1(K) is the Raviart-Thomas vector space of index 1 on K, P_1(K) is the space of linear polynomials over K, s_e and u_e are previously computed scalar functions, and n is the unitary normal vector.

Here is my attempt to modify your code to implement such mixed post-processing (I did not modify the LocalAssembler class, only the code right after).

# Copyright (C) 2023 Jørgen S. Dokken
#
# Local assembly in DOLFINx using numpy
#

import time

import cffi
import dolfinx.fem.petsc as petsc
import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import div, ds, dx, grad, inner, dot
import os

class LocalAssembler():
def __init__(self, form):
self.form = fem.form(form)
self.update_coefficients()
self.update_constants()
subdomain_ids = self.form.integral_ids(fem.IntegralType.cell)
assert(len(subdomain_ids) == 1)
assert(subdomain_ids[0] == -1)
is_complex = np.issubdtype(ScalarType, np.complexfloating)
nptype = "complex128" if is_complex else "float64"
self.kernel = getattr(self.form.ufcx_form.integrals(fem.IntegralType.cell)[0], f"tabulate_tensor_{nptype}")
self.active_cells = self.form.domains(fem.IntegralType.cell, -1)
assert len(self.form.function_spaces) == self.form.rank
self.local_shape = [0 for _ in range(self.form.rank)]
for i, V in enumerate(self.form.function_spaces):
self.local_shape[i] = V.dofmap.dof_layout.block_size * V.dofmap.dof_layout.num_dofs

e0 = self.form.function_spaces[0].element
if self.form.rank == 1:
e1 = self.form.function_spaces[0].element
elif self.form.rank == 2:
e1 = self.form.function_spaces[1].element
else:
raise RuntimeError(f"Assembly for tensor of tensor of {rank=} is not supported")

needs_transformation_data = e0.needs_dof_transformations or e1.needs_dof_transformations or \
self.form.needs_facet_permutations
if needs_transformation_data:
raise NotImplementedError("Dof transformations not implemented")

self.ffi = cffi.FFI()
V = self.form.function_spaces[0]
self.x_dofs = V.mesh.geometry.dofmap
self.local_tensor = np.zeros(self.local_shape, dtype=ScalarType)

def update_coefficients(self):
self.coeffs = fem.assemble.pack_coefficients(self.form)[(fem.IntegralType.cell, -1)]

def update_constants(self):
self.consts = fem.assemble.pack_constants(self.form)

def update(self):
self.update_coefficients()
self.update_constants()

def assemble(self, i:int):
x = self.form.function_spaces[0].mesh.geometry.x
geometry = np.zeros((len(x_dofs), 3), dtype=np.float64)
geometry[:, :] = x[x_dofs]

facet_index = np.zeros(0, dtype=np.intc)
facet_perm = np.zeros(0, dtype=np.uint8)
if self.coeffs.shape == (0, 0):
coeffs = np.zeros(0, dtype=ScalarType)
else:
coeffs = self.coeffs[i,:]
ffi_fb = self.ffi.from_buffer
self.local_tensor.fill(0)
self.kernel(ffi_fb(self.local_tensor), ffi_fb(coeffs), ffi_fb(self.consts), ffi_fb(geometry),
ffi_fb(facet_index), ffi_fb(facet_perm))
return self.local_tensor

# No modifications until here -------------------------------------------------

N = 30
msh = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0.0, 0.0]),
np.array([0.5, 0.5])], [N, N],
mesh.CellType.triangle)

x = ufl.SpatialCoordinate(msh) # Spatial coordinates
nn = ufl.FacetNormal(msh)      # Unitary normal vector

# Source terms
se_ = lambda x: ufl.sin(ufl.pi*x[0])*ufl.sin(ufl.pi*x[1])
ue_ = lambda x: 2*ufl.pi*ufl.pi*ufl.sin(ufl.pi*x[0])*ufl.sin(ufl.pi*x[1])

se = se_(x)
ue = ue_(x)

# Mixed approximation space
el_vd = ufl.FiniteElement('DRT', msh.ufl_cell(), 2)
el_p = ufl.FiniteElement('DG', msh.ufl_cell(), 1)
W = fem.FunctionSpace(msh, ufl.MixedElement([el_vd, el_p]))

# Test and trial functions
TrialF = ufl.TrialFunction(W)
TestF  = ufl.TestFunction(W)
(z,p) = ufl.split(TrialF)
(v,q)  = ufl.split(TestF)

a = inner(z, v)*dx
a += -inner(p, div(v))*dx
a += inner(div(z), q)*dx

L = -se*dot(v, nn)*ds
L += inner(ue, q)*dx

lhs_assembler = LocalAssembler(a)
rhs_assembler = LocalAssembler(L)

wh = fem.Function(W)
bs = W.dofmap.bs
for cell in range(msh.topology.index_map(msh.topology.dim).size_local):
A = lhs_assembler.assemble(cell)
b = rhs_assembler.assemble(cell)
cell_dofs = W.dofmap.cell_dofs(cell)
unrolled_dofs = np.zeros(len(cell_dofs)*bs, dtype=np.int32)
for i, dof in enumerate(cell_dofs):
for j in range(bs):
unrolled_dofs[i*bs+j] = dof*bs+j
wh.x.array[unrolled_dofs] = np.linalg.solve(A, b)
wh.x.scatter_forward()

(vdh, ph) = wh.split()


However, the solutions I got seem wrong. In particular, the y component of z_h has a noticeable problem in the top boundary, as shown in the figure below.

My guess is that the problem is in the term - \int_{\partial K} s_e \, (v_h \cdot n) \textrm{ ds}. For the post-processing to be correct, the normal vectors n must be outwards (relative to the element K), and I think that is not always satisfied when using ufl.FacetNormal.

Any ideas on how to fix it?

Sorry to take more of your time,
Giovanni

The code is missing the post-processing. How are you making that visualization?

Facetnormals always point out of a cell.

Sorry, I forgot to add that to the post. But I’m using

# Plot ---

with io.XDMFFile(msh.comm, 'LM-post.xdmf', 'w') as xdmf_file:
xdmf_file.write_mesh(msh)
xdmf_file.write_function(vdh)


at the end of the computation. Then I open the .xmdf file in Paraview to viausalize the y component.

Obs: In the code I’m calling the variable z_h by vdh and the variable p_h by ph.

If n is always outwards, is there something else that could go wrong with the boundary integrals \int_{\partial K} p_e \, (v_h \cdot n) \textrm{ ds}?

If you use XDMFFile, the fields are interpolated into first order Lagrange, thus you lose alot of information.

Interpolate the function into a suitable DG space and use the VTXWriter, as it can write higher ord er Lagrange/DG

I guess it would be something like that?

# Plot ---

V = fem.VectorFunctionSpace(msh,("DG",2))
zh = fem.Function(V)
zh.interpolate(vdh)

from dolfinx.io import VTXWriter

with VTXWriter(msh.comm, "z.bp", [zh]) as vtx:
vtx.write(0.0)


but I’m having troubles trying to open the resulting files on paraview (version 5.11.1).

Yes, something along those lines.
What version of adios do you have installed?
There are some issues with the latest adios as it defaults to bp5, while paraview only supports bp4: Add ADIOS2 engine as option by garth-wells · Pull Request #2636 · FEniCS/dolfinx · GitHub
Paraview opens ADIOS2-VTX BP5 files as BP4 (#22108) · Issues · ParaView / ParaView · GitLab

I’m still working on solving the adios problem. A colleague from my lab will help me with that later.

Meanwhile, I changed the code to print the source vector ‘b’ for each cell.

wh = fem.Function(W)
bs = W.dofmap.bs
for cell in range(msh.topology.index_map(msh.topology.dim).size_local):
A = lhs_assembler.assemble(cell)
b = rhs_assembler.assemble(cell)
print('Source vector cell', cell)
print(b)
print('----------------------------------------')
cell_dofs = W.dofmap.cell_dofs(cell)
unrolled_dofs = np.zeros(len(cell_dofs)*bs, dtype=np.int32)
for i, dof in enumerate(cell_dofs):
for j in range(bs):
unrolled_dofs[i*bs+j] = dof*bs+j
wh.x.array[unrolled_dofs] = np.linalg.solve(A, b)
wh.x.scatter_forward()


For all cells, the first six components, which correspond to the integral \int_{\partial K} s_e (v_h \cdot n) \textrm{ ds}, are zero. This should not be the case for a general function s_e, so I’m definitely doing something wrong in the implementation of that term.

Any ideas of what could it be? Maybe something to do with the integration measure ‘ds’?

For me to have time to look at this, you need to simplify your example.

Start by doing a projection of functions using one of the appropriate spaces (similar to what I did in my post). Then one can clarify if it is due to a single function space, the mixed space or a combination of the two.

This is possible but non-trivial. The most detailed example I can give you is here:

https://www.sciencedirect.com/science/article/pii/S0898122122004722

I need to update this code for v0.7.0 of DOLFINx.

It includes all of the steps you mention (local assembly including facet terms, local application of boundary conditions, local solve, assembly back into another function space), but for a different problem. It could probably be adapted to your problem.

1 Like

Thanks for your reply, Jack! I ended up using other tools to implement this local post-processing. I will look at the code you provided to try to implement it in FEniCS again.