Bad orientation

General Discussion
1

Dear colleages,

I am trying to solve a 3D problem on a thick slice with a pentagonal hole.

I have used gmsh to build up the mesh with the following code

cl1=0.5;
Mesh.Format=30;

Point(1)={-1,0,0,cl1};
Point(2)={1,0,0,cl1};
Point(3)={-1.8,-1,0,cl1};
Point(4)={1.8,-1,0,cl1};
Point(5)={0,-2,0,cl1};

Point(6)={-1,0,-5,cl1};
Point(7)={1,0,-5,cl1};
Point(8)={-1.8,-1,-5,cl1};
Point(9)={1.8,-1,-5,cl1};
Point(10)={0,-2,-5,cl1};

Point(11)={-5,3,0,cl1};
Point(12)={5,3,0,cl1};
Point(13)={-5,-15,0,cl1};
Point(14)={5,-15,0,cl1};

Point(15)={-5,3,-5,cl1};
Point(16)={5,3,-5,cl1};
Point(17)={-5,-15,-5,cl1};
Point(18)={5,-15,-5,cl1};

Line(1) = {3, 1};
Line(2) = {1, 2};
Line(3) = {2, 4};
Line(4) = {4, 5};
Line(5) = {5, 3};
Line(6) = {10, 9};
Line(7) = {9, 7};
Line(8) = {7, 6};
Line(9) = {6, 8};
Line(10) = {8, 10};
Line(11) = {10, 5};
Line(12) = {9, 4};
Line(13) = {7, 2};
Line(14) = {6, 1};
Line(15) = {8, 3};
Line(16) = {11, 12};
Line(17) = {12, 14};
Line(18) = {14, 13};
Line(19) = {13, 11};
Line(20) = {13, 17};
Line(21) = {17, 18};
Line(22) = {18, 16};
Line(23) = {16, 15};
Line(24) = {15, 17};
Line(25) = {18, 14};
Line(26) = {12, 16};
Line(27) = {15, 11};
Line Loop(1) = {5, -15, 10, 11};
Plane Surface(1) = {1};
Line Loop(2) = {11, -4, -12, -6};
Plane Surface(2) = {2};
Line Loop(3) = {12, -3, -13, -7};
Plane Surface(3) = {3};
Line Loop(4) = {13, -2, -14, -8};
Plane Surface(4) = {4};
Line Loop(5) = {14, -1, -15, -9};
Plane Surface(5) = {5};
Line Loop(10) = {23, 27, 16, 26};
Plane Surface(10) = {10};
Line Loop(11) = {22, -26, 17, -25};
Plane Surface(11) = {11};
Line Loop(12) = {21, 25, 18, 20};
Plane Surface(12) = {12};
Line Loop(13) = {27, -19, 20, -24};
Plane Surface(13) = {13};
Line Loop(14) = {16, 17, 18, 19};
Line Loop(15) = {2, 3, 4, 5, 1};
Plane Surface(14) = {14, 15};
Line Loop(16) = {22, 23, 24, 21};
Line Loop(17) = {6, 7, 8, 9, 10};
Plane Surface(15) = {16, 17};
Surface Loop(1) = {14, 10, 15, 11, 12, 13, 1, 5, 4, 3, 2};
Volume(1) = {1};

all look fine but when read the mesh with freefem

load “msh3”;
mesh3 Th(“pentagono.mesh”);

the error is "read mesh ok 0Mesh3, num Tetra:= 41384, num Vertice:= 8636 num boundary Triangles:= 7536 Mesh3::meshS, num Triangles:= 7536, num Vertice:= 3768 num boundary Edges:= 354
read mesh ok 0Mesh3, num Tetra:= 41384, num Vertice:= 8636 num boundary Triangles:= 7536
Mesh3::meshS, num Triangles:= 7536, num Vertice:= 3768 num boundary Edges:= 354
Bad orientation: The adj border element defined by [ 16 246 ] is oriented in the same direction in element 3922 and in the element 446 ****** bug in mesh construction? orientation parameter?
"

