Constrain assembled matrix in dolfinx

dolfinx
1

Hi,
My query is related to imposing constraints on assembly matrix which works in legacy-dolfin but, not the same in dolfinx. The MWE code below works for legacy-dolfin.

from dolfin import *
from scipy.sparse import csr_matrix
import numpy as np
import ufl
mesh=UnitSquareMesh(5,5)
Ve = VectorElement("CG", mesh.ufl_cell(), 1,dim=3) 
Re = VectorElement("R", mesh.ufl_cell(), 0,dim=4)
W = FunctionSpace(mesh, MixedElement([Ve, Re]))
V = FunctionSpace(mesh, Ve)

v_,lamb_ = TestFunctions(W)
dv, dlamb = TrialFunctions(W)
dx = Measure('dx')(domain=mesh)
C=as_tensor([[ 7.56043956e+08,  2.26813187e+08,  0.00000000e+00, 9.31322575e-10,  9.31322575e-10,  0.00000000e+00],
              [ 2.26813187e+08,  7.56043956e+08,  0.00000000e+00, 9.31322575e-10,  9.31322575e-10,  0.00000000e+00],
             [ 0.00000000e+00,  0.00000000e+00,  2.64600000e+08, 0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
            [-1.39698386e-09, -1.39698386e-09,  0.00000000e+00, 2.52014652e+06,  7.56043956e+05,  0.00000000e+00],
            [-1.39698386e-09, -1.39698386e-09,  0.00000000e+00, 7.56043956e+05,  2.52014652e+06,  0.00000000e+00],
            [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, 0.00000000e+00,  0.00000000e+00,  8.82000000e+05]])
def eps(v):
    return as_vector([0,v[1].dx(0),v[2].dx(1),v[1].dx(1)+v[2].dx(0),v[0].dx(1),v[0].dx(0)])

F2=dot(dot(C,eps(dv)),eps(v_))*dx
c1=lamb_[0]*dv[0]+lamb_[1]*dv[1]+lamb_[2]*dv[2]
c2=dlamb[0]*v_[0]+dlamb[1]*v_[1]+dlamb[2]*v_[2]
c3=lamb_[3]*(dv[1].dx(1)-dv[2].dx(0))+dlamb[3]*(v_[1].dx(1)-v_[2].dx(0))   

a2=F2+(c1+c2+c3)*dx
A2=assemble(a2)
ai, aj, av= as_backend_type(A2).mat().getValuesCSR()
A_csr=csr_matrix((av, aj, ai))
A=A_csr.toarray()  

print(A.shape)
print('Matrix rank after using constraints',np.linalg.matrix_rank(A))

# Without constraints
A2=assemble(F2)
ai, aj, av= as_backend_type(A2).mat().getValuesCSR()
A_csr=csr_matrix((av, aj, ai))
A=A_csr.toarray()  
print('Matrix rank without constraints',np.linalg.matrix_rank(A))

However, for implementing the constraints in dolfinx, I tried with the null constraints to restrict the rigid body motions. But, that doesn’t make affect. The below MWE in dolfinx does the same but, there is no affect in rank of amtrix. This shows the constraints are not implemented in coefficient matrix A_l. Kindly help me to resolve this issue.

import dolfinx
import basix
from dolfinx.fem import form, petsc, Function, functionspace
from ufl import TrialFunction, TestFunction, inner, lhs, rhs, dx, dot, eq, grad
from dolfinx import fem
import petsc4py.PETSc
from mpi4py import MPI
from contextlib import ExitStack
import numpy as np
import ufl
mesh_l=dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 5,5, dolfinx.mesh.CellType.triangle)

V_l = dolfinx.fem.functionspace(mesh_l, basix.ufl.element(
    "CG", mesh_l.topology.cell_name(), 1, shape=(3, )))

w_l = Function(V_l)
dv = TrialFunction(V_l)
v_= TestFunction(V_l)
dx_l= ufl.Measure('dx')(domain=mesh_l)

C=ufl.as_tensor([[ 7.56043956e+08,  2.26813187e+08,  0.00000000e+00, 9.31322575e-10,  9.31322575e-10,  0.00000000e+00],
              [ 2.26813187e+08,  7.56043956e+08,  0.00000000e+00, 9.31322575e-10,  9.31322575e-10,  0.00000000e+00],
             [ 0.00000000e+00,  0.00000000e+00,  2.64600000e+08, 0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
            [-1.39698386e-09, -1.39698386e-09,  0.00000000e+00, 2.52014652e+06,  7.56043956e+05,  0.00000000e+00],
            [-1.39698386e-09, -1.39698386e-09,  0.00000000e+00, 7.56043956e+05,  2.52014652e+06,  0.00000000e+00],
            [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, 0.00000000e+00,  0.00000000e+00,  8.82000000e+05]])
def eps(v):
    return ufl.as_vector([0,v[1].dx(0),v[2].dx(1),v[1].dx(1)+v[2].dx(0),v[0].dx(1),v[0].dx(0)])

F2=dot(dot(C,eps(dv)),eps(v_))*dx

#### Without applying constraints    ###############
A_l=dolfinx.fem.petsc.assemble_matrix(form(F2))
A_l.assemble()
print(A_l[:,:].shape)
print('Matrix rank without constraints',np.linalg.matrix_rank(A_l[:,:]))

# Nullspace Constraint implement

index_map = V_l.dofmap.index_map
nullspace_basis = [dolfinx.la.create_petsc_vector(index_map, V_l.dofmap.index_map_bs) for i in range(4)]
with ExitStack() as stack:
    vec_local = [stack.enter_context(x.localForm()) for x in nullspace_basis]
    basis = [np.asarray(x) for x in vec_local]

dofs = [V_l.sub(i).dofmap.list for i in range(3)]
# Build translational null space basis
for i in range(3):
    basis[i][dofs[i]] = 1.0
# Build rotational null space basis
x = V_l.tabulate_dof_coordinates()
dofs_block = V_l.dofmap.list
x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
basis[3][dofs[1]] = -x2
basis[3][dofs[2]] = x1
dolfinx.la.orthonormalize(nullspace_basis)
nullspace_l = petsc4py.PETSc.NullSpace().create(nullspace_basis, comm=MPI.COMM_WORLD)

A_l.setNullSpace(nullspace_l)   
print('Matrix rank after using constraints',np.linalg.matrix_rank(A_l[:,:]))

I don’t know where I am wrong in dolfinx to constrain A_l matrix for rigid body motion. Any help is greatly appreciated!

2

See GitHub - jorgensd/dolfinx_mpc: Extension for dolfinx to handle multi-point constraints.

3

Lets consider the following:
you want to remove rigid body motions to ensure that your PETSc solver can solve the problem. this is different from removing the nullspace from the matrix, see: MatSetNullSpace — PETSc 3.21.2 documentation

4

@dokken and @francesco-ballarin

Thanks for your response. Can you provide syntax of MatSetnullspace for dolfinx for any mwe? I find the reference syntax related to C language and no other related mwe.

Is it different from ‘’‘A_l.setnullspace(nullspace_l)’‘’ command as setnullspace is not working as shown in my mwe for dolfinx. My aim to apply the constraints and obtain assembled matrix where rigid body motion is removed.