Help with mesh compatibility

General Discussion
1

Sorry for doing a lot of questions. But I did the following:

// paraview plot
load "iovtk"
int NbTriangles = 8;

//------------------Building mesh--------------------
int localRef = 4;

//-----------------Global mesh-----------------------
mesh P = square(2,2);
plot(P, wait=true, fill=true);

//Local mesh - triangle 1
border a1(t=0, 0.5){x=t; y=0; label=1;};
border a2(t=0, 0.5){x=0.5; y=t; label=2;};
border a3(t=0, 0.5){x=t; y=x; label=3;};

mesh Th1 = buildmesh(a1(localRef)+ a2(localRef)+ a3(-localRef));
//plot(Th1, wait=true, fill=true);

//Local mesh - triangle 2
border b1(t=0, 0.5){x=0; y=t; label=2;};
border b2(t=0, 0.5){x=t; y=0.5; label=1;};
border b3(t=0, 0.5){x=t; y=x; label=3;};

mesh Th2 = buildmesh(b1(-localRef)+ b2(-localRef)+ b3(localRef));
//plot(Th2, wait=true, fill=true);

//Local mesh - triangle 3
border c1(t=0.5, 1.0){x=t; y=0; label=1;};
border c2(t=0, 0.5){x=1; y=t; label=2;};
border c3(t=0.5, 1.0){x=t; y=x-0.5; label=3;};

mesh Th3 = buildmesh(c1(localRef)+ c2(localRef)+ c3(-localRef));
//plot(Th3, wait=true, fill=true);

//Local mesh - triangle 4
border d1(t=0, 0.5){x=0.5; y=t; label=2;};
border d2(t=0.5, 1.0){x=t; y=0.5; label=1;};
border d3(t=0.5, 1.0){x=t; y=x-0.5; label=3;};

mesh Th4 = buildmesh(d1(-localRef)+ d2(-localRef)+ d3(localRef));
//plot(Th4, wait=true, fill=true);

//Local mesh - triangle 5
border e1(t=0.5, 1){x=t; y=0.5; label=1;};
border e2(t=0.5, 1){x=1; y=t; label=2;};
border e3(t=0.5, 1){x=t; y=x; label=3;};

mesh Th5 = buildmesh(e1(localRef)+ e2(localRef)+ e3(-localRef));
//plot(Th5, wait=true, fill=true);

//Local mesh - triangle 6
border f1(t=0.5, 1){x=t; y=1; label=1;};
border f2(t=0.5, 1){x=0.5; y=t; label=2;};
border f3(t=0.5, 1){x=t; y=x; label=3;};

mesh Th6 = buildmesh(f1(-localRef)+ f2(-localRef)+ f3(localRef));
//plot(Th6, wait=true, fill=true);

//Local mesh - triangle 7
border g1(t=0, 0.5){x=t; y=0.5; label=1;};
border g2(t=0.5, 1.0){x=0.5; y=t; label=2;};
border g3(t=0, 0.5){x=t; y=0.5+x; label=3;};

mesh Th7 = buildmesh(g1(localRef)+ g2(localRef) + g3(-localRef));
//plot(Th7, wait=true, fill=true);

//Local mesh - triangle 8
border h1(t=0.5, 1.0){x=0; y=t; label=2;};
border h2(t=0, 0.5){x=t; y=1; label=1;};
border h3(t=0, 0.5){x=t; y=0.5+x; label=3;};

mesh Th8 = buildmesh(h1(-localRef)+ h2(-localRef) + h3(localRef));
//plot(Th8, wait=true, fill=true);

//-------------Global+local mesh-----------------------
mesh Omega = Th1+Th2+Th3+Th4+Th5+Th6+Th7+Th8;
plot(Omega, wait=true, fill=true);

//----------------------------------------------------
//-------------Allocate local meshes------------------
mesh [int] Th(NbTriangles);