I do not know if freefem has the capability of handling with this.

Does someone deal with a problem like this?

Thanks in advance

Jose

2

Dear colleages,

I am still working on generating with GMSH and reading the mesh with FreeFem for a slice with a pentagonal hole:
cl1=0.5;
Mesh.Format=1;
// Select older mesh version, for compatibility of FF++
Mesh.MshFileVersion = 2.2;

Point(1)={-1,0,0,cl1};
Point(2)={1,0,0,cl1};
Point(3)={-1.8,-1,0,cl1};
Point(4)={1.8,-1,0,cl1};
Point(5)={0,-2,0,cl1};

Point(6)={-1,0,-5,cl1};
Point(7)={1,0,-5,cl1};
Point(8)={-1.8,-1,-5,cl1};
Point(9)={1.8,-1,-5,cl1};
Point(10)={0,-2,-5,cl1};

Point(11)={-5,3,0,cl1};
Point(12)={5,3,0,cl1};
Point(13)={-5,-15,0,cl1};
Point(14)={5,-15,0,cl1};

Point(15)={-5,3,-5,cl1};
Point(16)={5,3,-5,cl1};
Point(17)={-5,-15,-5,cl1};
Point(18)={5,-15,-5,cl1};

Line(1) = {3, 1};
Line(2) = {1, 2};
Line(3) = {2, 4};
Line(4) = {4, 5};
Line(5) = {5, 3};
Line(6) = {10, 9};
Line(7) = {9, 7};
Line(8) = {7, 6};
Line(9) = {6, 8};
Line(10) = {8, 10};
Line(11) = {10, 5};
Line(12) = {9, 4};
Line(13) = {7, 2};
Line(14) = {6, 1};
Line(15) = {8, 3};
Line(16) = {11, 12};
Line(17) = {12, 14};
Line(18) = {14, 13};
Line(19) = {13, 11};
Line(20) = {13, 17};
Line(21) = {17, 18};
Line(22) = {18, 16};
Line(23) = {16, 15};
Line(24) = {15, 17};
Line(25) = {18, 14};
Line(26) = {12, 16};
Line(27) = {15, 11};
Line Loop(1) = {5, -15, 10, 11};
Plane Surface(1) = {1};
Line Loop(2) = {11, -4, -12, -6};
Plane Surface(2) = {2};
Line Loop(3) = {12, -3, -13, -7};
Plane Surface(3) = {3};
Line Loop(4) = {13, -2, -14, -8};
Plane Surface(4) = {4};
Line Loop(5) = {14, -1, -15, -9};
Plane Surface(5) = {5};
Line Loop(10) = {23, 27, 16, 26};
Plane Surface(10) = {10};
Line Loop(11) = {22, -26, 17, -25};
Plane Surface(11) = {11};
Line Loop(12) = {21, 25, 18, 20};
Plane Surface(12) = {12};
Line Loop(13) = {27, -19, 20, -24};
Plane Surface(13) = {13};
Line Loop(14) = {16, 17, 18, 19};
Line Loop(15) = {2, 3, 4, 5, 1};
Plane Surface(14) = {14, 15};
Line Loop(16) = {22, 23, 24, 21};
Line Loop(17) = {6, 7, 8, 9, 10};
Plane Surface(15) = {16, 17};
Surface Loop(1) = {14, 10, 15, 11, 12, 13, 1, 5, 4, 3, 2};
Volume(1) = {1};

when I tried to read the resulting mesh with

// Level of output
verbosity = 10;
//
// Loading of libraries
//
load “msh3”
load “gmsh”;
load “iovtk”;
//load “medit”;
//
// Loading the mesh. It is hardly needed specifying the procedence with gmshload3
//
mesh3 Th3I=gmshload3(“pentagonV2.msh”);

the error message is

