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;