How to use Lagrange Multipliers in FreeFEM

General Discussion
#1

Hi all,

I am very new to the finite element field, and I want to try to solve the a variant of the fokker-plank equation with periodic boundary conditions. To do this, I need to find a way to impose an actual concentration in the system. Since I have a Robin boundary condition (no flux) at the surface of an inner wall, but periodic everywhere else, the problem does not have a unique solution. I’d like to impose some type of constraint that makes the integral over all space equal a constant or something like that to impose a unique solution. From my reading, it seems like the way to do that is through a lagrange multiplier, but I am not sure how to implement that in FreeFEM++.

For Example, here is how I would write and solve my equations when I had a neumann and a dirichlet boundary conditions and so it was unique (I still was periodic on just the z direction.)

fespace Vh(Th3 , P1, periodic=[[5,x,y],[6,x,y]]); 
Vh u, v // probability and test function.

problem FPeqn(prob3d,testfunc3d) = 
	int3d(Th3)( "bilinear term")
	+ int3d(Th3)("linear term")
	
            + int2d(Th3,innerlabel)("robin boundary term") // Robin BC
            + on(outerlabel, v = 1/(2*pi)) // Dirichlet BC 
;
   FPeqn

But from the brief example they have on the documentation https://doc.freefem.org/examples/finite-element.html#examplelagrangemultipliers , the way that they solve an equation with Lagrange Multipliers is very different than what I am used to. They use "varf"s instead of “problem” and they also make a block matrix.

I guess my question is, what is this block of code that they have. And how would I go about building that for a 3d system?

// Parameters
func f = 1 + x - y;

// Mesh
mesh Th = square(10, 10);

// Fespace
fespace Vh(Th, P1);
int n = Vh.ndof;
int n1 = n+1;
Vh uh, vh;

// Problem
varf va (uh, vh)
   = int2d(Th)(
          dx(uh)*dx(vh)
    + dy(uh)*dy(vh)
)
;

varf vL (uh, vh) = int2d(Th)(f*vh);
varf vb (uh, vh) = int2d(Th)(1.*vh);

matrix A = va(Vh, Vh);
real[int] b = vL(0, Vh);
real[int] B = vb(0, Vh);

// Block matrix
matrix AA = [ [ A, B ], [ B', 0 ] ];
set(AA, solver=sparsesolver);

real[int] bb(n+1), xx(n+1), b1(1), l(1);
b1 = 0;
// Builds the right hand side block
bb = [b, b1];

 // Solve
xx = AA^-1 * bb;

 // Set values
[uh[],l] = xx;

// Display
cout << " l = " << l(0) << " , b(u, 1) =" << B'*uh[] << endl;
#2

Maybe the documentation from FEniCS is cleaner? It solves exactly the same problem, so the variational formulations are the same. If you want to switch to 3D, switch mesh to mesh3 and int2d to int3d.

#3

That does help a bit thank you. It seems like the section

matrix A = va(Vh, Vh);
real[int] b = vL(0, Vh);
real[int] B = vb(0, Vh);
// Block matrix
  matrix AA = [ [ A, B ], [ B', 0 ] ];
set(AA, solver=sparsesolver);

real[int] bb(n+1), xx(n+1), b1(1), l(1);
b1 = 0;
// Builds the right hand side block
bb = [b, b1];

// Solve
xx = AA^-1 * bb;

// Set values
[uh[],l] = xx;

Is unique to FreeFEM. Is this block of code something that needs to be changed or is this just how I would extract the solutions to it?

#4

I don’t think you need to change that for your application.
The formulation in FreeFEM is really similar to what you would do in MATLAB, by defining a block matrix and concatenating vertically vectors or extracting scalars.

#5

Thank you for your help with this! After tinkering with it for a few days, I have another question if it isn’t too much of a bother.

I am trying to figure out what these sections are doing.

varf vL (uh, vh) = int2d(Th)(f*vh);
varf vb (uh, vh) = int2d(Th)(1.*vh);

I am under the impression here that f is the non homogeneous term in the original equation. If my equation is homogeneous, then can I set f = 0, or can I just delete this line with vL?
(if so, what do I do with b = vL(0, Vh); and bb = [b, b1] later on?)

Also, what is the vb term representing? Is that what the solution integrates to? For example, if my constraint is that int3d(u) = Constant, would I put that constant in the vb term above or would I put it something below like b1 = Constant?

matrix A = va(Vh, Vh);
real[int] b = vL(0, Vh);
real[int] B = vb(0, Vh);

// Block matrix
matrix AA = [ [ A, B ], [ B', 0 ] ];
set(AA, solver=sparsesolver);

real[int] bb(n+1), xx(n+1), b1(1), l(1);
b1 = 0;

I’m thinking I won’t have to change the block matrix stuff above, except maybe the b1 thing?

// Builds the right hand side block
bb = [b, b1];

// Solve
xx = AA^-1 * bb;

// Set values
[uh[],l] = xx;
#6

f is the source term, vb discretizes the constraint int(u) = b1.

#7

So if there is no source term, I can get away with f = 0?

Thank you!

#8

Why couldn’t you? Try it and you’ll see for yourself…

1 Like