Th[0] = Th1;
Th[1] = Th2;
Th[2] = Th3;
Th[3] = Th4;
Th[4] = Th5;
Th[5] = Th6;
Th[6] = Th7;
Th[7] = Th8;

I what to compute the Laplace problem in each local refined mesh and then assemble to the global coarse mesh. But I don’t know if this is right. Can anyone help me?

//----------------------------------------------------
//Starting to compute Local problems

//def gradient
macro grad(u) [dx(u), dy(u)]//

//parameters
real kappa = 1.0;
int boundLabel = 1;

//RHS 
func f = 8*pi*pi*sin(2*pi*x)*sin(2*pi*y);

//dirichlet boundary condition
func g = 0;

//polinomial order
func Pklocal = P1;
func Pkglobal = P1;

fespace Vhglobal(P,Pkglobal);
Vhglobal [int] Thathf(NbTriangles); 
Vhglobal [int] ThlambdaH(NbTriangles);

for(int i=0; i<NbTriangles; i++){

 fespace Vhlocal(Th[i], Pklocal);

 Vhlocal Thl, Thf, vh;

 int ndoflocal = Vhlocal.ndof;
 
 //Local problem Thf
  varf aThf(Thf,vh)
   = int2d(Th[i])(kappa*(grad(Thf)'*grad(vh)));

  varf fh(Thf,vh) 
   = int2d(Th[i])(f*vh);

  matrix Af = aThf(Vhlocal, Vhlocal);

  //----------------------------------------------------------
  //Local problem ThlambdaH
  varf aThlambdaH(ThlambdaH,vh)
   = int2d(Th[i])(kappa*(grad(ThlambdaH)'*grad(vh)));

  
  varf lambdaH(ThlambdaH,vh) 
   = int1d(Th[i],boundLabel)(1.*vh);

  matrix Al = aThlambdaH(Vhlocal, Vhlocal);
  //-------------------------------------------------
  //Lagrange multiplier
  varf bh(Thf,vh)
  = int2d(Th[i])(1.*vh);

  real[int] bf(ndoflocal), blambda(ndoflocal);
  bf = fh(0,Vhlocal);
  blambda = lambdaH(0,Vhlocal);

  real[int] S = bh(0,Vhlocal);

  //block matrix
  matrix AAf = [ [ Af ,  S ] ,
                 [ S', 0 ] ] ;

  //block matrix
   matrix AAl = [ [ Al ,  S ] ,
                  [ S', 0 ] ] ;
 
  real[int] bbf(ndoflocal+1),xxThf(ndoflocal+1);
  real[int] bbl(ndoflocal+1),xxThlH(ndoflocal+1); 
  real[int] b1(1),l(1);
  b1=0;
 
  // build the block rhs 
  bbf = [ bf, b1];
  bbl = [ blambda, b1];
 
  set(AAf,solver=sparsesolver);
  xxThf = AAf^-1*bbf; // solve the linear system

  [Thf[],l] = xxThf;  // set the value

  set(AAl,solver=sparsesolver);
  xxThlH = AAl^-1*bbl; // solve the linear system

  [Thl[],l] = xxThlH;  // set the value


  //saving each solution on a global vector  
  ThlambdaH[i] = Thl;
  Thathf[i] = Thf;
  
  //Save as praview file
  int[int] Order = [1]; //scalar field
  string DataNamef = "Thf";
  savevtk("Thf_"+i+".vtk", Th[i], Thathf[i], dataname=DataNamef, order=Order);
  string DataNamel = "Thl";
  savevtk("Thl_"+i+".vtk", Th[i], ThlambdaH[i], dataname=DataNamel, order=Order);


  

}

Vhglobal globalSolution=0;

for(int i=0; i<NbTriangles; i++)
 globalSolution = globalSolution + ThlambdaH[i] + Thathf[i];

 plot(globalSolution, wait=true, fill=true);

//Save as praview file
  int[int] Order = [1]; //scalar field
  string DataName = "globalSol";
  savevtk("globalSol.vtk", P, globalSolution, dataname=DataName, order=Order);