Constrains using null space in dolfinx

[bagla0]

Hi,
I am implementing nullspace to apply constrains in alternative to lagrange multiplier in dolfinx. The MWE is:

from mpi4py import MPI
from dolfinx import mesh
import numpy as np
from ufl import Jacobian, as_vector, dot, cross,sqrt, conditional, replace
from ufl import lt,SpatialCoordinate, as_tensor, VectorElement, FiniteElement, MixedElement, Measure
mesh = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
Em = 50e3
num = 0.2
Er = 210e3
nur = 0.3
material_parameters = [(Em, num), (Er, nur)]
nphases = len(material_parameters)
def eps_te(v):
    E1=as_vector([0,v[2].dx(0),0,0,v[1].dx(0),v[0].dx(1)])
    return as_tensor([(E1[0],0.5*E1[5],0.5*E1[4]),(0.5*E1[5],E1[1],0.5*E1[3]),(0.5*E1[4],0.5*E1[3],E1[2])]),E1

def sigma(v,i,Eps):     
        E,nu=material_parameters[i]          
        lmbda = E*nu/(1+nu)/(1-2*nu)
        mu = E/2/(1+nu)
        C1=lmbda+2*mu
        C=as_tensor([(C1,lmbda,lmbda,0,0,0),(lmbda,C1,lmbda,0,0,0),(lmbda,lmbda,C1,0,0,0),(0,0,0,mu,0,0),(0,0,0,0,mu,0),(0,0,0,0,0,mu)])        
        s1= dot(C,eps_te(v)[1]+Eps)
        return as_tensor([(s1[0],s1[5],s1[4]),(s1[5],s1[1],s1[3]),(s1[4],s1[3],s1[2])]), C
from dolfinx.fem import FunctionSpace, form, Function
import numpy as np
from ufl import TrialFunction, TestFunction, inner, form, lhs, rhs, dx
Ve = VectorElement("S", mesh.ufl_cell(), 2,dim=3)
V = FunctionSpace(mesh, Ve)
dv = TrialFunction(V)
v_ = TestFunction(V)
Eps= as_vector((1,0,0,0,0,0))
dx = Measure('dx')(domain=mesh)
F2 = sum([inner(sigma(dv, 0, Eps)[0], eps_te(v_)[0])*dx]) 

a=lhs(F2)
L=rhs(F2)

1. I am getting error to delete cache file to work lhs. After deleting cache file, it starts working. Do I need to do this everytime?
2. Getting error in making ufl.form

a_cpp = form(a)
f_cpp = form(L)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[20], line 3
      1 a=lhs(F2)
      2 L=rhs(F2)
----> 3 a_cpp = form(a)
      4 f_cpp = form(L)

TypeError: 'module' object is not callable

3.a Based on singular poisson example using null space where the constraint was \int u dx =0 . And
below code is used. Does, 1 shown in (C.array,1.0) represents any arbitrary constant coming after integration.

Can you explain in detail where does the 3 constraints are applied (for u1,u2,u3 averaged to 0 over entire dx) ?
3b. I got error in c.interpolate. I don’t understand its significance of using x.shape. What should be for this example?

c = Function(V)
c.interpolate(lambda x: np.ones(x.shape[1]))
C = c.vector
assert np.allclose(C.array, 1.0)  # because V is a Lagrange space

4. if I want a constraint like \int dv.dx(0) =0, how can I implement this ?

5. how can I implement periodic boundary condition. Previously, we give argument in defining function space of mixedelement constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10). I don’t know where should I implement it?
   Any help is appreciated in using nullspace.

[dokken]

@bagla0

I am getting error to delete cache file to work lhs. After deleting cache file, it starts working. Do I need to do this everytime?

Unfortunately, I cannot reproduce this behaviour.

@bagla0

Getting error in making ufl.form

You are using the wrong form.
You should use dolfinx.fem.form, not ufl.form.
Remove form from

@bagla0