5665 vertices
closing file 24009 6078
24009 tetrahedrons
6078 triangles
30459 numElements
– GMesh3 , n V: 5665 , n Elm: 24009 , n B Elm: 6078mes 877 674.198 , bb: (-5 -15 -5) , (5 3 0)
– BuildAdj:nva= 3 4 6078
– warning true boundary element 102 is no in correct orientation
– warning true boundary element 103 is no in correct orientation
– warning true boundary element 104 is no in correct orientation
– warning true boundary element 105 is no in correct orientation
– warning true boundary element 106 is no in correct orientation
– warning true boundary element 107 is no in correct orientation
– warning true boundary element 108 is no in correct orientation
– warning true boundary element 109 is no in correct orientation
– warning true boundary element 110 is no in correct orientation
Warning change orientation of 1938 faces
Warning error in boundary oriention 1938 faces

These are my last efforts after consulting to the GMSH community, but I must confess I am quite lost.

If you have any suggestion, or suspect this could be a bug, any idea will be welcomed

Thank you in advance,

Jose

3

Here is how you can have FreeFEM reorient the mesh for you.

load "msh3";
lockOrientation = false;
mesh3 Th("pentagono.mesh");
lockOrientation = true;
mesh3 renumbered = Th;
fespace Vh(renumbered, P1);
Vh u = x + y + z;
load "medit"
medit("Th", Th, u);
3 Likes
4

Hi

as I got the following error message:

Fatal Error in mesh : RFQ.mesh number badly oriented element 532

I tried your proposal for getting Freefem reorienting the mesh. Unfortunaltely, the code does not fix it. Any idea?

thank you for your valuable help,

5

Without anything to reproduce the problem: no.

6

OK. I thought that maybe the error message would be enough to ring a bell.

Here is the case. I am trying to build a 3D mesh with gmsh made of a volume with multiple plane surfaces as a start in order to solve Poisson equation with Freefem. I have by the way no problem to do that with the kernel “built-in” in Gmsh and exporting a INRIA .mesh. Also I need at some point to use “OpenCascade” mesher in Gmsh as the geometry will get more complex.

Here is the Gmsh code:

SetFactory(“OpenCASCADE”);
//---------------------------
// Parameters of the geometry
//---------------------------

Mesh.CharacteristicLengthMin = .5;
Mesh.CharacteristicLengthMax = 10;

DIM = 3;
L = 3.;
R0 = 1.;
rho = 0.85;
res = 0.1;
Mesh.Optimize = 1;
p1[] += newp ; Point(newp) = { 0, 0, 0, res} ;
p1[] += newp ; Point(newp) = { R0, 0, 0, res} ;
p1[] += newp ; Point(newp) = { 0, R0, 0, res} ;
p1[] += newp ; Point(newp) = { rho, 3R0, 0, res} ;
p1[] += newp ; Point(newp) = { 3R0, 3R0, 0, res} ;
p1[] += newp ; Point(newp) = { 3R0, rho, 0, res} ;
p1[] += newp ; Point(newp) = { rho, R0+rho, 0, res} ;
p1[] += newp ; Point(newp) = { R0+rho, rho, 0, res} ;
p1[] += newp ; Point(newp) = { R0+rho, 0, 0, res} ;
p1[] += newp ; Point(newp) = { 0., R0+rho, 0, res} ;

p2[] += newp ; Point(newp) = { 0, 0, L, res} ;
p2[] += newp ; Point(newp) = { R0, 0, L, res} ;
p2[] += newp ; Point(newp) = { 0, R0, L, res} ;
p2[] += newp ; Point(newp) = { rho, 3R0, L, res} ;
p2[] += newp ; Point(newp) = { 3R0, 3R0, L, res} ;
p2[] += newp ; Point(newp) = { 3R0, rho, L, res} ;
p2[] += newp ; Point(newp) = { rho, R0+rho, L, res} ;
p2[] += newp ; Point(newp) = { R0+rho, rho, L, res} ;
p2[] += newp ; Point(newp) = { R0+rho, 0, L, res} ;
p2[] += newp ; Point(newp) = { 0., R0+rho, L, res} ;

