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