Thank you for the reply;

Yes, I am truncating from a bigger mesh. I am attaching a test version here. Th is my global mesh and Th1 my domain 1 then Th2 my domain 3.

load “iovtk”;

load “msh3”

mesh3 Th(“3sections.mesh”);

Th = readmesh3(“3sections.mesh”);

int LayerMax=1;

int totalSteps=5;

real cellx=(0.05715-(-0.01905))/9;

real celly=(0.0189484)/29;

real x0=-0.01905;

real y2=0.0189484-4*celly;*

real xMax=x0+4cellx;

func part1=(y>y2)*(x<xMax);

func part2=(y<=y2);

mesh3 Th1 = trunc(Th, part1 || part2 ,flabel=1,fregion=5);

real y2New=0.0189484-4*celly;*

real xMaxNew=x0+8cellx;

func part1New=(y>y2)*(x<xMax);

func part2New=(y<=y2);

mesh3 Th2 = trunc(Th, part1New || part2New ,flabel=1,fregion=5); ////// trunc the mesh depending on the region of interest by the function at 2nd position of the bracket,

//flabel for changing the surface tag and fregion for changing the volume tag

macro normal [N.x, N.y, N.z]//

macro fluxVECTOR(L) [0, L,0 ]//

real Tamb=300;

real dt=0.08;

real rho=4430;

real Cp=500;

real cond=10;

real h=10;

for (int i = 0; i < 1; i++)

```
{
```

fespace Vh(Th1,P1);

Vh Told=Tamb, T=Tamb, TT;

/////////////////////////// Varf Function ////////////////////////

varf vA(T, TT)

= int3d(Th1)

(

rho*Cp*T*TT/dt*

+ cond(dx(T)*dx(TT) + dy(T)*dy(TT) + dz(T)*dz(TT))

)

varf vB(T, TT) = int3d(Th1)(rho*Cp*T*TT/dt) ;*

varf vRHS(T,TT) = int2d(Th1, 1)(hTamb*TT) + int2d(Th1,1)( (fluxVECTOR(1.2e8*0.45*exp((-2*((x-(xMax-2*cellx))^2+(z-(1.5e-3))^2))/(2*cellx)^2) ) *(1-part2)) ’* normal*TT ) ;

real[int] rhs0 = vRHS(0, Vh);

matrix A = vA(Vh, Vh, solver=GMRES); ////// Stiffness Matrix

matrix B = vB(Vh, Vh); ////// Mass Matrix

real[int] b = B*T[];

b += rhs0;

```
T[] = A^-1*b;
```

int[int] Order = [1];

savevtk(“Varf.vtu”, Th1,T, dataname=“Temp”, order=Order);

`}`