//+
Line(1) = {11, 12};
//+
Line(2) = {11, 13};
//+
Line(3) = {14, 15};
//+
Line(4) = {15, 16};
//+
Line(5) = {16, 18};
//+
Line(6) = {14, 17};
//+
Line(7) = {5, 15};
//+
Line(8) = {4, 5};
//+
Line(9) = {4, 14};
//+
Line(10) = {4, 7};
//+
Line(11) = {7, 17};
//+
Line(12) = {2, 1};
//+
Line(13) = {3, 1};
//+
Line(14) = {8, 6};
//+
Line(15) = {6, 16};
//+
Line(16) = {6, 5};
//+
Line(17) = {1, 11};
//+
Line(18) = {13, 3};
//+
Line(19) = {3, 7};
//+
Line(20) = {17, 13};
//+
Line(21) = {2, 8};
//+
Line(22) = {18, 12};
//+
Line(23) = {18, 8};
//+
Line(24) = {12, 2};
//+

Line(25) = {11, 1};
//+
Line(26) = {1, 2};
//+
Line(27) = {2, 12};
//+
Line(28) = {12, 11};
//+
Line(29) = {7, 3};
//+
Line(30) = {17, 7};
//+
Line(31) = {13, 17};
//+
Line(32) = {3, 13};
//+
Line(33) = {8, 18};
//+
Line(34) = {16, 6};
//+
Line(35) = {15, 5};
//+
Line(36) = {5, 4};
//+
Line(37) = {12, 11};
//+
Line(38) = {17, 14};
//+
Line(39) = {13, 17};
//+
Line(40) = {13, 11};
//+
Line(41) = {15, 14};
//+
Line(42) = {16, 15};
//+
Line(43) = {18, 16};
//+
Line(44) = {12, 18};
//+
Line(45) = {5, 6};
//+
Line(46) = {7, 4};
//+
Line(47) = {1, 3};
//+
Line(48) = {8, 2};
//+
Line(49) = {6 ,8};
//+
Line(50) = {14 ,4};

//+
Curve Loop(1) = {13, 17, 2, 18};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {25, 26, 27, 28};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {29, 30, 31, 32};
//+
Plane Surface(3) = {3};
//+
Curve Loop(40) = {21, 33, 22, 24};
//+
Plane Surface(4) = {40};
//+
Curve Loop(5) = {5, 15, 14, 23};
//+
Plane Surface(5) = {5};
//+
Curve Loop(6) = {34, 4, 7, 16};
//+
Plane Surface(6) = {6};
//+
Curve Loop(7) = {35, 3, 9, 36};
//+
Plane Surface(7) = {7};
//+
Curve Loop(80) = {1, 40, 20, 6, 41, 42, 43, 44};
//+
Plane Surface(8) = {80};
//+
Curve Loop(9) = {45,8,46,19,47,12,48,49};
//+
Plane Surface(9) = {9};
//+
Curve Loop(10) = {38, 50, 10, 11};
//+
Plane Surface(10) = {10};

//+
Surface Loop(1) = {1, 6, 7, 2, 10, 8, 3, 4, 9, 5} Using Sewing;
//+
Volume(1) = {1};

//+
Show “*”;

for building the geometry to be meshed with Gmsh. No error message is sent by Gmsh and I can export a file RFQ.mesh (INRIA) mesh as I did with the kernel “built-in”. The FreeFem code is at the moment just

load “msh3”
load “medit”
mesh3 Th(“RFQ.mesh”);

and I got first

read mesh ok 0Mesh3, num Tetra:= 1061, num Vertice:= 694 num boundary Triangles:= 1146
Mesh3::meshS, num Triangles:= 1146, num Vertice:= 637 num boundary Edges:= 305
Err Border element not in mesh
0 : 69 70 251
1 : 73 74 252
2 : 70 251 254
3 : 70 71 254
4 : 76 77 253
5 : 69 251 253
6 : 74 252 254
7 : 76 251 253
8 : 71 252 254

