In my code, I have a numerical solution that I computed inside a cube. I also have the exact solution. I want to calculate the error between the two over a subset region of the cube. To do this, I build a new mesh onto which I interpolate my solution. However, the numerically computed solution is not interpolating onto all the elements of my new mesh. Can someone see where I’m going wrong? (Please note, I’m not quite sure how to interpret the rmid, rup, rdown parameters)

I’ve included a screenshot of the plot of the interpolated solution and also a screen shot of a cut of the solution.

//Builds a mesh on a shell

////{

real cx1 = 0.5; // centre of big disc

real cy1 = 0.5;

real r1 = 0.7; // radius of big

real cx2 = 0.5; // centre of small disc

real cy2 = 0.5;

real r2 = myEpsPlus; // radius of small disc

real zmin=0.1,zmax=0.9;

int[int] rup=[1,1], rdown=[1,2], rmid=[1,3,2,4];

border C1(t1=0,2pi){x=cx1+r1cos(t1);y=cy1+r1sin(t1);label=1;};pi){x=cx1+r2

border C2(t1=0,2cos(t1);y=cy1+r2sin(t1);;label=2;};`ThAnulus= buildmesh(C1(20)+C2(-20)); //plot(Thcercle,wait =1); ThShell=buildlayers(ThAnulus, 20, zbound=[zmin,zmax], labelmid=rmid, labelup = rup, labeldown = rdown); plot(ThShell,wait=1); // fespace for solution and cut fespace VhShell(ThShell, P2); // u VhShell uu = u, uexact = exactSol; // u is the numerically computed solution plot(uu, wait = true); plot(uexact, wait = true); ////}`