Local solution on global mesh

Hi there,

can anyone help me with ploting local solution on a global mesh?

This is my global mesh

int globalH = 2;
mesh calP = square(globalH,globalH);

then, I store the nodes coordinates to build local meshes in each global triangle

for(int K = 0; K < NbTriangles; K++){
 real coordX0, coordX1, coordX2;
 real coordY0, coordY1, coordY2;

 for(int i = 0; i< NbVertices; i++){
   if(elemNodesGlobal(K,0) == i){
    coordX0 = calP(i).x;
    coordY0 = calP(i).y;
   }
   if(elemNodesGlobal(K,1) == i){
    coordX1 = calP(i).x;
    coordY1 = calP(i).y;
   }
   if(elemNodesGlobal(K,2) == i){
    coordX2 = calP(i).x;
    coordY2 = calP(i).y;
   }
 }


 real[int, int] nodes = [[coordX0, coordY0], [coordX1,coordY1],
                          [coordX2, coordY2]];
 
 //cout << nodes << "\n";

 border f0(t=0,1){P.x=nodes(1,0)*t + nodes(0,0)*(1-t);P.y=nodes(1,1)*t + nodes(0,1)*(1-t); label=11;};
 border f1(t=0,1){P.x=nodes(2,0)*t + nodes(1,0)*(1-t);P.y=nodes(2,1)*t + nodes(1,1)*(1-t); label=22;};
 border f2(t=0,1){P.x=nodes(0,0)*t + nodes(2,0)*(1-t);P.y=nodes(0,1)*t + nodes(2,1)*(1-t); label=33;};

 Th[K] = buildmesh(f0(1) + f1(1) + f2(1));
 Th[K] = trunc(Th[K], abs(u0 -K) < 1e-5, split=localh);

 int NbTrianglesLocal = Th[K].nt;
 int[int] localLabels = [33,22,11];

then I prooced to compute local problems that are needed for my global problem an then assemble it to global matrices to solve the global problem.

Last, I need to gather the final solution that consists from local problems solutions and global problem.

//--------------------------Gathering the final approx solution

real globalSumL2 = 0, globalSumH1 = 0, globalSumH1semi = 0;
real[int, int] uHhg(NbTriangles,1000);


for(int K = 0; K < NbTriangles; K++){
 
 real[int] ThhatfK = Thf(K,:);

 real[int] u0hK(ThhatfK.n);
 for(int i = 0; i < u0hK.n; i++)
  u0hK(i) = u0h[K];

 real[int,int] ThlambdaHK(nDoFl0K,ThhatfK.n);
 for(int i = 0; i < nDoFl0K; i++){
  for(int j = 0; j < ThlambdaHK.m; j++){
   ThlambdaHK(i,j) = ThlambdaH((K*nDoFl0K)+i,j);
  }
 }
 real[int] sumThlambdaHK(ThlambdaHK.m);
 for(int i = 0; i < nDoFl0K; i++){
  int ii = locationMatrix(K,i);
  real[int] aux = lambdaH[][ii]*ThlambdaHK(i,:); 

  sumThlambdaHK = sumThlambdaHK + aux;
  
 }
 
 real[int] uHhK(u0hK.n);
 
 for(int i = 0; i < u0hK.n; i++)
  uHhK(i) = u0hK(i) + sumThlambdaHK(i) + ThhatfK(i); 

 fespace VhKaux(Th[K],Pk);
 VhKaux auxuHhK;
 auxuHhK[] = uHhK;
 //plot(auxuHhK, wait=1, fill=1);

 uHhg(K,:) =auxuHhK[];
 uHhg.resize(NbTriangles, u0hK.n);

 //Error
 VhKaux error = u - auxuHhK;

  L2error[K] = sqrt(int2d(Th[K])((error)^2));
  H1semierror[K] = sqrt(int2d(Th[K])((dx(auxuHhK)- dxu)^2)+int2d(Th[K])((dy(auxuHhK) - dyu)^2));
  H1error[K] = sqrt(L2error[K]^2 + H1semierror[K]^2);

  globalSumL2 = globalSumL2 + (L2error[K])^2;
  globalSumH1 = globalSumH1 + (H1error[K])^2;
  globalSumH1semi = globalSumH1semi + (H1semierror[K])^2;

  //Save as praview file
  int[int] Order = [1];
  string DataName = "u";
  savevtk("uh_"+K+".vtk", Th[K], auxuHhK, dataname=DataName, order=Order);


 
}

I can plot each local solution individually, but I don’t know how to plot as a whole.
For the mesh I just do the following:

mesh Omega;
for(int i = 0; i < NbTriangles; i++)
  Omega = Omega + Th[i];

plot(Omega, wait=1);

Is there a way to do this for the solution uHh? Til now all I managed to do was store it in a vector of vectors.