I understood this was due to the presence of 4 points anticipated to be the centers for constructing circles later. These points are not in the mesh. I deleted these 4 vertices in RFQ.mesh directly (pt 9, 10, 19 and 20).

This solved the issue and I get now

the tet 0 is badly orienated tet -0.696641
the tet 2 is badly orienated tet -1.41456
the tet 4 is badly orienated tet -0.106155
the tet 6 is badly orienated tet -0.642099
the tet 8 is badly orienated tet -0.197834
the tet 11 is badly orienated tet -0.0435579
the tet 13 is badly orienated tet -0.0429028
the tet 15 is badly orienated tet -0.714966
the tet 16 is badly orienated tet -0.197144
the tet 20 is badly orienated tet -0.0885662
Fatal Error in mesh : RFQ.mesh number badly oriented element 661
nt 1292 / nbe 1146 / mes 2.72958 mesb 951.542
current line = 8
Assertion fail : (0)
line :988, in file …/femlib/Mesh3dn.cpp
Assertion fail : (0)
line :988, in file …/femlib/Mesh3dn.cpp
err code 6 , mpirank 0

Also, I did check the orientation of the surfaces with Gmsh. It sounds good as all surface normals point out:

RFQ_TEST-Line_OC
RFQ_TEST-Line_OC785Ă—681 21.1 KB

The same orientation went well when I used the kernel “built-in” in Gmsh.

Basically, you have all the story.

Thx

7

Surround your code between three single brackets

like this

otherwise we can’t copy/paste it easily, please.

8

Thx for the tip. Here are the 2 codes.

Gmsh code:

SetFactory(“OpenCASCADE”);
//---------------------------
// Parameters of the geometry
//---------------------------

Mesh.CharacteristicLengthMin = .5;
Mesh.CharacteristicLengthMax = 10;

DIM = 3;
L = 3.;
R0 = 1.;
rho = 0.85;
res = 0.1;
Mesh.Optimize = 1;
p1[] += newp ; Point(newp) = { 0, 0, 0, res} ;
p1[] += newp ; Point(newp) = { R0, 0, 0, res} ;
p1[] += newp ; Point(newp) = { 0, R0, 0, res} ;
p1[] += newp ; Point(newp) = { rho, 3 <em>R0, 0, res} ;
p1[] += newp ; Point(newp) = { 3</em> R0, 3 <em>R0, 0, res} ;
p1[] += newp ; Point(newp) = { 3</em> R0, rho, 0, res} ;
p1[] += newp ; Point(newp) = { rho, R0+rho, 0, res} ;
p1[] += newp ; Point(newp) = { R0+rho, rho, 0, res} ;
p1[] += newp ; Point(newp) = { R0+rho, 0, 0, res} ;
p1[] += newp ; Point(newp) = { 0., R0+rho, 0, res} ;

p2[] += newp ; Point(newp) = { 0, 0, L, res} ;
p2[] += newp ; Point(newp) = { R0, 0, L, res} ;
p2[] += newp ; Point(newp) = { 0, R0, L, res} ;
p2[] += newp ; Point(newp) = { rho, 3 <em>R0, L, res} ;
p2[] += newp ; Point(newp) = { 3</em> R0, 3 <em>R0, L, res} ;
p2[] += newp ; Point(newp) = { 3</em> R0, rho, L, res} ;
p2[] += newp ; Point(newp) = { rho, R0+rho, L, res} ;
p2[] += newp ; Point(newp) = { R0+rho, rho, L, res} ;
p2[] += newp ; Point(newp) = { R0+rho, 0, L, res} ;
p2[] += newp ; Point(newp) = { 0., R0+rho, L, res} ;

