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
#
# SPDX-License-Identifier: MIT
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
x_dofs = self.x_dofs.links(i)
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
