Good morning,
I am trying to solve the Stokes problem with full Dirichlet boundary condition. In this case the problem is ill posed, because the pressure is defined up to a constant. To resolve the issue I am trying to add the null space. In the code below, what I am trying to do is to attach to the matrix a vector composed of zeros for the dofs relative to the velocity and ones for the dofs relative to the pressure. The problem is that somehow the compiler says that this is not the right basis.
Am I missing something?
import numpy as np
import dolfinx
import ufl
from petsc4py import PETSc
from mpi4py import MPI
N = 2
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N,
cell_type=dolfinx.mesh.CellType.hexahedron)
P2 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = ufl.MixedElement([P2, P1])
W = dolfinx.fem.FunctionSpace(mesh, TH)
def u_ex(x):
return np.vstack((4*(x[1]-1)*x[1]*(x[0]-1)*(x[0]-1)*x[0]*x[0]*(x[1]-0.5), -4*(x[0]-0.5)*x[1]*x[1]*x[0]*(x[1]-1)*(x[1]-1)*(x[0]-1), np.zeros_like(x[0])))
def p_ex(x):
return x[0]+x[1]+x[2]-3./2
V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()
### Define the boundary condition -> which in this case are only of Dirichlet type
uD = dolfinx.fem.Function(V)
uD.interpolate(u_ex)
# Create facet to cell connectivity required to determine boundary facets
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets)
bcs = [dolfinx.fem.dirichletbc(uD, boundary_dofs)]
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
x = ufl.SpatialCoordinate(mesh)
f = ufl.as_vector([1-24*(-1 + x[0])*(x[0]-1)*x[0]*x[0]*(-0.5 + x[1])- 48*(0.166667 - x[0] + x[0]*x[0])*(-1 + x[1])*(-0.5 + x[1])*x[1],
1+24*(-0.5 + x[0])*(-1 + x[1])*(x[1]-1)*x[1]*x[1] + 8*(-1 + x[0])*(-0.5 + x[0])*x[0]*(1 - 6*x[1] + 6*x[1]*x[1]),
1.0+0*x[0]])
dx = ufl.Measure("dx", domain=mesh)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx + \
ufl.div(v)*p*dx + q*ufl.div(u)*dx
L = ufl.inner(f, v) * dx
a_f = dolfinx.fem.form(a)
L_f = dolfinx.fem.form(L)
A = dolfinx.fem.petsc.assemble_matrix(a_f, bcs=bcs)
A.assemble()
# Tp print the matrix
from scipy.sparse import csr_matrix
assert isinstance(A, PETSc.Mat)
ai, aj, av = A.getValuesCSR()
Asp = csr_matrix((av, aj, ai))
print(Asp)
b = dolfinx.fem.petsc.assemble_vector(L_f)
dolfinx.fem.petsc.apply_lifting(b, [a_f], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bcs)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
# Create vector that spans the null space
# To see dof: https://fenicsproject.discourse.group/t/how-to-obtain-dofs-of-functionspace-in-dolfinx/7090
num_dofs_global_u = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
num_dofs_global = W.dofmap.index_map.size_local * W.dofmap.index_map_bs
null_vec = np.zeros(num_dofs_global)
print(num_dofs_global, ", ", null_vec.size)
null_vec[num_dofs_global_u:] = 1.0
print(null_vec)
def array2petsc4py(g):
Xpt = PETSc.Vec().create()
print("*")
Xpt.setSizes(g.shape[0])
print("*")
# Allocate memory for the PETSc vector
Xpt.setUp()
print("*")
# Copy the values from the NumPy array to the PETSc vector
Xpt.array = g
Xpt.assemble()
return Xpt
null_vec_petsc = array2petsc4py(null_vec)
null_vec_petsc.normalize()
nullspace = PETSc.NullSpace().create(vectors=[null_vec_petsc])
assert nullspace.test(A)
A.setNearNullSpace(nullspace)