//+
Line(1) = {11, 12};
//+
Line(2) = {11, 13};
//+
Line(3) = {14, 15};
//+
Line(4) = {15, 16};
//+
Line(5) = {16, 18};
//+
Line(6) = {14, 17};
//+
Line(7) = {5, 15};
//+
Line(8) = {4, 5};
//+
Line(9) = {4, 14};
//+
Line(10) = {4, 7};
//+
Line(11) = {7, 17};
//+
Line(12) = {2, 1};
//+
Line(13) = {3, 1};
//+
Line(14) = {8, 6};
//+
Line(15) = {6, 16};
//+
Line(16) = {6, 5};
//+
Line(17) = {1, 11};
//+
Line(18) = {13, 3};
//+
Line(19) = {3, 7};
//+
Line(20) = {17, 13};
//+
Line(21) = {2, 8};
//+
Line(22) = {18, 12};
//+
Line(23) = {18, 8};
//+
Line(24) = {12, 2};
//+

Line(25) = {11, 1};
//+
Line(26) = {1, 2};
//+
Line(27) = {2, 12};
//+
Line(28) = {12, 11};
//+
Line(29) = {7, 3};
//+
Line(30) = {17, 7};
//+
Line(31) = {13, 17};
//+
Line(32) = {3, 13};
//+
Line(33) = {8, 18};
//+
Line(34) = {16, 6};
//+
Line(35) = {15, 5};
//+
Line(36) = {5, 4};
//+
Line(37) = {12, 11};
//+
Line(38) = {17, 14};
//+
Line(39) = {13, 17};
//+
Line(40) = {13, 11};
//+
Line(41) = {15, 14};
//+
Line(42) = {16, 15};
//+
Line(43) = {18, 16};
//+
Line(44) = {12, 18};
//+
Line(45) = {5, 6};
//+
Line(46) = {7, 4};
//+
Line(47) = {1, 3};
//+
Line(48) = {8, 2};
//+
Line(49) = {6 ,8};
//+
Line(50) = {14 ,4};

//+
Curve Loop(1) = {13, 17, 2, 18};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {25, 26, 27, 28};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {29, 30, 31, 32};
//+
Plane Surface(3) = {3};
//+
Curve Loop(40) = {21, 33, 22, 24};
//+
Plane Surface(4) = {40};
//+
Curve Loop(5) = {5, 15, 14, 23};
//+
Plane Surface(5) = {5};
//+
Curve Loop(6) = {34, 4, 7, 16};
//+
Plane Surface(6) = {6};
//+
Curve Loop(7) = {35, 3, 9, 36};
//+
Plane Surface(7) = {7};
//+
Curve Loop(80) = {1, 40, 20, 6, 41, 42, 43, 44};
//+
Plane Surface(8) = {80};
//+
Curve Loop(9) = {45,8,46,19,47,12,48,49};
//+
Plane Surface(9) = {9};
//+
Curve Loop(10) = {38, 50, 10, 11};
//+
Plane Surface(10) = {10};

//+
Surface Loop(1) = {1, 6, 7, 2, 10, 8, 3, 4, 9, 5} Using Sewing;
//+
Volume(1) = {1};

//+
Show “*”;

and here is the Freefem code for just reading the INRIA mesh:

load “msh3”
load “medit”
mesh3 Th(“RFQ.mesh”);

I also tried:

load “msh3”
load “medit”

lockOrientation = false;
mesh3 Th("RFQ.mesh");
lockOrientation = true;
mesh3 renumbered = Th;
fespace Vh(renumbered, P1);
Vh u = x + y + z;
medit("Th", Th, u);

but got the same result.

Thx

9

Can you try the code you copy/paste, please? I still can’t parse your .geo.
At least the first line should be:

SetFactory("OpenCASCADE");

and not:

SetFactory(“OpenCASCADE”);
10

The bad quoting is coming the forum editor as I did not first use the code insertion properly as you pointed out. In my .geo file it is properly quoted.