from ufl import TrialFunction, TestFunction, inner, form, lhs, rhs, dx

@bagla0

3b. I got error in c.interpolate. I don’t understand its significance of using x.shape. What should be for this example?

You have a function space V that is a vector function space.
The input to interpolate is a set of coordinates, x, of shape (3, number_of_points), where x[0] is the list of all x-coordinates, x[1] all y coordinates etc.
The input is designed this way so that you can used vectorized functions in numpy.
As your space is a vector space of dim=3, the output should be (3, number_of_points), not just number_of_points.

@bagla0

if I want a constraint like \int dv.dx(0) =0∫dv.dx(0)=0\int dv.dx(0) =0, how can I implement this ?

Not sure why you would want such a constraint, what is the physical interpretation of this.

@bagla0

how can I implement periodic boundary condition. Previously, we give argument in defining function space of mixedelement constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10). I don’t know where should I implement it?
Any help is appreciated in using nullspace.

This is a topic covered in many other posts on discourse, see for instance: Periodic BC in dolfinx - #2 by nate

[francesco-ballarin]

@bagla0

Hi,
I am implementing nullspace to apply constrains in alternative to lagrange multiplier in dolfinx. The MWE is:

IMHO, @bagla0 you should have led with a bit more of context. By previous interactions, I guess that these questions are related to Mixed function space with dolfinx - how to implement an integral constraint - #6 by bagla0, but any other reader may not be aware of this.

@dokken

@bagla0

3b. I got error in c.interpolate. I don’t understand its significance of using x.shape. What should be for this example?

You have a function space V that is a vector function space.
The input to interpolate is a set of coordinates, x, of shape (3, number_of_points), where x[0] is the list of all x-coordinates, x[1] all y coordinates etc.
The input is designed this way so that you can used vectorized functions in numpy.
As your space is a vector space of dim=3, the output should be (3, number_of_points), not just number_of_points.

In addition to what @dokken is saying, keep in mind that should have two rigid body motions (in 2D) associated to a translation. A basis could be (1, 0) for displacement in the x-direction and (0, 1) for displacement in the y-direction. Therefore, rather than

c.interpolate(lambda x: np.ones(x.shape[1]))

you will have to come up with a lambda function that is non-zero only on one component of the displacement. To do that, follow the suggestions above.

@dokken

@bagla0

if I want a constraint like \int dv.dx(0) =0∫dv.dx(0)=0\int dv.dx(0) =0, how can I implement this ?

Not sure why you would want such a constraint, what is the physical interpretation of this.

Not sure either. Instead, @bagla0 you definitely need to add a third (in 2D) vector in null space to prevent rigid rotations. That’s not really my field, but I guess the third element in the null space could be something like (-y, x) in 2D. Maybe @bleyerj can confirm, or fix my suggestion if it is wrong. Depending on what the third vector looks like, you may need to carry out an orthonormalization with dolfinx.la.orthonormalize before passing the three vectors to PETSc.

[bagla0]

Thanks @dokken and @francesco-ballarin for your response. Sorry, I am still having query from interpolation step. I am trying to implement nullspace method of applying constraints.

1. I am not able to understand why we are interpolating if we want to apply constraint over entire mesh.
   Why do we create c function and its interpolation? In general, interpolation means, we give expression as input.

@bagla0

c = Function(V)
c.interpolate(lambda x: np.ones(x.shape[1]))
C = c.vector
assert np.allclose(C.array, 1.0)  # because V is a Lagrange space

2. Can you give the steps in algorithm form from creating interpolation to using it for nullspace and tell the reason for using assert(C.array,1.0) ?
   I have some understanding that we give solving function u as u_o +Z v for null space where u_o is solution of constraint applied, Z as basis of null space and v as vector of dimension (no of constraints X1). I might be wrong, but I am trying to relate theory with process.

@dokken

Not sure why you would want such a constraint, what is the physical interpretation of this.

