Lagrange multiplier and Dirichlet boundary condition

General Discussion
1

Hi all, I am trying to solve poroelasticity in a rectangular domain [0,1]\times[0,1]. The four sides are labeled as 1,2,3,4 and at the top boundary \Gamma (label 4), I have term that includes Lagrange Multiplier as <g, v> and at \Gamma I want to weakly impose u_p = u_s, where u_s is known. I wrote the following code:

problem poroelasticity([u1p,u2p,p,g1,g2], [v1p,v2p,q,t1,t2]) =
			int2d(Lh) ( 
						2*nup*((dx(u1p)*dx(v1p)) + 0.5*(dx(u2p) + dy(u1p))*(dx(v2p) + dy(v1p)) + (dy(u2p)*dy(v2p))) // 2nu(D(u),D(v))
					+ lambdap*(dx(u1p) + dy(u2p))*(dx(v1p) + dy(v2p)) // lambda(div(u),div(v))	
					- alphap*(p*dx(v1p) + p*dy(v2p)) // alpha(p,div(v))	
					+ (spp/dt)*p*q // s/dt(p,q)
					+ (alphap/dt)*(dx(u1p)*q+ dy(u2p)*q) // alpha/dt(div(u),q)
					+ kappinvp*(dx(p)*dx(q) + dy(p)*dy(q)) //kappa(\grad p, \grad q)
					)

		-	int2d(Lh)(
						f1p*v1p + f2p*v2p + f3p*q 
					+ (spp/dt)*pold*q 
					+ (alphap/dt)*(dx(u1pold)*q + dy(u2pold)*q)
					)
 	

		+	int1d(Lh,4)	(g1*v1p + g2*v2p)


		+	int1d(Lh,4)	(u1p*t1 + u2p*t2)

		-	int1d(Lh,4)	(u1s*t1 + u2s*t2)
 

		+ on(1,2,3, u1p=u1ptrue, u2p=u2ptrue, p = ptrue);// BC

However it keeps giving me the following error message:

 kk 2 6 : umfpack_di_numeric failed
  current line = 263
Exec error : umfpack_di_numeric failed
   -- number :1
Exec error : umfpack_di_numeric failed
   -- number :1

I would be grateful if someone could help me fix this issue. Thank you!

2

Dear Hemanta,
Your matrix is singular if g1,g2,t1,t2 are in a fespace over the whole mesh Lh (I don’t have this information from your post), because the value of g1,g2 is not determined inside the domain. A way to cure this is to add a small term
+int2d(Lh)(1.e-8*(g1*t1+g2*t2)) making the problem well-posed inside the domain, without perturbing your boundary formulation.
Otherwise a more clear formulation is to take g1,g2,t1,t2 defined only on the boundary \Gamma. An example doing this is here

as stokes_lagrange.edp
There the matrix is built by blocks. Probably you can do also by using a composite space.

3

Thank you Francois. Yes, I had defined g1, g2, t1, t2 over the whole mesh Lh. Your suggestion of using a small term is helping me now. I had gone through the discussion you mentioned (even before posting my approach here), but I felt very difficult to follow it. I will try to implement using composite space. Thank you again.