Thanks @francesco-ballarin for the explanation. I moved forward with the tutorial for my problem implementation.
I made update in above MWE with mesh as mesh = mesh.create_unit_square(MPI.COMM_WORLD, 1, 1, mesh.CellType.quadrilateral)
and function space as Ve = VectorElement("Lagrange", mesh.ufl_cell(), 1,dim=2)
for validating the theory. This gives 4 noded quad element mesh with 8 dof total. v is unknown vector in the problem.
Suppose I have 3 constraints \int v[0] dx =0, \int v[1] dx =0 and \int (x[1]*v[0]-x[0]*v[1])dx=0. I mapped the dof using interpolating function.
c1 = dolfinx.fem.Function(V)
c1.interpolate(lambda x: np.array([[1],[0]])) # Constraint $v1$ averaged to zero.
C1 = c1.vector
C1.scale(1 / C1.norm())
#
c2 = dolfinx.fem.Function(V)
c2.interpolate(lambda x: np.array([[0],[1]])) # Constraint $v[1]$ averaged to zero.
C2 = c2.vector
C2.scale(1 / C2.norm())
#
c3 = dolfinx.fem.Function(V)
c3.interpolate(lambda x: np.array((-x[1],x[0])) # Constraint $-x2*v[0]+v[1]*x1$ averaged to zero.
C3 = c3.vector
C3.scale(1 / C3.norm())
Error1a: On creating the nullspace, it is taking only vectors as argument, not even array. So, I made individual nulspaces as:
nullspace1 = petsc4py.PETSc.NullSpace().create(vectors=[C1], comm=mesh.comm)
nullspace2 = petsc4py.PETSc.NullSpace().create(vectors=[C2], comm=mesh.comm)
nullspace3 = petsc4py.PETSc.NullSpace().create(vectors=[C3], comm=mesh.comm)
Since, vectors argument is not taking (vectors = [[C1],[C2],[C3]]).
Error 1b: To make A.setNullspace
nullspace=np.hstack((nullspace1,nullspace2,nullspace3,nullspace4))
# Set the nullspace
A.setNullSpace(nullspace)
TypeError Traceback (most recent call last)
Cell In[292], line 2
1 # Set the nullspace
----> 2 A.setNullSpace(nullspace)
TypeError: Argument 'nsp' has incorrect type (expected petsc4py.PETSc.NullSpace, got numpy.ndarray)
Q: How can I give nullspace matrix (set of vectors) for A? (Since, tutorial is a scalar valued function)
Error2: I am not getting C1 vector as expected. For 4 nodes mesh, the first constraint of \int v[0] =0 should give the following interpolated vector as ([1 0 1 0 1 0 1 0 ]). However I am getting a different version.
C1.array
array([1.0e+000, 0.0e+000, 0.0e+000, 1.6e-322, 1.6e-322, 2.4e-322,
2.4e-322, 5.9e-323])
The strange thing is happening for C2. array as the 1 is coming for 2 components. whereas my expectations were [01 01 01 01].
C2.array
array([0.00000000e+000, 1.00000000e+000, 1.00000000e+000, 1.58101007e-322,
1.58101007e-322, 1.63041663e-322, 1.63041663e-322, 5.04443758e+180])
Can you clarify this weird dof mapping and comment on above procedure of implementing constraints?
I am greatly thankful for any help.