Region and labels in imported mesh


I made a 3D geometry using GMSH and exported it using .mesh format with it’s Physical entities.

load “msh3”
load “medit”
mesh3 Th = readmesh3(“MCore.mesh”);
func co = 1.0*(region==2);
fespace Vh(Th,P0);
Vh pco = co;

My question is why pco can not detect the half of the nodes on the boundary of entity. but co does.
How make pco to include the right and bottom boundary of (region==2).



seph.rar (190.8 KB)


Medit only knows how to plot P1 functions properly. If you use ParaView, e.g.:

int[int] fforder(1);
fforder = 0;
savevtk("test.vtu", Th, pco, order = fforder);

You’ll see that you have the correct region.
Screenshot 2020-06-17 at 9.45.46 AM

1 Like

Dear Professor

Thank you for your prompt and informative response.
I will use ParaView more often from now on.

but i am still a little bit confused.
for example on the same mesh, P1 and P2 refer to the points on the left and right side of the box with region=2.

cout << “Region P1 =” << Th(0.0055,-0.0025,0.0005).region <<endl;
cout << “Region P2 =” << Th(0.0095,-0.0025,0.0005).region <<endl;

will result in:
Region P1= 1
Region P2= 2

I thought both of them should be equal to 2

cout << “co P1 =” << co(0.0055,-0.0025,0.0005) <<endl;
cout << “co P2 =” << co(0.0095,-0.0025,0.0005) <<endl;

cout << “pco P1 =” << pco(0.0055,-0.0025,0.0005) <<endl;
cout << “pco P2 =” << pco(0.0095,-0.0025,0.0005) <<endl;

will result in

co P1 =0
co P2 =0

pco P1 =0
pco P2 =1

I can’t understand why pco and co are different although they look to be same when they are visualized with ParaView.

Even when I took a close look at ParaView by plotting the points:

It seems that pco returns different values on the different sides of the box.

The same result was obtained by GMSH.

ofstream ff(“graph.pos”);
ff << “View “pots” {” << endl;
for (int i = 0; i < Th.nv; i++)
ff << “SP(” << Th(i).x << “,” << Th(i).y << “,” << Th(i).z << “){” << pco(Th(i).x,Th(i).y,Th(i).z) << “};” << endl;
ff << “};” << endl;


The right side of the box has red dots. equivalent to pco=1
and the left side has pco=0

how can i set pco=1 on all the points belonging to region 2?

seph.rar (218.7 KB)

I’m guessing this comes from the interpolator. If you change to:

cout << "pco P2 =" << pco(0.00950001,-0.0025,0.0005) <<endl;

I get 0. There is notion of point here, you are using P0 functions, the unknowns are the values on the tetrahedra.