1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4 
5     This program is free software: you can redistribute it and/or modify
6     it under the terms of the GNU General Public License as published by
7     the Free Software Foundation, either version 3 of the License, or
8     (at your option) any later version.
9     This program is distributed in the hope that it will be useful,
10     but WITHOUT ANY WARRANTY; without even the implied warranty of
11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12     GNU General Public License for more details.
13     You should have received a copy of the GNU General Public License
14     along with this program.  If not, see <http://www.gnu.org/licenses/>.
15  */
16 
17 /*!
18   \file unit_Mesh_gmsh.cpp
19   \author N. Kielbasiewicz
20   \since 13 dec 2013
21   \date 5 mar 2014
22 
23   Low level tests of Mesh class using system call to gmsh
24 */
25 
26 
27 //===============================================================================
28 //library dependences
29 #include "xlife++-libs.h"
30 #include "testUtils.hpp"
31 
32 //===============================================================================
33 //stl dependences
34 #include <iostream>
35 #include <fstream>
36 #include <vector>
37 
38 //===============================================================================
39 
40 using namespace xlifepp;
41 
42 namespace unit_Mesh_gmsh {
43 
unit_Mesh_gmsh(bool check)44 String unit_Mesh_gmsh(bool check)
45 {
46   String rootname = "unit_Mesh_gmsh";
47 #ifdef XLIFEPP_WITH_GMSH
48   trace_p->push(rootname);
49   std::stringstream out;                  //string stream receiving results
50   out.precision(testPrec);
51   char s[200];
52   Number it_l;
53   std::ifstream fin;
54   verboseLevel(2);
55   setDefaultCharacteristicLength(0.5);
56 
57   // Mesh built from a .geo file, input file for Gmsh, which is called to create the Gmsh format file
58   std::cout << "--- Testing loadGeo" << std::endl;
59   Mesh xlifepp2D(inputsPathTo("xlifepp2D.geo"), "xlifepp2D", _geo);
60   Mesh xlifepp2D3D(inputsPathTo("xlifepp2D3D.geo"), "xlifepp2D3D", _geo);
61   Mesh xlifepp3D(inputsPathTo("xlifepp3D.geo"), "xlifepp3D", _geo);
62 
63   out << "*********************************************" << std::endl;
64 
65   std::cout << "--- Segment with sides" << std::endl;
66   out << "--- Segment with sides" << std::endl << std::endl;
67   Mesh mesh1dFromGmsh(Segment(_v1=Point(0.,0.,0.), _v2=Point(1.,2.,3.), _nnodes=6, _side_names=Strings("Gamma_1","Gamma_2")),
68                       1, _gmsh);
69   fin.open("xlifepp_script.geo");
70   it_l=0;
71   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
72   fin.close();
73 
74   out << "*********************************************" << std::endl;
75 
76   std::cout << "--- Rectangle with triangles" << std::endl;
77   out << "--- Rectangle with triangles" << std::endl << std::endl;
78   Mesh mesh2dP1(Rectangle(_xmin=0.,_xmax=2.,_ymin=0.,_ymax=4.,_nnodes=6,
79                           _side_names=Strings("Gamma_1","Gamma_2","Gamma_3","Gamma_4")), _triangle, 1, _gmsh);
80   fin.open("xlifepp_script.geo");
81   it_l=0;
82   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
83   fin.close();
84 
85   out << std::endl << "*********************************************" << std::endl;
86 
87   std::cout << "--- Rectangle with structured triangles" << std::endl;
88   out << "--- Rectangle with structured triangles" << std::endl << std::endl;
89   Mesh mesh2dP12(Rectangle(_xmin=0.,_xmax=2.,_ymin=0.,_ymax=4.,_nnodes=6,
90                           _side_names=Strings("Gamma_1","Gamma_2","Gamma_3","Gamma_4")), _triangle, 1, _gmsh, _structuredMesh);
91   fin.open("xlifepp_script.geo");
92   it_l=0;
93   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
94   fin.close();
95 
96   out << std::endl << "*********************************************" << std::endl;
97 
98   std::cout << "--- Rectangle with quadrangles" << std::endl;
99   out << "--- Rectangle with quadrangles" << std::endl << std::endl;
100   Mesh mesh2dQ1(Rectangle(_xmin=0.,_xmax=2.,_ymin=0.,_ymax=4.,_nnodes=6,
101                           _side_names=Strings("Gamma_1","Gamma_2","Gamma_3","Gamma_4")), _quadrangle, 1, _gmsh);
102   fin.open("xlifepp_script.geo");
103   it_l=0;
104   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
105   fin.close();
106 
107   out << std::endl << "*********************************************" << std::endl;
108 
109   std::cout << "--- Rectangle with segments" << std::endl;
110   out << "--- Rectangle with segments" << std::endl << std::endl;
111   Mesh mesh2dSeg(Rectangle(_xmin=0.,_xmax=2.,_ymin=0.,_ymax=4.,_nnodes=6,
112                            _side_names=Strings("Gamma_1","Gamma_2","Gamma_3","Gamma_4")), _segment, 1, _gmsh);
113   fin.open("xlifepp_script.geo");
114   it_l=0;
115   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
116   fin.close();
117 
118   out << std::endl << "*********************************************" << std::endl;
119 
120   std::cout << "--- Cuboid with tetrahedra" << std::endl;
121   out << "--- Cuboid with tetrahedra" << std::endl << std::endl;
122   Mesh mesh3dP1(Cuboid(_xmin=0.,_xmax=2.,_ymin=0.,_ymax=4.,_zmin=0.,_zmax=1.,_nnodes=6,
123                        _side_names=Strings("Gamma_1","Gamma_2","Gamma_3","Gamma_4","Gamma_5","Gamma_6")), _tetrahedron, 1, _gmsh);
124   fin.open("xlifepp_script.geo");
125   it_l=0;
126   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
127   fin.close();
128 
129   out << std::endl << "*********************************************" << std::endl;
130 
131   std::cout << "--- Cuboid with hexahedra" << std::endl;
132   out << "--- Cuboid with hexahedra" << std::endl << std::endl;
133   Mesh mesh3dQ1(Cuboid(_xmin=0.,_xmax=2.,_ymin=0.,_ymax=4.,_zmin=0.,_zmax=1.,_nnodes=6,
134                        _side_names=Strings("Gamma_1","Gamma_2","Gamma_3","Gamma_4","Gamma_5","Gamma_6")), _hexahedron, 1, _gmsh);
135   fin.open("xlifepp_script.geo");
136   it_l=0;
137   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
138   fin.close();
139 
140   out << std::endl << "*********************************************" << std::endl;
141 
142   std::cout << "--- Cuboid with triangles" << std::endl;
143   out << "--- Cuboid with triangles" << std::endl << std::endl;
144   Mesh mesh3dTri(Cuboid(_xmin=0.,_xmax=2.,_ymin=0.,_ymax=4.,_zmin=0.,_zmax=1.,_nnodes=6,
145                         _side_names=Strings("Gamma_1","Gamma_2","Gamma_3","Gamma_4","Gamma_5","Gamma_6")), _triangle, 1, _gmsh);
146   fin.open("xlifepp_script.geo");
147   it_l=0;
148   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
149   fin.close();
150 
151   out << std::endl << "*********************************************" << std::endl;
152 
153   std::cout << "--- pentagon with triangles" << std::endl;
154   out << "--- pentagon with triangles" << std::endl << std::endl;
155   std::vector<Point> vertices(5), vertices2(4);
156   vertices[0]=Point(0.,0.,0.);
157   vertices[1]=Point(2.,0.,0.);
158   vertices[2]=Point(3.,1.,0.);
159   vertices[3]=Point(1.,4.,0.);
160   vertices[4]=Point(-1.,2.,0.);
161   Polygon pg1(_vertices=vertices, _nnodes=3);
162   Mesh meshPolygon(pg1, _triangle, 1, _gmsh);
163   fin.open("xlifepp_script.geo");
164   it_l=0;
165   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
166   fin.close();
167 
168   out << std::endl << "*********************************************" << std::endl;
169 
170   std::cout << "--- Polyhedron with triangles" << std::endl;
171   out << "--- Polyhedron with triangles" << std::endl << std::endl;
172   vertices[0]=Point(0.,0.,1.);
173   vertices[1]=Point(2.,0.,1.);
174   vertices[2]=Point(3.,1.,1.);
175   vertices[3]=Point(1.,4.,1.);
176   vertices[4]=Point(-1.,2.,1.);
177   Polygon pg2(_vertices=vertices, _nnodes=3);
178   vertices2[0]=Point(0.,0.,0.);
179   vertices2[1]=Point(2.,0.,0.);
180   vertices2[2]=Point(2.,0.,1.);
181   vertices2[3]=Point(0.,0.,1.);
182   Polygon pg3(_vertices=vertices2, _nnodes=3);
183   vertices2[0]=Point(2.,0.,0.);
184   vertices2[1]=Point(3.,1.,0.);
185   vertices2[2]=Point(3.,1.,1.);
186   vertices2[3]=Point(2.,0., 1.);
187   Polygon pg4(_vertices=vertices2, _nnodes=3);
188   vertices2[0]=Point(3.,1.,0.);
189   vertices2[1]=Point(1.,4.,0.);
190   vertices2[2]=Point(1.,4.,1.);
191   vertices2[3]=Point(3.,1., 1.);
192   Polygon pg5(_vertices=vertices2, _nnodes=3);
193   vertices2[0]=Point(1.,4.,0.);
194   vertices2[1]=Point(-1.,2.,0.);
195   vertices2[2]=Point(-1.,2.,1.);
196   vertices2[3]=Point(1.,4., 1.);
197   Polygon pg6(_vertices=vertices2, _nnodes=3);
198   vertices2[0]=Point(-1.,2.,0.);
199   vertices2[1]=Point(0.,0.,0.);
200   vertices2[2]=Point(0.,0.,1.);
201   vertices2[3]=Point(-1.,2., 1.);
202   Polygon pg7(_vertices=vertices2, _nnodes=3);
203   std::vector<Polygon> faces(7);
204   faces[0]=pg1;
205   faces[1]=pg2;
206   faces[2]=pg3;
207   faces[3]=pg4;
208   faces[4]=pg5;
209   faces[5]=pg6;
210   faces[6]=pg7;
211   Polyhedron ph(_faces=faces);
212   Mesh meshPolyhedron(ph, _triangle, 1, _gmsh);
213   fin.open("xlifepp_script.geo");
214   it_l=0;
215   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
216   fin.close();
217 
218   out << std::endl << "*********************************************" << std::endl;
219 
220   std::cout << "--- Cylinder with circular basis, with tetrahedron" << std::endl;
221   out << "--- Cylinder with circular basis, with tetrahedron" << std::endl << std::endl;
222   std::vector<Number> ncyl(12,3);
223   Cylinder cyl(_center1=Point(0.,0.,0.), _v1=Point(2.,0.,0.), _v2=Point(0.,2.,0.), _center2=Point(-1.,0.,4.), _nnodes=ncyl,
224                _domain_name="Omega", _side_names="Gamma");
225   Mesh meshCylinder(cyl, _tetrahedron, 1, _gmsh);
226   fin.open("xlifepp_script.geo");
227   it_l=0;
228   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
229   fin.close();
230 
231   out << std::endl << "*********************************************" << std::endl;
232 
233   std::cout << "--- Prism with triangular basis, with tetrahedron" << std::endl;
234   out << "--- Prism with triangular basis, with tetrahedron" << std::endl << std::endl;
235   std::vector<Real> dir(3);
236   dir[0]=-1.;
237   dir[2]=4.;
238   std::vector<Number> nprism1(9,3);
239   Prism prism1(_v1=Point(0.,0.,0.), _v2=Point(3.,0.,0.), _v3=Point(0.,2.,0.), _direction=dir, _nnodes=nprism1, _domain_name="Omega",
240                _side_names="Gamma");
241   Mesh meshPrism1(prism1, _tetrahedron, 1, _gmsh);
242   fin.open("xlifepp_script.geo");
243   it_l=0;
244   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
245   fin.close();
246 
247   out << std::endl << "*********************************************" << std::endl;
248 
249   std::cout << "--- Prism with squared basis, with tetrahedron" << std::endl;
250   out << "--- Prism with squared basis, with tetrahedron" << std::endl << std::endl;
251   std::vector<Number> nprism2(12,3);
252   Prism prism2(_basis=Square(_origin=Point(0.,0.,0.), _length=2.), _direction=dir, _nnodes=nprism2, _domain_name="Omega",
253                _side_names=Strings(6,"Gamma"));
254   Mesh meshPrism2(prism2, _tetrahedron, 1, _gmsh);
255   fin.open("xlifepp_script.geo");
256   it_l=0;
257   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
258   fin.close();
259 
260   out << std::endl << "*********************************************" << std::endl;
261 
262   std::cout << "--- RevCylinder with triangles" << std::endl;
263   out << "--- RevCylinder with triangles" << std::endl << std::endl;
264   RevCylinder revcyl(_center1=Point(0.,0.,0.), _radius=2., _center2=Point(0.,0.,4.), _nnodes=6, _side_names="Gamma");
265   Mesh meshRevCyl(revcyl, _triangle, 1, _gmsh);
266   fin.open("xlifepp_script.geo");
267   it_l=0;
268   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
269   fin.close();
270 
271   out << std::endl << "*********************************************" << std::endl;
272 
273   // a cone with elliptical basis does not work with gmsh engine !!!!!!
274   // std::cout << "--- Cone" << std::endl;
275   // Cone cone(_center1=Point(0.,0.,0.), _v1=Point(2.,0.,0.), _v2=Point(0.,2.,0.), _apex=Point(0.,0.,4.), _nnodes=11,
276   //           _side_names="Gamma");
277   // Mesh meshCone(cone, _tetrahedron, 1, _gmsh);
278   // fin.open("xlifepp_script.geo");
279   // it_l=0;
280   // while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
281   // fin.close();
282 
283   // out << std::endl << "*********************************************" << std::endl;
284 
285   std::cout << "--- RevCone with triangles" << std::endl;
286   out << "--- RevCone with triangles" << std::endl << std::endl;
287   RevCone revcone(_center=Point(0.,0.,0.), _radius=2., _apex=Point(0.,0.,4.), _nnodes=6, _side_names="Gamma");
288   Mesh meshRevCone(revcone, _triangle, 1, _gmsh);
289   fin.open("xlifepp_script.geo");
290   it_l=0;
291   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
292   fin.close();
293 
294   out << std::endl << "*********************************************" << std::endl;
295 
296   Ellipse ell1(_center=Point(0.,0.,0.), _v1=Point(4,0.,0.), _v2=Point(0.,5.,0.), _nnodes=6, _domain_name="Omega1",
297                _side_names=Strings("Gamma_1","Gamma_2","Gamma_3","Gamma_4"));
298   Rectangle rect1(_xmin=-2., _xmax=2., _ymin=-4., _ymax=4., _nnodes=6, _domain_name="Omega2",
299                   _side_names=Strings("Gamma_5","Gamma_6","Gamma_7","Gamma_8"));
300   Ellipse ell2(_center=Point(1.,2.,0.), _v1=Point(1.5,2.,0.), _v2=Point(1.,3.,0.), _nnodes=6, _domain_name="Omega3",
301                _side_names=Strings("Gamma_9","Gamma_10","Gamma_11","Gamma_12"));
302   Ellipse ell3(_center=Point(0.,0.,0.), _v1=Point(0.5,0.,0.), _v2=Point(0.,1.,0.), _nnodes=6, _domain_name="Omega4",
303                _side_names=Strings("Gamma_13","Gamma_14","Gamma_15","Gamma_16"));
304   Rectangle rect2(_xmin=5., _xmax=6., _ymin=0., _ymax=1., _nnodes=6, _domain_name="Omega5",
305                   _side_names=Strings("Gamma_17","Gamma_18","Gamma_19","Gamma_20"));
306   Disk disk1(_center=Point(5.5,0.5,0.), _v1=Point(5.7,0.5,0.), _v2=Point(5.5,0.7,0.), _nnodes=6, _domain_name="Omega6",
307              _side_names=Strings("Gamma_21","Gamma_22","Gamma_23","Gamma_24"));
308 
309   // case canonical - canonical
310   std::cout << "--- Ellipse - Ellipse with triangles" << std::endl;
311   out << "--- Ellipse - Ellipse with triangles" << std::endl << std::endl;
312   Mesh mesh2dP1Composite1(ell1-ell2, _triangle, 1, _gmsh);
313   fin.open("xlifepp_script.geo");
314   it_l=0;
315   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
316   fin.close();
317 
318   out << std::endl << "*********************************************" << std::endl;
319 
320   // case canonical + canonical
321   std::cout << "--- Ellipse + Ellipse with triangles" << std::endl;
322   out << "--- Ellipse + Ellipse with triangles" << std::endl << std::endl;
323   Mesh mesh2dP1Composite2(ell1+ell2, _triangle, 1, _gmsh);
324   fin.open("xlifepp_script.geo");
325   it_l=0;
326   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
327   fin.close();
328 
329   out << std::endl << "*********************************************" << std::endl;
330 
331   // case composite - canonical
332   std::cout << "--- (Ellipse + Rectangle) - Ellipse with triangles" << std::endl;
333   out << "--- (Ellipse + Rectangle) - Ellipse with triangles" << std::endl << std::endl;
334   Mesh mesh2dP1Composite3((ell1+rect1)-ell2, _triangle, 1, _gmsh);
335   fin.open("xlifepp_script.geo");
336   it_l=0;
337   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
338   fin.close();
339 
340   out << std::endl << "*********************************************" << std::endl;
341 
342   // case composite + canonical
343   std::cout << "--- (Ellipse + Rectangle) + Ellipse with triangles" << std::endl;
344   out << "--- (Ellipse + Rectangle) + Ellipse with triangles" << std::endl << std::endl;
345   Mesh mesh2dP1Composite4((ell1+rect1)+ell2, _triangle, 1, _gmsh);
346   fin.open("xlifepp_script.geo");
347   it_l=0;
348   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
349   fin.close();
350 
351   out << std::endl << "*********************************************" << std::endl;
352 
353   // case canonical - composite
354   std::cout << "--- Rectangle - (Ellipse + Ellipse) with triangles" << std::endl;
355   out << "--- Rectangle - (Ellipse + Ellipse) with triangles" << std::endl << std::endl;
356   Mesh mesh2dP1Composite5(rect1-(ell2+ell3), _triangle, 1, _gmsh);
357   fin.open("xlifepp_script.geo");
358   it_l=0;
359   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
360   fin.close();
361 
362   out << std::endl << "*********************************************" << std::endl;
363 
364   // case canonical + composite
365   std::cout << "--- Rectangle + (Ellipse + Ellipse) with triangles" << std::endl;
366   out << "--- Rectangle + (Ellipse + Ellipse) with triangles" << std::endl << std::endl;
367   Mesh mesh2dP1Composite6(rect1+(ell2+ell3), _triangle, 1, _gmsh);
368   fin.open("xlifepp_script.geo");
369   it_l=0;
370   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
371   fin.close();
372 
373   out << std::endl << "*********************************************" << std::endl;
374 
375   // case composite - composite
376   std::cout << "--- (Ellipse + Rectangle) - (Ellipse + Ellipse) with triangles" << std::endl;
377   out << "--- (Ellipse + Rectangle) - (Ellipse + Ellipse) with triangles" << std::endl << std::endl;
378   Mesh mesh2dP1Composite7((ell1+rect1)-(ell2+ell3), _triangle, 1, _gmsh);
379   fin.open("xlifepp_script.geo");
380   it_l=0;
381   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
382   fin.close();
383 
384   out << std::endl << "*********************************************" << std::endl;
385 
386   // case composite + composite
387   std::cout << "--- (Ellipse + Rectangle) + (Ellipse + Ellipse) with triangles" << std::endl;
388   out << "--- (Ellipse + Rectangle) + (Ellipse + Ellipse) with triangles" << std::endl << std::endl;
389   Mesh mesh2dP1Composite8((ell1+rect1)+(ell2+ell3), _triangle, 1, _gmsh);
390   fin.open("xlifepp_script.geo");
391   it_l=0;
392   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
393   fin.close();
394 
395   out << std::endl << "*********************************************" << std::endl;
396 
397   // case multi-composite
398   std::cout << "--- (Ellipse + Rectangle) - (Ellipse + Ellipse) + Rectangle - Disk with triangles" << std::endl;
399   out << "--- (Ellipse + Rectangle) - (Ellipse + Ellipse) + Rectangle - Disk with triangles" << std::endl << std::endl;
400   Mesh mesh2dP1Composite9((ell1+rect1)-(ell2+ell3)+rect2-disk1, _triangle, 1, _gmsh);
401   fin.open("xlifepp_script.geo");
402   it_l=0;
403   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
404   fin.close();
405 
406   out << std::endl << "*********************************************" << std::endl;
407 
408   // case composite 3d
409   std::cout << "--- Ellipsoid - Parallelepiped with tetrahedra" << std::endl;
410   out << "--- Ellipsoid - Parallelepiped with tetrahedra" << std::endl << std::endl;
411   Ellipsoid elld1(_center=Point(0.,0.,0.), _v1=Point(3.,0.,0.), _v2=Point(0.,2.,0.), _v6=Point(0.,0.,1.), _nnodes=5,
412                   _domain_name="Omega");
413   Parallelepiped para1(_v1=Point(-0.5,-0.5,-0.5), _v2=Point(0.5,-0.5,-0.5), _v4=Point(-0.5,0.5,-0.5), _v5=Point(-0.5,-0.5,0.5),
414                        _nnodes=3);
415   Mesh mesh3dP1Composite((elld1-para1), _tetrahedron, 1, _gmsh);
416   fin.open("xlifepp_script.geo");
417   it_l=0;
418   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
419   fin.close();
420 
421   out << std::endl << "*********************************************" << std::endl;
422 
423   Point a1(-1.5,-4.,0.);
424   Point b1(1.5,-4.,0.);
425   Point c1(2.,-3.5,0.);
426   Point d1(2.,3.5,0.);
427   Point e1(1.5,4.,0.);
428   Point f1(-1.5,4.,0.);
429   Point g1(-2.,3.5,0.);
430   Point h1(-2.,-3.5,0.);
431   Segment seg1(_v1=a1, _v2=b1, _nnodes=11, _domain_name="AB");
432   CircArc circ1(_center=Point(1.5,-3.5,0.), _v1=b1, _v2=c1, _nnodes=3, _domain_name="BC");
433   Segment seg2(_v1=c1, _v2=d1, _nnodes=6, _domain_name="CD");
434   CircArc circ2(_center=Point(1.5,3.5,0.), _v1=d1, _v2=e1, _nnodes=3, _domain_name="DE");
435   Segment seg3(_v1=e1, _v2=f1, _nnodes=11, _domain_name="EF");
436   CircArc circ3(_center=Point(-1.5,3.5,0.), _v1=f1, _v2=g1, _nnodes=3, _domain_name="FG");
437   Segment seg4(_v1=g1, _v2=h1, _nnodes=6, _domain_name="GH");
438   CircArc circ4(_center=Point(-1.5,-3.5,0.), _v1=h1, _v2=a1, _nnodes=3, _domain_name="HA");
439 
440   std::cout << "--- Rounded Rectangle (loop) with triangles" << std::endl;
441   out << "--- Rounded Rectangle (loop) with triangles" << std::endl << std::endl;
442   Mesh mesh2dP1Loop(surfaceFrom(seg1+circ1+seg2+circ2+seg3+circ3+seg4+circ4), _triangle, 1, _gmsh);
443   fin.open("xlifepp_script.geo");
444   it_l=0;
445   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
446   fin.close();
447 
448   out << std::endl << "*********************************************" << std::endl;
449 
450   std::cout << "--- Rounded Rectangle with triangles (copy)" << std::endl;
451   out << "--- Rounded Rectangle with triangles (copy)" << std::endl << std::endl;
452   Geometry sf1=surfaceFrom(seg1+circ1+seg2+circ2+seg3+circ3+seg4+circ4, "Omega2");
453   Mesh mesh2dP1Loop2(sf1, _triangle, 1, _gmsh);
454   fin.open("xlifepp_script.geo");
455   it_l=0;
456   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
457   fin.close();
458 
459   out << std::endl << "*********************************************" << std::endl;
460 
461   Point a2(0,0,0);
462   Point b2(2,0,0);
463   Point c2(2,1,0);
464   Point d2(0,1,0);
465   Point e2(0,0,1);
466   Point f2(2,0,1);
467   Point g2(2,1,1);
468   Point h2(0,1,1);
469   Segment seg5(_v1=a2, _v2=b2, _nnodes=6, _domain_name="AB");
470   Segment seg6(_v1=b2, _v2=c2, _nnodes=6, _domain_name="BC");
471   Segment seg7(_v1=c2, _v2=d2, _nnodes=6, _domain_name="CD");
472   Segment seg8(_v1=d2, _v2=a2, _nnodes=6, _domain_name="DA");
473   Rectangle rrect1(_v1=a2, _v2=b2, _v4=d2, _nnodes=6, _domain_name="R1");
474   Rectangle rrect2(_v1=e2, _v2=f2, _v4=h2, _nnodes=6, _domain_name="R2");
475   Rectangle rrect3(_v1=a2, _v2=b2, _v4=e2, _nnodes=6, _domain_name="R3");
476   Rectangle rrect4(_v1=b2, _v2=c2, _v4=f2, _nnodes=6, _domain_name="R4");
477   Rectangle rrect5(_v1=c2, _v2=d2, _v4=g2, _nnodes=6, _domain_name="R5");
478   Rectangle rrect6(_v1=d2, _v2=a2, _v4=h2, _nnodes=6, _domain_name="R6");
479   Geometry sf4=surfaceFrom(seg5+seg6+seg7+seg8, "R1");
480 
481   std::cout << "--- Cube (loop) with tetrahedra" << std::endl;
482   out << "--- Cube (loop) with tetrahedra" << std::endl << std::endl;
483   Mesh mesh3dP1Loop1(volumeFrom(rrect1+rrect2+rrect3+rrect4+rrect5+rrect6), _tetrahedron, 1, _gmsh);
484   fin.open("xlifepp_script.geo");
485   it_l=0;
486   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
487   fin.close();
488 
489   out << std::endl << "*********************************************" << std::endl;
490 
491   std::cout << "--- Cube (loop with border loop) with tetrahedra" << std::endl;
492   out << "--- Cube (loop with border loop) with tetrahedra" << std::endl << std::endl;
493   Mesh mesh3dP1Loop2(volumeFrom(sf4+rrect2+rrect3+rrect4+rrect5+rrect6), _tetrahedron, 1, _gmsh);
494   fin.open("xlifepp_script.geo");
495   it_l=0;
496   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
497   fin.close();
498 
499   out << std::endl << "*********************************************" << std::endl;
500 
501   // half-sphere
502   std::cout << "--- Half-sphere with ruledSurfaceFrom/VolumeFrom" << std::endl;
503   out << "--- Half-sphere with ruledSurfaceFrom/VolumeFrom" << std::endl << std::endl;
504   Real csconeheight=8, csbaseradius=2, cssphereheight=csbaseradius*csbaseradius/csconeheight;
505   Real cssphereradius=std::sqrt(csbaseradius*csbaseradius+cssphereheight*cssphereheight);
506   Number nnodes=6;
507   Point csbasecenter(0.,0.,0.);
508   Point cscenter(0.,0.,-cssphereheight);
509   Point cstop(0.,0.,-cssphereheight-cssphereradius);
510   Point csb1(csbaseradius,0.,0.), csb2(0.,csbaseradius,0.), csb3(-csbaseradius,0.,0.), csb4(0.,-csbaseradius,0.);
511   CircArc csarc1(_center=cscenter, _v1=cstop, _v2=csb1, _nnodes=nnodes);
512   CircArc csarc2(_center=cscenter, _v1=cstop, _v2=csb2, _nnodes=nnodes);
513   CircArc csarc3(_center=cscenter, _v1=cstop, _v2=csb3, _nnodes=nnodes);
514   CircArc csarc4(_center=cscenter, _v1=cstop, _v2=csb4, _nnodes=nnodes);
515   CircArc csarcb1(_center=csbasecenter, _v1=csb1, _v2=csb2, _nnodes=nnodes);
516   CircArc csarcb2(_center=csbasecenter, _v1=csb2, _v2=csb3, _nnodes=nnodes);
517   CircArc csarcb3(_center=csbasecenter, _v1=csb3, _v2=csb4, _nnodes=nnodes);
518   CircArc csarcb4(_center=csbasecenter, _v1=csb4, _v2=csb1, _nnodes=nnodes);
519   Geometry fs1=ruledSurfaceFrom(csarcb1+(~csarc2)+csarc1,"Sphere1");
520   Geometry fs2=ruledSurfaceFrom(csarcb2+(~csarc3)+csarc2,"Sphere2");
521   Geometry fs3=ruledSurfaceFrom(csarcb3+(~csarc4)+csarc3,"Sphere3");
522   Geometry fs4=ruledSurfaceFrom(csarcb4+(~csarc1)+csarc4,"Sphere4");
523   Geometry fs5=planeSurfaceFrom(csarcb1+csarcb2+csarcb3+csarcb4,"Base");
524   Mesh mesh3dP1Halfsphere(volumeFrom(fs1+fs2+fs3+fs4+fs5,"Omega"), _triangle, 1, _gmsh);
525   fin.open("xlifepp_script.geo");
526   it_l=0;
527   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
528   fin.close();
529 
530   out << std::endl << "*********************************************" << std::endl;
531 
532   // false conesphere
533   std::cout << "--- Conesphere with ruledSurfaceFrom/VolumeFrom" << std::endl;
534   out << "--- Conesphere with ruledSurfaceFrom/VolumeFrom" << std::endl << std::endl;
535   Point csapex(0.,0.,csconeheight);
536   Segment csseg1(_v1=csapex, _v2=csb1, _nnodes=3*nnodes);
537   Segment csseg2(_v1=csapex, _v2=csb2, _nnodes=3*nnodes);
538   Segment csseg3(_v1=csapex, _v2=csb3, _nnodes=3*nnodes);
539   Segment csseg4(_v1=csapex, _v2=csb4, _nnodes=3*nnodes);
540   Geometry fc1=ruledSurfaceFrom(csarcb1+(~csseg2)+csseg1,"Cone1");
541   Geometry fc2=ruledSurfaceFrom(csarcb2+(~csseg3)+csseg2,"Cone2");
542   Geometry fc3=ruledSurfaceFrom(csarcb3+(~csseg4)+csseg3,"Cone3");
543   Geometry fc4=ruledSurfaceFrom(csarcb4+(~csseg1)+csseg4,"Cone4");
544 
545   Mesh mesh3dP1Conesphere(volumeFrom(fc1+fc2+fc3+fc4+fs1+fs2+fs3+fs4,"Omega"), _triangle, 1, _gmsh);
546   fin.open("xlifepp_script.geo");
547   it_l=0;
548   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
549   fin.close();
550 
551   out << std::endl << "*********************************************" << std::endl;
552 
553   Point i1(3,0,0);
554   Point j1(3,1,0);
555   Point k1(3,0,1);
556   Point l1(3,1,1);
557   Rectangle rrect7(_v1=b2, _v2=i1, _v4=c2, _nnodes=6, _domain_name="R7");
558 
559   std::cout << "--- Adjacent Rectangles with triangles" << std::endl;
560   out << "--- Adjacent Rectangles with triangles" << std::endl << std::endl;
561   Mesh mesh2dP1CompositeShared(rrect1+rrect7, _triangle, 1, _gmsh);
562   fin.open("xlifepp_script.geo");
563   it_l=0;
564   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
565   fin.close();
566 
567   out << std::endl << "*********************************************" << std::endl;
568 
569   Parallelepiped ppara1(_v1=a2, _v2=b2, _v4=d2, _v5=e2, _nnodes=6, _domain_name="P1");
570   Parallelepiped ppara2(_v1=b2, _v2=i1, _v4=c2, _v5=f2, _nnodes=6, _domain_name="P2");
571 
572   std::cout << "--- Adjacent Parallelepipeds with tetrahedra" << std::endl;
573   out << "--- Adjacent Parallelepipeds with tetrahedra" << std::endl << std::endl;
574   Mesh mesh3dP1CompositeShared(ppara1+ppara2, _triangle, 1, _gmsh);
575 
576   fin.open("xlifepp_script.geo");
577   it_l=0;
578   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
579   fin.close();
580 
581   out << std::endl << "*********************************************" << std::endl;
582 
583   // canonical - loop
584   std::cout << "--- Ellipse - Rounded Rectangle (loop) with triangles" << std::endl;
585   out << "--- Ellipse - Rounded Rectangle (loop) with triangles" << std::endl << std::endl;
586   Geometry g=ell1-sf1;
587   g.domName("Omega"); // BUG is ignored in this case, it does not replace the name of ell1 when analyzed
588   Mesh mesh2dP1CompositeLoopHole1(g, _triangle, 1, _gmsh);
589   fin.open("xlifepp_script.geo");
590   it_l=0;
591   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
592   fin.close();
593 
594   out << std::endl << "*********************************************" << std::endl;
595 
596   // canonical + loop (loop inside canonical)
597   std::cout << "--- Ellipse + Rounded Rectangle (loop) with triangles" << std::endl;
598   out << "--- Ellipse + Rounded Rectangle (loop) with triangles" << std::endl << std::endl;
599   Mesh mesh2dP1CompositeLoopHole2(ell1+sf1, _triangle, 1, _gmsh);
600   fin.open("xlifepp_script.geo");
601   it_l=0;
602   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
603   fin.close();
604 
605   out << std::endl << "*********************************************" << std::endl;
606 
607   // loop - canonical
608   std::cout << "--- Rounded Rectangle (loop) - Ellipse with triangles" << std::endl;
609   out << "--- Rounded Rectangle (loop) - Ellipse with triangles" << std::endl << std::endl;
610   Mesh mesh2dP1CompositeLoopHole3(sf1-ell2, _triangle, 1, _gmsh);
611   fin.open("xlifepp_script.geo");
612   it_l=0;
613   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
614   fin.close();
615 
616   out << std::endl << "*********************************************" << std::endl;
617 
618   // loop + canonical (canonical inside loop : undetermined inclusion -> need to force)
619   std::cout << "--- Rounded Rectangle (loop) + Ellipse with triangles (need to force inclusion)" << std::endl;
620   out << "--- Rounded Rectangle (loop) + Ellipse with triangles (need to force inclusion)" << std::endl << std::endl;
621   Mesh mesh2dP1CompositeLoopHole4(sf1+(+ell2), _triangle, 1, _gmsh);
622   fin.open("xlifepp_script.geo");
623   it_l=0;
624   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
625   fin.close();
626 
627   out << std::endl << "*********************************************" << std::endl;
628 
629   // loop + canonical (loop inside canonical)
630   std::cout << "--- Rounded Rectangle (loop) + Ellipse with triangles" << std::endl;
631   out << "--- Rounded Rectangle (loop) + Ellipse with triangles" << std::endl << std::endl;
632   Mesh mesh2dP1CompositeLoopHole5(sf1+ell1, _triangle, 1, _gmsh);
633   fin.open("xlifepp_script.geo");
634   it_l=0;
635   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
636   fin.close();
637 
638   out << std::endl << "*********************************************" << std::endl;
639 
640   // composite - loop
641   std::cout << "--- (Ellipse + Rectangle) - Rounded Rectangle (loop) with triangles" << std::endl;
642   out << "--- (Ellipse + Rectangle) - Rounded Rectangle (loop) with triangles" << std::endl << std::endl;
643   Mesh mesh2dP1CompositeLoopHole6((ell1+rect2)-sf1, _triangle, 1, _gmsh);
644   fin.open("xlifepp_script.geo");
645   it_l=0;
646   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
647   fin.close();
648 
649   out << std::endl << "*********************************************" << std::endl;
650 
651   // composite + loop (loop inside one component)
652   std::cout << "--- (Ellipse + Rectangle) + Rounded Rectangle (loop) with triangles" << std::endl;
653   out << "--- (Ellipse + Rectangle) + Rounded Rectangle (loop) with triangles" << std::endl << std::endl;
654   Mesh mesh2dP1CompositeLoopHole7((ell1+rect2)+sf1, _triangle, 1, _gmsh);
655   fin.open("xlifepp_script.geo");
656   it_l=0;
657   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
658   fin.close();
659 
660   out << std::endl << "*********************************************" << std::endl;
661 
662   // loop - composite
663   std::cout << "--- Rounded Rectangle (loop) - (Ellipse + Ellipse) with triangles" << std::endl;
664   out << "--- Rounded Rectangle (loop) - (Ellipse + Ellipse) with triangles" << std::endl << std::endl;
665   Mesh mesh2dP1CompositeLoopHole8(sf1-(ell2+ell3), _triangle, 1, _gmsh);
666   fin.open("xlifepp_script.geo");
667   it_l=0;
668   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
669   fin.close();
670 
671   out << std::endl << "*********************************************" << std::endl;
672 
673   // loop + composite (composite inside loop : undetermined inclusions -> need to force)
674   std::cout << "--- Rounded Rectangle (loop) + (Ellipse + Ellipse) with triangles (need to force inclusion)" << std::endl;
675   out << "--- Rounded Rectangle (loop) + (Ellipse + Ellipse) with triangles (need to force inclusion)" << std::endl << std::endl;
676   Mesh mesh2dP1CompositeLoopHole9(sf1+(+(ell2+ell3)), _triangle, 1, _gmsh);
677   fin.open("xlifepp_script.geo");
678   it_l=0;
679   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
680   fin.close();
681 
682   out << std::endl << "*********************************************" << std::endl;
683 
684   // composite with loop - composite
685   std::cout << "--- (Ellipse + Rounded Rectangle (loop)) - (Ellipse + Ellipse) with triangles" << std::endl;
686   out << "--- (Ellipse + Rounded Rectangle (loop) - (Ellipse + Ellipse) with triangles" << std::endl << std::endl;
687   Mesh mesh2dP1CompositeLoopHole10((ell1+sf1)-(ell2+ell3), _triangle, 1, _gmsh);
688   fin.open("xlifepp_script.geo");
689   it_l=0;
690   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
691   fin.close();
692 
693   out << std::endl << "*********************************************" << std::endl;
694 
695   // composite with loop + composite (composite inside loop : undetermined inclusions -> need to rewrite and force)
696   // (ell1+sf1)+(ell2+ell3) should be written as follows
697   // (sf1+(+(ell2+ell3)))+ell1 or ell1+(sf1+(+(ell2+ell3)))
698   std::cout << "--- (Rounded Rectangle (loop) + (Ellipse + Ellipse)) + Ellipse with triangles (need to force inclusion)" << std::endl;
699   out << "--- (Rounded Rectangle (loop) + (Ellipse + Ellipse)) + Ellipse with triangles (need to force inclusion)" << std::endl
700       << std::endl;
701   Mesh mesh2dP1CompositeLoopHole12((sf1+(+(ell2+ell3)))+ell1, _triangle, 1, _gmsh);
702   fin.open("xlifepp_script.geo");
703   it_l=0;
704   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
705   fin.close();
706 
707   out << std::endl << "*********************************************" << std::endl;
708 
709   std::cout << "--- Ellipse + (Rounded Rectangle (loop) + (Ellipse + Ellipse)) with triangles (need to force inclusion)" << std::endl;
710   out << "--- Ellipse + (Rounded Rectangle (loop) + (Ellipse + Ellipse) with triangles (need to force inclusion)" << std::endl
711       << std::endl;
712   Mesh mesh2dP1CompositeLoopHole13(ell1+(sf1+(+(ell2+ell3))), _triangle, 1, _gmsh);
713   fin.open("xlifepp_script.geo");
714   it_l=0;
715   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
716   fin.close();
717 
718   out << std::endl << "*********************************************" << std::endl;
719 
720   // loop - loop
721   Point a3(0,0,0.);
722   Point b3(1,0,0.);
723   Point c3(0,1,0.);
724   Point d3(-1,0,0.);
725   Segment seg9(_v1=d3, _v2=b3, _nnodes=6);
726   CircArc circ5(_center=a3, _v1=b3, _v2=c3, _nnodes=3);
727   CircArc circ6(_center=a3, _v1=c3, _v2=d3, _nnodes=3);
728   Geometry sf2=surfaceFrom(seg9+circ5+circ6,"Omega7");
729 
730   std::cout << "--- Rounded Rectangle (loop) - Half Disk (loop) with triangles" << std::endl;
731   out << "--- Rounded Rectangle (loop) - Half Disk (loop) with triangles" << std::endl << std::endl;
732   Mesh mesh2dP1CompositeLoopHole14(sf1-sf2, _triangle, 1, _gmsh);
733   fin.open("xlifepp_script.geo");
734   it_l=0;
735   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
736   fin.close();
737 
738   out << std::endl << "*********************************************" << std::endl;
739 
740   // loop + loop (undetermined inclusion -> need to force)
741   std::cout << "--- Rounded Rectangle (loop) + Half Disk (loop) with triangles (need to force inclusion)" << std::endl;
742   out << "--- Rounded Rectangle (loop) + Half Disk (loop) with triangles (need to force inclusion)" << std::endl << std::endl;
743   Mesh mesh2dP1CompositeLoopHole15(sf1+(+sf2), _triangle, 1, _gmsh);
744   fin.open("xlifepp_script.geo");
745   it_l=0;
746   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
747   fin.close();
748 
749   out << std::endl << "*********************************************" << std::endl;
750 
751   // loop + loop (no inclusion)
752   Point a4(5.5,0.5,0.);
753   Point b4(5.7,0.5,0.);
754   Point c4(5.5,0.7,0.);
755   Point d4(5.3,0.5,0.);
756   Segment seg10(_v1=d4, _v2=b4, _nnodes=6);
757   CircArc circ7(_center=a4, _v1=b4, _v2=c4, _nnodes=3);
758   CircArc circ8(_center=a4, _v1=c4, _v2=d4, _nnodes=3);
759   Geometry sf3=surfaceFrom(seg10+circ7+circ8, "Omega6");
760 
761   std::cout << "--- Rounded Rectangle (loop) + Half Disk (loop) with triangles" << std::endl;
762   out << "--- Rounded Rectangle (loop) + Half Disk (loop) with triangles" << std::endl << std::endl;
763   Mesh mesh2dP1CompositeLoopHole16(sf1+sf3, _triangle, 1, _gmsh);
764   fin.open("xlifepp_script.geo");
765   it_l=0;
766   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
767   fin.close();
768 
769   out << std::endl << "*********************************************" << std::endl;
770 
771   // composite with loop + loop (no inclusion)
772   std::cout << "--- (Rounded Rectangle (loop) + (Ellipse + Ellipse)) + Half Disk (loop) with triangles" << std::endl;
773   out << "--- (Rounded Rectangle (loop) + (Ellipse + Ellipse)) + Half Disk (loop) with triangles" << std::endl << std::endl;
774   Mesh mesh2dP1CompositeLoopHole17((sf1+(+(ell2+ell3)))+sf3, _triangle, 1, _gmsh);
775   fin.open("xlifepp_script.geo");
776   it_l=0;
777   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
778   fin.close();
779 
780   out << std::endl << "*********************************************" << std::endl;
781 
782   std::cout << "--- (Ellipse + Rounded Rectangle (loop)) - (Ellipse + Ellipse) + Rectangle - Half Disk (loop) with triangles"
783             << std::endl;
784   out << "--- (Ellipse + Rounded Rectangle (loop)) - (Ellipse + Ellipse) + Rectangle - Half Disk (loop) with triangles"
785       << std::endl << std::endl;
786   Mesh mesh2dP1CompositeLoopHole18((ell1+sf1)-(ell2+ell3)+rect2-sf3, _triangle, 1, _gmsh);
787   fin.open("xlifepp_script.geo");
788   it_l=0;
789   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
790   fin.close();
791 
792   out << std::endl << "*********************************************" << std::endl;
793 
794   // composite + loop (loop inside one component + test += operator)
795   g=ell1+rect2;
796   g+=sf1;
797   std::cout << "--- Ellipse + Rectangle += Rounded Rectangle (loop) with triangles" << std::endl;
798   out << "--- Ellipse + Rectangle += Rounded Rectangle (loop) with triangles" << std::endl << std::endl;
799   Mesh mesh2dP1CompositeLoopHole19(g, _triangle, 1, _gmsh);
800   fin.open("xlifepp_script.geo");
801   it_l=0;
802   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
803   fin.close();
804 
805   out << std::endl << "*********************************************" << std::endl;
806 
807   // canonical - loop with use of operator -=
808   // to preserve ell1, we have to convert ell1 int a composite geometry with one component
809   std::cout << "--- Ellipse -= Rounded Rectangle (loop) with triangles" << std::endl;
810   out << "--- Ellipse -= Rounded Rectangle (loop) with triangles" << std::endl << std::endl;
811   g=toComposite(ell1);
812   g-=sf1;
813   Mesh mesh2dP1CompositeLoopHole20(g, _triangle, 1, _gmsh);
814   fin.open("xlifepp_script.geo");
815   it_l=0;
816   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
817   fin.close();
818 
819   out << std::endl << "*********************************************" << std::endl;
820 
821   std::cout << "--- extrusion of Ellipse by Translation with tetrahedra" << std::endl;
822   out << "--- extrusion of Ellipse by Translation with tetrahedra" << std::endl << std::endl;
823   Geometry extr1=extrude(ell1, Translation(0.,0.,4.), 5, "ExtrOmega", "Gamma");
824   Mesh mesh3dP1CompositeExtrusion1(extr1, _tetrahedron, 1, _gmsh);
825   fin.open("xlifepp_script.geo");
826   it_l=0;
827   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
828   fin.close();
829 
830   out << std::endl << "*********************************************" << std::endl;
831 
832   std::cout << "--- extrusion of Rounded Rectangle (loop) by Rotation (pi_/2) with tetrahedra" << std::endl;
833   out << "--- extrusion of Rounded Rectangle (loop) by Rotation (pi_/2) with tetrahedra" << std::endl << std::endl;
834   Geometry extr2=extrude(sf1, Rotation3d(Point(5.,0.,0.),0.,5.,0.,pi_/2), 10, "ExtrOmega2");
835   Mesh mesh3dP1CompositeExtrusion2(extr2, _tetrahedron, 1, _gmsh);
836   fin.open("xlifepp_script.geo");
837   it_l=0;
838   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
839   fin.close();
840 
841   out << std::endl << "*********************************************" << std::endl;
842 
843   std::cout << "--- extrusion of Ellipse - Ellipse by Rotation (pi_/2) with tetrahedra" << std::endl;
844   out << "--- extrusion of Ellipse - Ellipse by Rotation (pi_/2) with tetrahedra" << std::endl << std::endl;
845   Geometry extr3=extrude(ell1-ell2, Rotation3d(Point(5.,0.,0.),0.,5.,0.,pi_/2), 10, "ExtrOmega3", "Gamma");
846   Mesh mesh3dP1CompositeExtrusion3(extr3, _tetrahedron, 1, _gmsh);
847   fin.open("xlifepp_script.geo");
848   it_l=0;
849   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
850   fin.close();
851 
852   out << std::endl << "*********************************************" << std::endl;
853 
854   std::cout << "--- extrusion of Rounded Rectangle (loop) by Rotation (pi_) with tetrahedra" << std::endl;
855   out << "--- extrusion of Rounded Rectangle (loop) by Rotation (pi_) with tetrahedra" << std::endl << std::endl;
856   Geometry extr4=extrude(sf1, Rotation3d(Point(5.,0.,0.),0.,5.,0.,pi_), 10, "ExtrOmega4");
857   Mesh mesh3dP1CompositeExtrusion4(extr4, _tetrahedron, 1, _gmsh);
858   fin.open("xlifepp_script.geo");
859   it_l=0;
860   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
861   fin.close();
862 
863   out << std::endl << "*********************************************" << std::endl;
864 
865   std::cout << "--- extrusion of Ellipse - Ellipse by Rotation (pi_) with tetrahedra" << std::endl;
866   out << "--- extrusion of Ellipse - Ellipse by Rotation (pi_) with tetrahedra" << std::endl << std::endl;
867   Geometry extr5=extrude(ell1-ell2, Rotation3d(Point(5.,0.,0.),0.,5.,0.,pi_), 10, "ExtrOmega5", "Gamma");
868   Mesh mesh3dP1CompositeExtrusion5(extr5, _tetrahedron, 1, _gmsh);
869   fin.open("xlifepp_script.geo");
870   it_l=0;
871   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
872   fin.close();
873 
874   out << std::endl << "*********************************************" << std::endl;
875 
876   std::cout << "--- extrusion of Rounded Rectangle (loop) by Rotation (3*pi_/2) with tetrahedra" << std::endl;
877   out << "--- extrusion of Rounded Rectangle (loop) by Rotation (3*pi_/2) with tetrahedra" << std::endl << std::endl;
878   Geometry extr6=extrude(sf1, Rotation3d(Point(5.,0.,0.),0.,5.,0.,3.*pi_/2), 10, "ExtrOmega6");
879   Mesh mesh3dP1CompositeExtrusion6(extr6, _tetrahedron, 1, _gmsh);
880   fin.open("xlifepp_script.geo");
881   it_l=0;
882   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
883   fin.close();
884 
885   out << std::endl << "*********************************************" << std::endl;
886 
887   std::cout << "--- extrusion of Ellipse - Ellipse by Rotation (3*pi_/2) with tetrahedra" << std::endl;
888   out << "--- extrusion of Ellipse - Ellipse by Rotation (3*pi_/2) with tetrahedra" << std::endl << std::endl;
889   Geometry extr7=extrude(ell1-ell2, Rotation3d(Point(5.,0.,0.),0.,5.,0.,3.*pi_/2), 10, "ExtrOmega7", "Gamma");
890   Mesh mesh3dP1CompositeExtrusion7(extr7, _tetrahedron, 1, _gmsh);
891   fin.open("xlifepp_script.geo");
892   it_l=0;
893   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
894   fin.close();
895 
896   out << std::endl << "*********************************************" << std::endl;
897 
898   std::cout << "--- extrusion of Rounded Rectangle (loop) by Rotation (2*pi_) with tetrahedra" << std::endl;
899   out << "--- extrusion of Rounded Rectangle (loop) by Rotation (2*pi_) with tetrahedra" << std::endl << std::endl;
900   Geometry extr8=extrude(sf1, Rotation3d(Point(5.,0.,0.),0.,5.,0.,2.*pi_), 10, "ExtrOmega8");
901   Mesh mesh3dP1CompositeExtrusion8(extr8, _tetrahedron, 1, _gmsh);
902   fin.open("xlifepp_script.geo");
903   it_l=0;
904   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
905   fin.close();
906 
907   out << std::endl << "*********************************************" << std::endl;
908 
909   std::cout << "--- extrusion of Ellipse - Ellipse by Rotation (2*pi_) with tetrahedra" << std::endl;
910   out << "--- extrusion of Ellipse - Ellipse by Rotation (2*pi_) with tetrahedra" << std::endl << std::endl;
911   Geometry extr9=extrude(ell1-ell2, Rotation3d(Point(5.,0.,0.),0.,5.,0.,2.*pi_), 10, "ExtrOmega9", "Gamma");
912   Mesh mesh3dP1CompositeExtrusion9(extr9, _tetrahedron, 1, _gmsh);
913   fin.open("xlifepp_script.geo");
914   it_l=0;
915   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
916   fin.close();
917 
918   out << std::endl << "*********************************************" << std::endl;
919 
920   std::cout << "--- extrusion of CircArc by Translation with triangles" << std::endl;
921   out << "--- extrusion of CircArc by Translation with triangles" << std::endl << std::endl;
922   Geometry extr10=extrude(circ1, Translation(0.,0.,4.), 5, "ExtrOmega10");
923   Mesh mesh2dP1CompositeExtrusion10(extr10, _triangle, 1, _gmsh);
924   fin.open("xlifepp_script.geo");
925   it_l=0;
926   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
927   fin.close();
928 
929   out << std::endl << "*********************************************" << std::endl;
930 
931   // true conesphere
932   std::cout << "--- Conesphere with extrusion by rotation (2.*pi_)" << std::endl;
933   out << "--- Conesphere with extrusion by rotation (2.*pi_)" << std::endl << std::endl;
934   Real csbaseradiusint=2., csconeheightint=8.;
935   Real cssphereheightint=csbaseradiusint*csbaseradiusint/csconeheightint;
936   Real cssphereradiusint=sqrt(csbaseradiusint*csbaseradiusint+cssphereheightint*cssphereheightint);
937   Real diff=2.;
938   Real csbaseradiusext=2.+diff, csconeheightext=8.+2.*diff;
939   Real cssphereheightext=csbaseradiusext*csbaseradiusext/csconeheightext;
940   Real cssphereradiusext=sqrt(csbaseradiusext*csbaseradiusext+cssphereheightext*cssphereheightext);
941   nnodes=10;
942 
943   Point csapexint(0.,0.,csconeheightint), csp1int(csbaseradiusint,0.,0.), csp2int(0.,0.,-cssphereheightint-cssphereradiusint);
944   Point csapexext(0.,0.,csconeheightext), csp1ext(csbaseradiusext,0.,0.), csp2ext(0.,0.,-cssphereheightext-cssphereradiusext);
945 
946   Segment css1(_v1=csp1ext, _v2=csapexext, _nnodes=2*nnodes);
947   Segment css2(_v1=csapexext, _v2=csapexint, _nnodes=2*nnodes);
948   Segment css3(_v1=csapexint,_v2=csp1int, _nnodes=2*nnodes);
949   CircArc csc1(_center=Point(0.,0.,-cssphereheightint), _v1=csp1int, _v2=csp2int, _nnodes=nnodes);
950   Segment css4(_v1=csp2int, _v2=csp2ext, _nnodes=2*nnodes);
951   CircArc csc2(_center=Point(0.,0.,-cssphereheightext),_v1=csp2ext,_v2=csp1ext, _nnodes=nnodes);
952 
953   Geometry csbase=planeSurfaceFrom(css1+css2+css3+csc1+css4+csc2, "Gamma");
954   Geometry csextrude=extrude(csbase, Rotation3d(Point(0.,0.,0.),0.,0.,1.,2.*pi_), "Omega1",
955                              Strings("Gamma1", "Gamma2", "Gamma3", "Gamma4", "Gamma5", "Gamma6", "Gamma7", "Gamma8",
956                                      "Gamma9", "Gamma10", "Gamma11", "Gamma12", "Gamma13", "Gamma14", "Gamma15", "Gamma16"));
957 
958   Mesh m1(csextrude,_tetrahedron,1,_gmsh);
959   fin.open("xlifepp_script.geo");
960   it_l=0;
961   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
962   fin.close();
963 
964   out << std::endl << "*********************************************" << std::endl;
965 
966   // crack 2D 1st version
967   Point x1(0,0,0);
968   Point x2(1,0,0);
969   Point x3(1,1,0);
970   Point x4(0,1,0);
971   Point x5(0.2,0.2,0);
972   Point x6(0.8,0.8,0);
973   Point x7(0.2,0,0);
974   Point x8(0.8,1,0);
975   Segment ss1(_v1=x1, _v2=x7, _domain_name="Gamma");
976   Segment ss2(_v1=x2, _v2=x7, _domain_name="Gamma");
977   Segment ss3(_v1=x3, _v2=x2, _domain_name="Gamma");
978   Segment ss4(_v1=x8, _v2=x3, _domain_name="Gamma");
979   Segment ss5(_v1=x8, _v2=x4, _domain_name="Gamma");
980   Segment ss6(_v1=x4, _v2=x1, _domain_name="Gamma");
981   Segment ss7(_v1=x7, _v2=x5);
982   Segment ss8(_v1=x5, _v2=x6, _nnodes=3, _domain_name="Crack");
983   Segment ss9(_v1=x6, _v2=x8);
984   crack(ss8);
985   Geometry gsf1=surfaceFrom(ss7+ss8+ss9+ss5+ss6+ss1, "Omega1");
986   Geometry gsf2=surfaceFrom(ss7+ss8+ss9+ss4+ss3+ss2, "Omega2");
987 
988   std::cout << "--- cracked Segment (loop version) with triangles" << std::endl;
989   out << "--- cracked Segment (loop version) with triangles" << std::endl << std::endl;
990   Mesh mesh2dP1CompositeCrack(gsf1+gsf2, _triangle, 1, _gmsh);
991   fin.open("xlifepp_script.geo");
992   it_l=0;
993   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
994   fin.close();
995 
996   out << std::endl << "*********************************************" << std::endl;
997 
998   // crack 2D 2st version
999   Rectangle rrect8(_v1=x1, _v2=x2, _v4=x4, _nnodes=2, _domain_name="Omega", _side_names="Gamma");
1000   Segment scrack(_v1=x5, _v2=x6, _nnodes=3, _domain_name="Crack", _side_names="Sigma");
1001   openCrack(scrack,String("Sigma"));
1002 
1003   std::cout << "--- open cracked Segment (operator+ version) with triangles" << std::endl;
1004   out << "--- open cracked Segment (operator+ version) with triangles" << std::endl << std::endl;
1005   Mesh mesh2dP1CompositeCrack2(rrect8+scrack, _triangle, 1, _gmsh);
1006   fin.open("xlifepp_script.geo");
1007   it_l=0;
1008   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
1009   fin.close();
1010 
1011   out << std::endl << "*********************************************" << std::endl;
1012 
1013   // crack 3D
1014   std::cout << "--- cracked Rectangle with tetrahedra" << std::endl;
1015   out << "--- cracked Rectangle with tetrahedra" << std::endl << std::endl;
1016   Rectangle rcrack(_v1=x1, _v2=x2, _v4=x4, _nnodes=2, _domain_name="Crack");
1017   crack(rcrack);
1018   Mesh mesh3dP1CompositeCrack1(elld1+rcrack, _tetrahedron, 1, _gmsh);
1019   fin.open("xlifepp_script.geo");
1020   it_l=0;
1021   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
1022   fin.close();
1023 
1024   out << std::endl << "*********************************************" << std::endl;
1025 
1026   // crack 3D
1027   std::cout << "--- cracked Segment with triangles (Arnaud's mesh)" << std::endl;
1028   out << "--- cracked Segment with triangles (Arnaud's mesh)" << std::endl << std::endl;
1029 
1030   Point p1ar(-0.3,0.), p2ar(-0.1,0.), p3ar(0.,0.), p4ar(0.3,0.);
1031   Point p5ar(-0.3,0.49), p6ar(-0.1,0.49), p7ar(0.,0.49), p8ar(0.3,0.49);
1032   Point p9ar(-0.3,0.51), p10ar(-0.1,0.51), p11ar(0.,0.51), p12ar(0.3,0.51);
1033   Point p13ar(-0.3,1.), p14ar(-0.1,1.), p15ar(0.,1.), p16ar(0.3,1.);
1034 
1035   Real h = 1., epColle = 0.02;
1036   Number MeshSize = 20;
1037   Segment s1ar(_v1=p1ar, _v2=p2ar, _nnodes=(Number)floor(0.2*MeshSize+1e-5)+1, _domain_name="GAMMAm");
1038   Segment s2ar(_v1=p2ar, _v2=p3ar, _nnodes=(Number)floor(0.1*MeshSize+1e-5)+1, _domain_name="GAMMAm");
1039   Segment s3ar(_v1=p3ar, _v2=p4ar, _nnodes=(Number)floor(0.3*MeshSize+1e-5)+1, _domain_name="GAMMAm");
1040   Segment s4ar(_v1=p5ar, _v2=p6ar, _nnodes=(Number)floor(0.2*MeshSize+1e-5)+1);
1041   Segment s5ar(_v1=p6ar, _v2=p7ar, _nnodes=(Number)floor(0.1*MeshSize+1e-5)+1);
1042   Segment s6ar(_v1=p7ar, _v2=p8ar, _nnodes=(Number)floor(0.3*MeshSize+1e-5)+1);
1043   Segment s7ar(_v1=p9ar, _v2=p10ar, _nnodes=(Number)floor(0.2*MeshSize+1e-5)+1);
1044   Segment s8ar(_v1=p10ar, _v2=p11ar, _nnodes=(Number)floor(0.1*MeshSize+1e-5)+1, _domain_name="DEFAUT"); crack(s8ar);
1045   Segment s9ar(_v1=p11ar, _v2=p12ar, _nnodes=(Number)floor(0.3*MeshSize+1e-5)+1);
1046   Segment s10ar(_v1=p13ar, _v2=p14ar, _nnodes=(Number)floor(0.2*MeshSize+1e-5)+1, _domain_name="GAMMAp");
1047   Segment s11ar(_v1=p14ar, _v2=p15ar, _nnodes=(Number)floor(0.1*MeshSize+1e-5)+1, _domain_name="GAMMAp");
1048   Segment s12ar(_v1=p15ar, _v2=p16ar, _nnodes=(Number)floor(0.3*MeshSize+1e-5)+1, _domain_name="GAMMAp");
1049 
1050   Segment s13ar(_v1=p1ar, _v2=p5ar, _nnodes=(Number)floor((h-epColle)/2*MeshSize+1e-5)+1, _domain_name="SIGMAm");
1051   Segment s14ar(_v1=p5ar, _v2=p9ar, _nnodes=(Number)floor(epColle*MeshSize+1e-5)+1, _domain_name="SIGMAm");
1052   Segment s15ar(_v1=p9ar, _v2=p13ar, _nnodes=(Number)floor((h-epColle)/2*MeshSize+1e-5)+1, _domain_name="SIGMAm");
1053   Segment s16ar(_v1=p2ar, _v2=p6ar, _nnodes=(Number)floor((h-epColle)/2*MeshSize+1e-5)+1);
1054   Segment s17ar(_v1=p6ar, _v2=p10ar, _nnodes=(Number)floor(epColle*MeshSize+1e-5)+1);
1055   Segment s18ar(_v1=p10ar, _v2=p14ar, _nnodes=(Number)floor((h-epColle)/2*MeshSize+1e-5)+1);
1056   Segment s19ar(_v1=p3ar, _v2=p7ar, _nnodes=(Number)floor((h-epColle)/2*MeshSize+1e-5)+1);
1057   Segment s20ar(_v1=p7ar, _v2=p11ar, _nnodes=(Number)floor(epColle*MeshSize+1e-5)+1);
1058   Segment s21ar(_v1=p11ar, _v2=p15ar, _nnodes=(Number)floor((h-epColle)/2*MeshSize+1e-5)+1);
1059   Segment s22ar(_v1=p4ar, _v2=p8ar, _nnodes=(Number)floor((h-epColle)/2*MeshSize+1e-5)+1, _domain_name="SIGMAp");
1060   Segment s23ar(_v1=p8ar, _v2=p12ar, _nnodes=(Number)floor(epColle*MeshSize+1e-5)+1, _domain_name="SIGMAp");
1061   Segment s24ar(_v1=p12ar, _v2=p16ar, _nnodes=(Number)floor((h-epColle)/2*MeshSize+1e-5)+1, _domain_name="SIGMAp");
1062 
1063   Geometry g1ar=surfaceFrom(s1ar+s16ar+(~s4ar)+(~s13ar), "OMEGA1");
1064   Geometry g2ar=surfaceFrom(s2ar+s19ar+(~s5ar)+(~s16ar), "OMEGA1");
1065   Geometry g3ar=surfaceFrom(s3ar+s22ar+(~s6ar)+(~s19ar), "OMEGA1");
1066   Geometry g4ar=surfaceFrom(s4ar+s17ar+(~s7ar)+(~s14ar), "OMEGA2");
1067   Geometry g5ar=surfaceFrom(s5ar+s20ar+(~s8ar)+(~s17ar), "OMEGA2");
1068   Geometry g6ar=surfaceFrom(s6ar+s23ar+(~s9ar)+(~s20ar), "OMEGA2");
1069   Geometry g7ar=surfaceFrom(s7ar+s18ar+(~s10ar)+(~s15ar), "OMEGA1");
1070   Geometry g8ar=surfaceFrom(s8ar+s21ar+(~s11ar)+(~s18ar), "OMEGA1");
1071   Geometry g9ar=surfaceFrom(s9ar+s24ar+(~s12ar)+(~s21ar), "OMEGA1");
1072   Geometry gar=g1ar+g2ar+g3ar+g4ar+g5ar+g6ar+g7ar+g8ar+g9ar;
1073   Mesh mesh2dar(gar,_triangle);
1074   out << "// domains:";
1075   for (number_t i=0; i < mesh2dar.nbOfDomains(); ++i)
1076   {
1077     out << " " << mesh2dar.domainName(i);
1078   }
1079   out << std::endl;
1080   fin.open("xlifepp_script.geo");
1081   it_l=0;
1082   while (fin.getline(s,200)) { if (it_l > 1) { out << s << std::endl; } ++it_l; }
1083   fin.close();
1084 
1085   //------------------------------------------------------------------------------------
1086   //   save results in a file or compare results with some references value in a file
1087   //------------------------------------------------------------------------------------
1088   trace_p->pop();
1089   if (check) { return diffResFile(out, rootname); }
1090   else { return saveResToFile(out, rootname); }
1091 #else
1092   return String(rootname+" skipped");
1093 #endif
1094 }
1095 
1096 }
1097