Lagrange multipliers

Is there an example to solve the Laplace equantion with Neumann boundary conditions using lagrange multipliers but in a scenario that it’s not 1D? I did the follwing but I’m getting errors.

``````fespace VhK(Th[i], Pk);
fespace V0K(Th[i], P0);
VhK etaK, etaf, phiK, varphi;
V0K vj, vf;
int ndofVhK = VhK.ndof;
int ndofV0K = V0K.ndof;

varf ah(etaf,varphi)

varf fh(etaf,varphi)
= int2d(Th[i])(f * varphi);

varf bh(vf, varphi)
= int2d(Th[i])(1.*varphi);

matrix A = ah(VhK,VhK); //stiffness matrix

//real[int] B = bh(0,VhK);
matrix B = bh(V0K,VhK);

real[int] F(ndofVhK);
F = fh(0,VhK);

//the block matrix
matrix AA = [ [ A ,  B ] ,
[ B', 0 ] ] ;

real[int]  bb(ndofVhK+ndofV0K),xx(ndofVhK+ndofV0K),b1(ndofV0K),l(ndofV0K);
b1=0;

// build the block rhs
bb = [ F, b1];

set(AA,solver=UMFPACK);
xx = AA^-1*bb; // solve the linear system

[etaf[],l] = xx;  // set the value

plot(etaf,wait=1, fill=1);
``````

Before to build a FreeFem++ script write the mathematical formulation.

In yours exemple your matrix B == 0 because bh is a linear form.
To help you I need the Lagrangian with find the matrix AA

if the use the vector B it is OK , the number of unknown is ndofVhK+1 because a vector is see a a matrix with one column.

Those are the problems I want to implement on freefem

I are no problem, this like

``````mesh Th=square(10,10);
fespace VhK(Th, P2);
fespace V0K(Th, P0);
VhK eta,phi;
V0K v,psi;
func f= 16.*x^(1-x)*x^(1-y)*y;
solve Pb([eta,v],[phi,psi]) =
- int2d(Th)(f*phi) + on(1,2,3,4, eta = 0);
plot(eta,wait=1);
plot(v,wait=1);
``````

So, I have a global mesh

``````int globalH = 2;
mesh P = square(globalH,globalH);

int NbTriangles = P.nt; //number of triangles

``````

and in each triangle I created a local mesh, currently I’m not refining the local mesh, so it’s only one triangle

``````mesh[int] Th(NbTriangles);
int localh = 1;

for(int i = 0; i < NbTriangles; i++){
Th[i] = trunc(P, abs(u0 -i) < 1e-5, split=localh);

``````

and used this to number the edges

`````` fespace Lambda0(Th[i], P0edge);
int[int,int] edgeNumberLocal(Th[i].nt, Lambda0.ndofK);
for (int t=0; t < Th[i].nt; ++t){
for(int i=0; i< 3; ++i){ // in 2d
int nu = Lambda0(t,i) ;
edgeNumberLocal(t,i) = nu;

``````

now in each local mesh I want to compute the problems that I sent you. I want so solve for eta_j^K and eta_f^K. I did what you suggested

`````` fespace VhK(Th[i], Pk);
fespace V0K(Th[i], P0);
VhK etaK, etafK, phiK, varphi;
V0K vj, vf;
int ndofVhK = VhK.ndof;
int ndofV0K = V0K.ndof;

solve localProblemetafK([etafK, vf], [varphi, phiK])
- int2d(Th[i])(f*varphi) + on(0,1,2, etafK = 0);

plot(etafK, wait=1);
plot(vf, wait = 1);
``````

and when I do this, I get the following error

`````` current line = 186
Internal error : Methode de Galerkine (a faire)
line  :11400, in file problem.cpp
Internal error : Methode de Galerkine (a faire)
line  :11400, in file problem.cpp
err code 7 ,  mpirank 0
``````

Also, is it right to do this?

`````` int[int] labelChange = [0,0, 1,1, 2,2];
Th[i] = change(Th[i], label = labelChange);

``````

Sorry, some fonction are no implemented.
I will try to do rapidly, but ii is not sURE;

1 Like