Lagrange multipliers

General Discussion
1

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)
    = int2d(Th[i])(kappa*grad(etaf)'*grad(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);
3

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.

4

Those are the problems I want to implement on freefem

formulation

5

I are no problem, this like

mesh Th=square(10,10);
macro grad(u) [dx(u),dy(u)]//
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)( grad(eta)'*grad(phi) + v*psi + eta*phi )
- int2d(Th)(f*phi) + on(1,2,3,4, eta = 0);
plot(eta,wait=1);
plot(v,wait=1);
6

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])(kappa*grad(etafK)'*grad(varphi) + vf*phiK + etafK*varphi)
      - 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
7

Also, is it right to do this?

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

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

1 Like