I have a domain like the one shown in the picture. I need to solve an equation where k = 1 inside the squares and 0 outside (ignore the red border). I know it is possible to solve this in multiple ways and I actually solved it defining the value of k using k=a*(x<b)*(x>c)… for each square.

My question is, is it possible to somehow define the region inside the squares like only one region (something like k=1*(region==squares)) and apply the value of k only one time to the hole area? I also tried to build a mesh with the squares and use that mesh to define the value of k, but i couldn’t make it work.

Remark , by default the region a different on each region.

So if you put all centre of square in a array cs(2,nsquare)

you can get the region numler of each square

here the mesh is call TH

int[int] regs= regions(Th);
int[int] labs= labels(Th);
int mxreg = regs.max+1;
assert(merge < 1000); // no too big
int[int] regs(nsquare);
for(int i=0; i<nsquare;++i)
regs[i] = Th(cs(0,i),cs(1,i)).region ; // get the region number of each square;
real[int] K(mxreg);
K= 2; // out of square;
for(int i=0; i<nsquare;++i)
K[regs[i]] =1.; // put 1 on each square
// new to use K
fespace Ph(Th,P0); // function constant per Triangle
Ph Kk=K[region];
plo(Kk, fill=1, wait=1);
Vh u,v;
solve Lap(u,v) ( K[region]*grad(u)'*grad(v)- +on(labs,u=1);
// region depend directly of the place of the current point .
//or
solve Lap1(u,v) ( Kk*grad(u)'*grad(v)- +on(labs,u=1);