I have no apparent problem with .geo file at least with Gmsh that reads it and meshes it in 3D without error message.

Thx

11

Just in case I copy/paste again the .geo code:

SetFactory("OpenCASCADE");

//---------------------------
// Parameters of the geometry
//---------------------------

// ERROR: The mesh file doesn't contain vertices 
// => forgot to mesh

// Bad orientation: The adj border element  defined by [  89 91 ] 
// car hérité tel que la version bien orienté avec BI
// Correction de l'orientation par définition de lignes supplém. OpenCascade ne reconnait pas les -
// pour changer l'orientation. Les lignes supplémentaires sont des inverses de premières pour palier.
// Dans ce cas, les faces sont bien orientées mais le maillage 3D dans GMSH plante => PLC error.
// PLC Error: A segment and a facet intersect at point
// 

Mesh.CharacteristicLengthMin = .5;
Mesh.CharacteristicLengthMax = 10;

DIM = 3; 

L = 3.; 
R0 = 1.; 
rho = 0.85; 
res = 0.1;


Mesh.Optimize = 1;

p1[] += newp ; Point(newp) = {  0,  0, 0, res} ;
p1[] += newp ; Point(newp) = { R0,  0, 0, res} ;
p1[] += newp ; Point(newp) = { 0, R0, 0, res} ;
p1[] += newp ; Point(newp) = { rho, 3*R0, 0, res} ;
p1[] += newp ; Point(newp) = { 3*R0, 3*R0, 0, res} ;
p1[] += newp ; Point(newp) = { 3*R0,  rho, 0, res} ;
p1[] += newp ; Point(newp) = { rho, R0+rho, 0, res} ;
p1[] += newp ; Point(newp) = {  R0+rho, rho, 0, res} ;
p1[] += newp ; Point(newp) = {  R0+rho, 0, 0, res} ;
p1[] += newp ; Point(newp) = { 0., R0+rho, 0, res} ;


p2[] += newp ; Point(newp) = {  0,  0, L, res} ;
p2[] += newp ; Point(newp) = { R0,  0, L, res} ;
p2[] += newp ; Point(newp) = { 0, R0, L, res} ;
p2[] += newp ; Point(newp) = { rho, 3*R0, L, res} ;
p2[] += newp ; Point(newp) = { 3*R0, 3*R0, L, res} ;
p2[] += newp ; Point(newp) = { 3*R0,  rho, L, res} ;
p2[] += newp ; Point(newp) = { rho, R0+rho, L, res} ;
p2[] += newp ; Point(newp) = {  R0+rho, rho, L, res} ;
p2[] += newp ; Point(newp) = {  R0+rho, 0, L, res} ;
p2[] += newp ; Point(newp) = { 0., R0+rho, L, res} ;



//+
Line(1) = {11, 12};
//+
Line(2) = {11, 13};
//+
Line(3) = {14, 15};
//+
Line(4) = {15, 16};
//+
Line(5) = {16, 18};
//+
Line(6) = {14, 17};
//+
Line(7) = {5, 15};
//+
Line(8) = {4, 5};
//+
Line(9) = {4, 14};
//+
Line(10) = {4, 7};
//+
Line(11) = {7, 17};
//+
Line(12) = {2, 1};
//+
Line(13) = {3, 1};
//+
Line(14) = {8, 6};
//+
Line(15) = {6, 16};
//+
Line(16) = {6, 5};
//+
Line(17) = {1, 11};
//+
Line(18) = {13, 3};
//+
Line(19) = {3, 7};
//+
Line(20) = {17, 13};
//+
Line(21) = {2, 8};
//+
Line(22) = {18, 12};
//+
Line(23) = {18, 8};
//+
Line(24) = {12, 2};
//+

