Hi all,

The following code is a very slightly modified code taken from the Stokes equations tutorial. I want to check how “well” the constrain \nabla \cdot \mathbf{u}=\partial _xu_x+\partial _yu_y=0 holds. I need it for a 3d example I am trying to solve…

How can I calculate this quantity in fenicsx? I tried it with no success I guess due to the type of the finite elements `ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)`

?

Any help will be greatly appreciated!!!

```
import numpy as np
import matplotlib.pyplot as plt
import ufl
import dolfinx
from dolfinx import cpp as _cpp
from dolfinx import fem
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
extract_function_spaces, form,
locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_rectangle,
locate_entities_boundary)
from ufl import div, dx, grad, inner
from mpi4py import MPI
from petsc4py import PETSc
import math
# Create mesh
msh = create_rectangle(MPI.COMM_WORLD,
[np.array([0, 0]), np.array([1, 1])],
[32, 32],
CellType.triangle, GhostMode.none)
# Generate experimental data
mesh=msh
# Create the function space
P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V, Q = FunctionSpace(msh, P2), FunctionSpace(msh, P1)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, boundary_facets)
dof_coordinates = V.tabulate_dof_coordinates()[boundary_dofs]
points1 = mesh.geometry.x
# Export the mesh points and the boundary points
length=len(points1)
length_BC=len(boundary_dofs)
x_points=np.zeros(length)
y_points=np.zeros(length)
z_points=np.zeros(length)
x_points_boundary=np.zeros(length_BC)
y_points_boundary=np.zeros(length_BC)
z_points_boundary=np.zeros(length_BC)
x_points = [points1[i][0] for i in range(length)]
y_points = [points1[i][1] for i in range(length)]
z_points = [points1[i][2] for i in range(length)]
x_points_boundary = [dof_coordinates[i][0] for i in range(length_BC)]
y_points_boundary = [dof_coordinates[i][1] for i in range(length_BC)]
z_points_boundary = [dof_coordinates[i][2] for i in range(length_BC)]
points_boundary=np.zeros((length_BC,3))
points_boundary=[(x_points_boundary[i],y_points_boundary[i],z_points_boundary[i]) for i in range(length_BC)]
points = np.zeros((2, len(points1)))
points=[[x_points[i],y_points[i],z_points[i]] for i in range(length)]
# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
return np.logical_or(np.logical_or(np.isclose(x[0], 0.0),
np.isclose(x[0], 1.0)),
np.isclose(x[1], 0.0))
# Function to mark the lid (y = 1)
def lid(x):
return np.isclose(x[1], 1.0)
# Lid velocity
def lid_velocity_expression(x):
return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))
# No-slip boundary condition for velocity field (`V`) on boundaries
# where x = 0, x = 1, and y = 0
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)
# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = Function(V)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
bc1 = dirichletbc(lid_velocity, locate_dofs_topological(V, 1, facets))
# Collect Dirichlet boundary conditions
bcs = [bc0, bc1]
# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))
a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx],
[inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx])
a_p11 = form(inner(p, q) * dx)
a_p = [[a[0][0], None],
[None, a_p11]]
########## SOLUTION ##########
A = fem.petsc.assemble_matrix_nest(a, bcs=bcs)
A.assemble()
P11 = fem.petsc.assemble_matrix(a_p11, [])
P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]])
P.assemble()
b = fem.petsc.assemble_vector_nest(L)
# Modify ('lift') the RHS for Dirichlet boundary conditions
fem.petsc.apply_lifting_nest(b, a, bcs=bcs)
# Sum contributions from ghost entries on the owner
for b_sub in b.getNestSubVecs():
b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
# Set Dirichlet boundary condition values in the RHS
bcs0 = fem.bcs_by_block(extract_function_spaces(L), bcs)
fem.petsc.set_bc_nest(b, bcs0)
# Create nullspace vector
null_vec = fem.petsc.create_vector_nest(L)
# Set velocity part to zero and the pressure part to a non-zero constant
null_vecs = null_vec.getNestSubVecs()
null_vecs[0].set(0.0), null_vecs[1].set(1.0)
# Normalize the vector, create a nullspace object, and attach it to the
# matrix
null_vec.normalize()
nsp = PETSc.NullSpace().create(vectors=[null_vec])
assert nsp.test(A)
A.setNullSpace(nsp)
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A, P)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
# Define the matrix blocks in the preconditioner with the velocity and
# pressure matrix index sets
nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(
("u", nested_IS[0][0]),
("p", nested_IS[0][1]))
# Set the preconditioners for each block
ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")
# Monitor the convergence of the KSP
ksp.setFromOptions()
u, p = Function(V), Function(Q)
x = PETSc.Vec().createNest([_cpp.la.petsc.create_vector_wrap(u.x), _cpp.la.petsc.create_vector_wrap(p.x)])
ksp.solve(b, x)
div_u = ufl.grad(u)
print(np.shape(points))
u_values = []
p_values = []
div_u_values = []
from dolfinx import geometry
bb_tree = geometry.BoundingBoxTree(msh, msh.topology.dim)
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions(bb_tree, points)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(msh, cell_candidates, points)
for i, point in enumerate(points):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_values = u.eval(points_on_proc, cells)
p_values = p.eval(points_on_proc, cells)
#div_u_values = div_u.eval(points_on_proc, cells)
uu=np.array(u_values)
p1=np.array(p_values)
ux=np.zeros(length)
uy=np.zeros(length)
pp=np.zeros(length)
ux=[uu[i][0] for i in range(length)]
uy=[uu[i][1] for i in range(length)]
pp=[p1[i] for i in range(length)]
norm_u_0 = u.x.norm()
norm_p_0 = p.x.norm()
if MPI.COMM_WORLD.rank == 0:
print("(A) Norm of velocity coefficient vector (nested, iterative): {}".format(norm_u_0))
print("(A) Norm of pressure coefficient vector (nested, iterative): {}".format(norm_p_0))
```