in beam theory, we require displacement derivatives averaged to zero. I have doubt of do we need to apply this constraint through c.interpolate and where does this constraint applied. If possible can you give in example expression form.

c.interpolate(mesh.geometry.x) 

    397     _interpolate(u, cells)
    398 except TypeError:
    399     # u is callable
--> 400     assert callable(u)


AssertionError: 

In Mixed function space with dolfinx - how to implement an integral constraint - #6 by bagla0, the domain without decomposition uses u_average and u_substract to implement \int u dx=0 after solving u through solver.

4.I couldn’t understand if we are applying constraints using null space, why do we need to subtract the average. Does that mean, no constrains was applied in entire code and we are just subtracting average at end.

[francesco-ballarin]

I am not able to understand why we are interpolating if we want to apply constraint over entire mesh.
Why do we create c function and its interpolation? In general, interpolation means, we give expression as input.

In my original tutorial, where the solution was a scalar-valued function, we are trying to obtain the FE degrees of freedom associate to the function which is identically equal to 1. Therefore, I use c.interpolate(lambda x: np.ones(x.shape[1])) where lambda x: np.ones(x.shape[1]) is just a way of saying that the function is constant and equal to one.
In your case, you will have two cases: (1,0) and (0,1), since the solution is now a vector-valued function (rather than a scalar-valued one).

Can you give the steps in algorithm form from creating interpolation to using it for nullspace

That’s what the code above is doing

tell the reason for using assert(C.array,1.0)

That’s just me being overzealous. For a Lagrange space, you could have avoided interpolation, and realized that to get the scalar-valued function identically equal to 1 it would have been enough to set all DOFs equal to one.

in beam theory, we require displacement derivatives averaged to zero

Sorry, as discussed in the other post, you can’t implement that. If that is a way of recovering uniqueness of the solution, drop that the constraint and set the appropriate null space instead.

c.interpolate(mesh.geometry.x)

Interpolate requires a lambda function, e.g.

c.interpolate(lambda x: x) 

I couldn’t understand if we are applying constraints using null space, why do we need to subtract the average. Does that mean, no constrains was applied in entire code and we are just subtracting average at end.

I think I have explained this in my tutorial: it’s not necessary to subtract the average, if you want (e.g., because you are comparing to a FEniCS legacy code that has that constraint) you can do it, otherwise just skip it.

[bagla0]

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.

@bagla0

dv = TrialFunction(V)

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.

[francesco-ballarin]

Try with something like the following

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())

dolfinx.la.orthogonalize([C1, C2, C3])

nullspace1 = petsc4py.PETSc.NullSpace().create(vectors=[C1, C2, C3], comm=mesh.comm)

I still believe you are thinking about this the wrong way. You must stop thinking about this using the word “constraints”, and start thinking of it with “null space”. The comments next to the calls to .interpolate don’t make any mathematical sense: for instance, the first component of the solution (comment next to c1) will not necessary average to zero.

Please discuss this with your supervisor, a professor at your department, or a senior colleague, to make sure you truly understand what you are doing.

[bagla0]

Thanks @francesco-ballarin for your response. I just find a small typo of orthonormalize: dolfinx.la.orthonormalize([C1, C2, C3]) and above section worked.

Can you suggest some references to have better understanding of nullspace procedure if possible. I will try to dig deep into its theoretical implement where I am going wrong.

@bagla0

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]) 

@bagla0

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 tell briefly about the C1.array and C2.array that I got based on my (1,0) and (0,1) input?

[bagla0]

I am getting difficulty to understand mathematics behind constraint implementation using nullspace. I tried looking for resources but couldn’t get that. Can you refer some resources where I could find some examples? Any help is greatly appreciated.

Thanks

[dokken]

I would suggest reading: KSP: Linear System Solvers — PETSc v3.20.3-458-ge4cb84850f3 documentation

[francesco-ballarin]

I tried explaining it as best as I could in tutorial_create_nullspace . Note that until you get to “With domain decomposition” you don’t really need any additional library (i.e., multiphenicsx).