How to work with regions on imported meshes (Gmsh)

Hi,

To define regions, I create them separately in a cad tool and mesh them in the Gmsh. By importing the mesh into FreeFem, I can distinguish the regions but FreeFem does not consider a unified domain and solve just one region or a part of the domain. Could you please provide a hint on this issue? How should we define regions by cad tools?
The code is written in the following lines. It simulates elastic deformation. The external mesh file “Part3.mesh” could not be uploaded since I am a new user. Thank you in advance.

load “msh3”
load “iovtk”
load “medit”

mesh3 Omega=readmesh3(“Part3.mesh”);

fespace Ph(Omega, P0);
Ph reg=region;
plot (reg,fill=true,value=true);

fespace Fes(Omega,P2);
Fes u1,u2,u3;
Fes v1,v2,v3;

func Fx=0.0;
func Fy=0.0;
func Fz=0.0;

func Tx=0.0;
func Ty=-0.01;
func Tz=0.0;// N/mm2

real Esoil=100;// N/mm2
real Edoor=100;// N/mm2
real nuSoil=0.3;
real nuDoor=0.3;

//========= Region ==============//
int soil =50;
int door = 49;
//==============================//
Ph E = Esoil * (region == soil) + Edoor * (region == door);
Ph nu = nuSoil * (region == soil) + nuDoor * (region == door);

int Top=51;
int Bottom=52;
int Sym=54;
int sides=53;
int intrface=56;

Ph miu=E/(2+2nu);
Ph lambda=E
nu/((1+nu)(1-2nu));

macro div(q1,q2,q3) (dx(q1)+dy(q2)+dz(q3))//
macro epsilon (q1,q2,q3) [dx(q1), dy(q2), dz(q3),sqrt(2)(dy(q1)+dx(q2)), sqrt(2)(dx(q3)+dz(q1))
,sqrt(2)(dz(q2)+dy(q3))]//
solve strLinear([u1,u2,u3],[v1,v2,v3])=
int3d(Omega)(lambda
div(u1,u2,u3)div(v1,v2,v3))
+int3d(Omega)(2.0
miu*(epsilon(u1,u2,u3)‘epsilon(v1,v2,v3)))
-int3d(Omega)([0.0,0.0,0.0]'
[v1,v2,v3])
-int2d(Omega,Top)([Tx,Ty,Tz]’*[v1,v2,v3])
+on(Bottom,u1=0.0,u2=0.0,u3=0.0)+on(sides,u1=0,u3=0)+on(Sym,u3=0);

savevtk(“Result.vtu”,Omega,u1,u2,u3,nu,E);
cout<<“==== End of the Simulation ===” <<endl;

I think you should try to define regular/physical regions that can be accessed in FreeFem. These depend on the format that you are using, see e.g., DMplex can not read region number