How to use Lagrange Multipliers in FreeFEM

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;``````
1 Like

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`.

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?

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.

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;``````

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

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

Thank you!

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

1 Like