Line(25) = {11, 1};
//+
Line(26) = {1, 2};
//+
Line(27) = {2, 12};
//+
Line(28) = {12, 11};
//+
Line(29) = {7, 3};
//+
Line(30) = {17, 7};
//+
Line(31) = {13, 17};
//+
Line(32) = {3, 13};
//+
Line(33) = {8, 18};
//+
Line(34) = {16, 6};
//+
Line(35) = {15, 5};
//+
Line(36) = {5, 4};
//+
Line(37) = {12, 11};
//+
Line(38) = {17, 14};
//+
Line(39) = {13, 17};
//+
Line(40) = {13, 11};
//+
Line(41) = {15, 14};
//+
Line(42) = {16, 15};
//+
Line(43) = {18, 16};
//+
Line(44) = {12, 18};
//+
Line(45) = {5, 6};
//+
Line(46) = {7, 4};
//+
Line(47) = {1, 3};
//+
Line(48) = {8, 2};
//+
Line(49) = {6 ,8};
//+
Line(50) = {14 ,4};


//+
Curve Loop(1) = {13, 17, 2, 18};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {25, 26, 27, 28};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {29, 30, 31, 32};
//+
Plane Surface(3) = {3};
//+
Curve Loop(40) = {21, 33, 22, 24};
//+
Plane Surface(4) = {40};
//+
Curve Loop(5) = {5, 15, 14, 23};
//+
Plane Surface(5) = {5};
//+
Curve Loop(6) = {34, 4, 7, 16};
//+
Plane Surface(6) = {6};
//+
Curve Loop(7) = {35, 3, 9, 36};
//+
Plane Surface(7) = {7};
//+
Curve Loop(80) = {1, 40, 20, 6, 41, 42, 43, 44};
//+
Plane Surface(8) = {80};

//+
Curve Loop(9) = {45,8,46,19,47,12,48,49};
//+
Plane Surface(9) = {9};
//+
Curve Loop(10) = {38, 50, 10, 11};
//+
Plane Surface(10) = {10};



//+
Surface Loop(1) = {1, 6, 7, 2, 10, 8, 3, 4, 9, 5} Using Sewing;
//+
Volume(1) = {1};


//+
Show "*";
12

Copy the code you have in your .geo, not the one from the first post above, which is wrong.

Thanks, it’s OK now.

13

Naive question, but wouldn’t it be best to fix the .geo when not using OpenCASCADE? What the FreeFEM error says is that there are triangles in the 3D mesh which are not part of the boundary, which usually happens when there is an issue during mesh generation.

14

When I use the kernel “Built-in” instead of “OpenCASCADE”, I do not meet any problem. Gmsh meshes the geometry and the INRI .mesh that I export is read without problem by FreeFEM. I can solve the Poisson equation and get coherent results then.

15

OK. To be more accurate, if you just use the code I sent but use “built-in” kernel instead, Freefem still sees an orientation problem. If I mentioned to not have a problem it was with the older file where a good orientation was obtained with +/- sign in the command “Curve loop” in the .geo file. In this case, Freefem accepted the exported mesh.
It is when I tried to use OpenCASCADE instead that I noticed the following:

  • The obtained orientation with “Built-in” was changed by OpenCASCADE. OpenCASCADE in Gmsh apparently ignores orientation ajustment with +/-,

  • the only way I found to get the orientation I want was to create new lines with opposite directions. Then OpenCASCADE was able to see the same orientation. This is the file I sent. But this one is rejected by Freefem with an orientation issue even when I plot the normals, they are similar to the older case the Freefem accepted.

Later, I discovered the OpenCASCADE option “Using sewing” that creates a shell even if geometrical elements do not have the required topology. This helped me obtaining the mesh in Gmsh. What I have not tried yet is too use this option with the older file (without additional lines with opposite direction). I try this asap.

I hope it helps understanding.

Thx

16

OK. I just did the test. I still have an orientation issue:

Bad orientation: The adj border element defined by [ 89 91 ] is oriented in the same direction in element 300 and in the element 124 ****** bug in mesh construction? orientation parameter?