1 // ********** DO NOT REMOVE THIS BANNER **********
2 //
3 // SUMMARY: Bamg: Bidimensional Anisotrope Mesh Generator
4 // RELEASE: 0
5 //
6 // AUTHOR: F. Hecht,
7 // ORG : UMPC
8 // E-MAIL : Frederic.Hecht@Inria.fr
9 //
10 // ORIG-DATE: Dec 97
11 // Modif March 98
12 /*
13
14 This file is part of Bamg or FreeFEm++
15
16 Freefem++ is free software; you can redistribute it and/or modify
17 it under the terms of the GNU Lesser General Public License as published by
18 the Free Software Foundation; either version 2.1 of the License, or
19 (at your option) any later version.
20
21 Freefem++ is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU Lesser General Public License for more details.
25
26 You should have received a copy of the GNU Lesser General Public License
27 along with Freefem++; if not, write to the Free Software
28 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
29 */
30 // E-MAIL : Frederic.Hecht@Inria.fr
31 //
32 // ORIG-DATE: Dec 97
33 #include <unistd.h>
34 #include <cstdlib>
35 #include <cstdio>
36 #include <cstring>
37 #include <setjmp.h>
38 #include <new>
39 #include <cassert>
40 #include "Meshio.h"
41 #include <iomanip>
42 #include "Mesh2.h"
43 #include "QuadTree.h"
44 using namespace bamg;
45 using namespace std;
46 #include <fstream>
47 #ifdef __MWERKS__
48 #define NBVMAX 10000
49 #else
50 #define NBVMAX 50000
51 #endif
52 // 0 1 2 3
53 // 0123456789012345678901234567890
54 const char whatbamg[] = "@(#)LJLL, feb. 20, 2006, Bamg v1.02";
55 const char *bamgversion = whatbamg + 28;
56 int initgraph = 0;
57 // long verbosity=2;
58
59 #ifdef __MWERKS__
60 #define R_OK 0
61 #define F_OK 0
62 #define W_OK 0
63 #include <unix.h>
access(char * fileName,int notUsed)64 inline int access(char *fileName, int notUsed) {
65 struct stat statRec;
66 return stat(fileName, &statRec);
67 }
68 #endif
69 #ifdef _MSC_VER
70 #define R_OK 04
71 #define F_OK 0
72 #define W_OK 02
73 #include <io.h>
access(char * fileName,int mode)74 inline int access(char *fileName, int mode) { return _access(fileName, mode); }
75 #endif
76
77 void NewHandler( );
NewHandler()78 void NewHandler( ) {
79 cerr << " Not enought memory" << endl;
80 exit(1);
81 }
82 void forDebug( );
83
MeshErrorIO(ios &)84 void MeshErrorIO(ios &) {
85 MeshError(999);
86 exit(1);
87 }
88
89 #ifdef DRAWING
90 bool withrgraphique = initgraphique;
91 #else
92 bool withrgraphique = false;
93 #endif
94
forDebug()95 void forDebug( ) {
96 #ifdef DRAWING
97 if (initgraph) {
98 if (CurrentTh) CurrentTh->inquire( );
99 closegraphique( );
100 initgraph = 0;
101 }
102 #endif
103 }
104
main(int argc,char ** argv)105 int main(int argc, char **argv) {
106 #ifdef __MWERKS__
107 cout << " version " << bamgversion << " for initialisation of the toolbox " << endl;
108 #endif
109 // JUST for HP CC bug
110 #ifdef LONGLONG
111 long long bidon = 2;
112 if ((double)bidon != 2.0) {
113 cerr << "Compiler bug Pb of cast of long long to a double => remove the LONGLONG define"
114 << endl;
115 exit(10);
116 }
117 #endif
118 double time0, time2, time3, time4;
119 // 2 way for uses ---
120 // 1 first mesh
121 // or adaptation
122 MeshIstreamErrorHandler = MeshErrorIO;
123 // double *sol=0;
124 // int adaptation = 0 ;
125 Int4 i;
126 hinterpole = 1; // by def interpolation a h
127 // int Err=0;
128 int fileout = 0;
129 int nbvx = NBVMAX;
130 int iso = 0, AbsError = 0, nbjacoby = 1, allquad = 0;
131 int NoMeshReconstruction = 0;
132 int Rescaling = 1;
133 double costheta = 2;
134 Real8 cutoffradian = -1;
135 double anisomax = 1e6;
136 double err = 0.01, errg = 0.1, coef = 1, hmin = 1.e-100, hmax = 1.e17, raison = 0, CutOff = 1e-5;
137 int KeepBackVertices = 1;
138 double hminaniso = 1e-100;
139 const double boundmaxsubdiv = 10;
140 double maxsubdiv = boundmaxsubdiv;
141 double omega = 1.8;
142 int NbSmooth = 3;
143 double *solMbb = 0, *solMBB = 0;
144 int *typesolsBB = 0;
145 Int4 nbsolbb = 0, lsolbb = 0;
146 Int4 nbsolBB = 0, lsolBB = 0;
147 int SplitEdgeWith2Boundary = 0;
148 int rbbeqMbb = 0, rBBeqMBB = 0;
149 int ChoiseHessien = 0;
150 double power = 1;
151 Triangles *Thr = 0, *Thb = 0;
152
153 char *fgeom = 0, *fmeshback = 0, *fmeshout = 0, *fmeshr = 0, *fmetrix = 0, *famfmt = 0, *fmsh = 0,
154 *fnopo = 0, *fftq = 0, *fam = 0, *famdba = 0, *rbb = 0, *rBB = 0, *wbb = 0, *wBB = 0,
155 *fMbb = 0, *foM = 0, *fMBB = 0;
156 verbosity = 2;
157 char *datargv[128];
158 int datargc = 1;
159 datargv[0] = datargv[1] = 0; // for create a error if no parameter
160 const char *datafile = argc == 2 ? argv[1] : "DATA_bamg";
161
162 atexit(forDebug);
163 set_new_handler(&NewHandler);
164
165 if (argc < 3) { // for test on the mac
166 if (!access(datafile, R_OK)) {
167 MeshIstream data(datafile);
168 while (data.cm( ) && datargc < 128) datargv[datargc++] = data.ReadStr( );
169
170 if (argc == 1) datargv[0] = argv[0]; // copy of the program name
171 argc = datargc;
172 argv = datargv;
173 }
174 };
175
176 for (i = 1; i < argc; i++)
177 if (!strcmp(argv[i], "-v") && ++i < argc)
178 verbosity = atoi(argv[i]);
179 else if (!strcmp(argv[i], "-NoRescaling"))
180 Rescaling = 0;
181 else if (!strcmp(argv[i], "-Rescaling"))
182 Rescaling = 1;
183 else if (!strcmp(argv[i], "-H") && ++i < argc)
184 ChoiseHessien = atoi(argv[i]);
185 else if (!strcmp(argv[i], "-power") && ++i < argc)
186 power = atof(argv[i]);
187 else if (!strcmp(argv[i], "-nbv") && ++i < argc)
188 nbvx = atoi(argv[i]);
189 else if (!strcmp(argv[i], "-nbs") && ++i < argc)
190 nbvx = atoi(argv[i]);
191 else if (!strcmp(argv[i], "-g") && ++i < argc)
192 fgeom = argv[i];
193 else if ((!strcmp(argv[i], "-thetaquad") || (!strcmp(argv[i], "-ThetaQuad"))) && ++i < argc)
194 costheta = 1 - Abs(cos(Pi * atof(argv[i]) / 180.0));
195 else if (!strcmp(argv[i], "-b") && ++i < argc)
196 fmeshback = argv[i];
197 else if (!strcmp(argv[i], "-r") && ++i < argc)
198 fmeshr = argv[i];
199 else if (!strcmp(argv[i], "-M") && ++i < argc)
200 fmetrix = argv[i];
201 else if (!strcmp(argv[i], "-rbb") && ++i < argc)
202 rbb = argv[i];
203 else if (!strcmp(argv[i], "-rBB") && ++i < argc)
204 rBB = argv[i];
205 else if (!strcmp(argv[i], "-Mbb") && ++i < argc)
206 fMbb = argv[i];
207 else if (!strcmp(argv[i], "-MBB") && ++i < argc)
208 fMBB = argv[i];
209 else if (!strcmp(argv[i], "-wBB") && ++i < argc)
210 wBB = argv[i];
211 else if (!strcmp(argv[i], "-wbb") && ++i < argc)
212 wbb = argv[i];
213 else if (!strcmp(argv[i], "-o") && ++i < argc)
214 fmeshout = argv[i];
215 else if (!strcmp(argv[i], "-intm") && i < argc)
216 hinterpole = 0;
217 else if (!strcmp(argv[i], "-oam_fmt") && ++i < argc)
218 famfmt = argv[i];
219 else if (!strcmp(argv[i], "-oam") && ++i < argc)
220 fam = argv[i];
221 else if (!strcmp(argv[i], "-omsh") && ++i < argc)
222 fmsh = argv[i];
223 else if (!strcmp(argv[i], "-oftq") && ++i < argc)
224 fftq = argv[i];
225 else if (!strcmp(argv[i], "-onopo") && ++i < argc)
226 fnopo = argv[i];
227 else if (!strcmp(argv[i], "-oamdba") && ++i < argc)
228 famdba = argv[i];
229 else if (!strcmp(argv[i], "-oM") && ++i < argc)
230 foM = argv[i];
231 else if (!strcmp(argv[i], "-coef") && ++i < argc)
232 coef = atof(argv[i]);
233 else if (!strcmp(argv[i], "-err") && ++i < argc)
234 err = atof(argv[i]);
235 else if (!strcmp(argv[i], "-errg") && ++i < argc)
236 errg = atof(argv[i]);
237 else if (!strcmp(argv[i], "-raison") && ++i < argc)
238 raison = atof(argv[i]);
239 else if (!strcmp(argv[i], "-ratio") && ++i < argc)
240 raison = atof(argv[i]);
241 else if (!strcmp(argv[i], "-hmin") && ++i < argc)
242 hmin = atof(argv[i]);
243 else if (!strcmp(argv[i], "-hminaniso") && ++i < argc)
244 hminaniso = atof(argv[i]);
245 else if (!strcmp(argv[i], "-hmax") && ++i < argc)
246 hmax = atof(argv[i]);
247 else if (!strcmp(argv[i], "-thetamax") && ++i < argc)
248 cutoffradian = atof(argv[i]) * Pi / 180.0;
249 else if (!strcmp(argv[i], "-anisomax") && ++i < argc)
250 anisomax = atof(argv[i]);
251 else if (!strcmp(argv[i], "-maxsubdiv") && ++i < argc)
252 maxsubdiv = atof(argv[i]);
253 else if (!strcmp(argv[i], "-iso"))
254 iso = 1;
255 else if (!strcmp(argv[i], "-KeepBackVertices"))
256 KeepBackVertices = 1;
257 else if (!strcmp(argv[i], "-noKeepBackVertices"))
258 KeepBackVertices = 0;
259 else if (!strcmp(argv[i], "-2q") || !strcmp(argv[i], "-2Q"))
260 allquad = 1;
261 else if (!strcmp(argv[i], "-2"))
262 allquad = 2;
263 else if (!strcmp(argv[i], "-aniso"))
264 iso = 0;
265 else if (!strcmp(argv[i], "-RelError"))
266 AbsError = 0;
267 else if (!strcmp(argv[i], "-AbsError"))
268 AbsError = 1;
269 else if (!strcmp(argv[i], "-splitpbedge"))
270 SplitEdgeWith2Boundary = 1;
271 else if (!strcmp(argv[i], "-nosplitpbedge"))
272 SplitEdgeWith2Boundary = 0;
273 else if (!strcmp(argv[i], "-NbJacobi") && ++i < argc)
274 nbjacoby = atoi(argv[i]);
275 else if (!strcmp(argv[i], "-CutOff") && ++i < argc)
276 CutOff = atof(argv[i]), AbsError = 0;
277 else if (!strcmp(argv[i], "-NbSmooth") && ++i < argc)
278 NbSmooth = atoi(argv[i]);
279 else if (!strcmp(argv[i], "-omega") && ++i < argc)
280 omega = atof(argv[i]);
281 else {
282 cout << " Error in the arguments " << i << " of " << bamgversion << endl;
283 cout << " argument [" << i << "]= `" << ((i < argc) ? argv[i] : " The Lastest ") << "`"
284 << endl;
285 cout << " Usage :" << endl;
286 cout << " Mesh INPUT: The 2 arguments are exclusives." << endl;
287 cout << "" << endl;
288 cout << " -g filename Set the geometry for mesh generation. " << endl;
289 cout << " DB mesh file." << endl;
290 cout << " -b filename Set the background mesh for mesh adaption " << endl;
291 cout << " (require {\tt -M} or {-Mbb} arguments). " << endl;
292 cout << " - am_fmt file if filename match *.am_fmt " << endl;
293 cout << " - ambda file if filename match *.amdba " << endl;
294 cout << " - nopo file if filename match *.nopo " << endl;
295 cout << " - am file if filename match *.am " << endl;
296 cout << " - ftq file if filename match *.ftq " << endl;
297 cout << " - msh file if filename match *.msh " << endl;
298 cout << " - otherwise the file is a BD mesh file " << endl;
299 cout << " Remark: the geometry is the background geometry." << endl;
300 cout << " DB mesh file." << endl;
301 cout << " -r filename Set the mesh for modification mesh " << endl;
302 cout << " no reconstruction " << endl;
303 cout << " same as in the case of -b argument" << endl;
304 cout << "" << endl;
305 cout << " -thetamax (double) change the angular limit for a corner in degre " << endl;
306 cout << " the angle is defined from 2 normals of 2 concecutives edges "
307 << endl;
308 cout << " if no geomtry cf reading an am_fmt, amdba, nopo .. file "
309 << endl;
310
311 cout << "" << endl;
312 cout << " METRIC definition or mesh size definition, one of the" << endl;
313 cout << " 2 next arguments is need in case of mesh adaption." << endl;
314 cout << "" << endl;
315 cout << " -M filename Set the metric which is defined on the background mesh "
316 << endl;
317 cout << " or on the geometry. Metric file." << endl;
318 cout << " -Mbb filename Set the solution defined on the background mesh for" << endl;
319 cout << " metric construction, the solutions was FE P1 defined on "
320 << endl;
321 cout << " the background mesh. bb file. " << endl;
322 cout << " -MBB filename same with -Mbb but with BB file " << endl;
323 cout << "" << endl;
324 cout << " -errg (double) Set the level of error on geometry (0.1) by default" << endl;
325 cout << " -coef (double) Set the value of mutiplicative " << endl;
326 cout << " coef on the mesh size (1 by default)." << endl;
327 cout << " -power (double) Set the power parameter of hessien to construct " << endl;
328 cout << " the metric (1 by default)" << endl;
329 cout << " -maxsubdiv (double) Change the metric such that the maximal subdivision "
330 << endl;
331 cout << " of a background's edge is bound by the " << endl;
332 cout << " given number (always limited by 10) " << endl;
333 cout << " -ratio (double) Set the ratio for a smoothing of the metric." << endl;
334 cout << " If ratio is 0 then no smoothing" << endl;
335 cout << " and if ratio is in [1.1,10] then the smoothing " << endl;
336 cout << " change the metrix such that the greatest geometrical" << endl;
337 cout << " progression (speed of mesh size variation) " << endl;
338 cout << " in a mesh is bounded by ratio (by default no smoothing)."
339 << endl;
340 cout << " -hmin (double) Set the value of the minimal edge size." << endl;
341 cout << " -hminaniso (double) Set the value of the minimal edge size and save aniso."
342 << endl;
343 cout << " -hmax (double) Set the value of the maximal edge size." << endl;
344 cout << " -NbSmooth (int) Number of Smoothing iteration " << endl;
345 cout << " (3 by default if the metric is set otherwise 0)." << endl;
346 cout << " -omega (double) relaxation parameter for Smoothing " << endl;
347 cout << " -splitpbedge split in 2 all internal edges with 2 vertex on boundary"
348 << endl;
349 cout << " -nosplitpbedge d'ont cut internal edges with 2 vertex on boundary (default)"
350 << endl;
351 cout << "" << endl;
352 cout << "" << endl;
353 cout << " the next arguments are used with the -Mbb argument" << endl;
354 cout << "" << endl;
355 cout << " -KeepBackVertices Try to Keep old vertices (default) " << endl;
356 cout << " -noKeepBackVertices all vertices are create from scratch " << endl;
357 cout << " -err (double) Set the level of error (default 0.01)" << endl;
358 cout << " -iso The constructed metric must be isotropic" << endl;
359 cout << " -aniso The constructed metric can be anisotropic " << endl;
360 cout << " -anisomax (double) Set maximum ratio of anisotropy " << endl;
361 cout << " 1 => isotropic (1e6 by default) " << endl;
362 cout << " -RelError Construct metric with relative error " << endl;
363 cout << " -AbsError Construct metric with with abs error " << endl;
364 cout << " -CutOff (double) Set the limit of the relative error (1.e-5)" << endl;
365 cout << " -NbJacobi (int) Set the number of Jacobi iteration for smoothing " << endl;
366 cout << " the construction of metric (1 by default)." << endl;
367 cout << " -NoRescaling Don't rescaling of all solution between [0..1] "
368 << "before metric computation " << endl;
369 cout << " -Rescaling Do rescaling of all solution between [0..1] "
370 << "before metric computation ((default case)" << endl;
371 cout << " -H (int) choices for computing the hessein (test)" << endl;
372 cout << "" << endl;
373 cout << "" << endl;
374 cout << " Definition of some internal variable and limitation." << endl;
375 cout << "" << endl;
376 cout << " -v (int) Set the level of impression (verbosity) " << endl;
377 cout << " between 0 and 10 (1 by default)." << endl;
378 cout << " -nbv (int) Set the maximal of vertices (" << nbvx << " by default)."
379 << endl;
380 cout << "" << endl;
381 cout << "" << endl;
382 cout << " To interpoled a solution form background mesh to generated " << endl;
383 cout << " mesh (in case of adpatation)" << endl;
384 cout << "" << endl;
385 cout << " -rbb filename Read solution file defined on the background " << endl;
386 cout << " mesh for interpolation on created mesh." << endl;
387 cout << " bb file. (by default the -Mbb filename)" << endl;
388 cout << " -wbb filename Write the file of interpolation of the solutions" << endl;
389 cout << " read with -rbb argument. " << endl;
390 cout << " bb file." << endl;
391 cout << " -rBB filename same -rbb but with BB file ";
392 cout << " -wBB filename same -wbb but with BB file ";
393 cout << "" << endl;
394 cout << " Output mesh file for adpation or generation." << endl;
395 cout << " Remark: one of output mesh file is require " << endl;
396 cout << endl;
397 cout << " -o filename Create a DB Mesh file." << endl;
398 cout << " -oamdba filename Create a amdba file." << endl;
399 cout << " -oftq filename Create a ftq file." << endl;
400 cout << " -omsh filename Create a msh file (freefem3 file)." << endl;
401 cout << " -oam_fmt filename Create a am_fmt file. " << endl;
402 cout << " -oam filename Create a am file. " << endl;
403 cout << " -onopo filename Create a nopo file. " << endl;
404 cout << " -oM filename Create a metric file. " << endl;
405
406 cout << endl;
407 cout << endl;
408 cout << " -thetaquad (double) minimal angle of a quadrangle " << endl;
409 cout << " -2q split triangles in 3 quad and quad in 4 quad " << endl;
410 cout << " -2 split triangles in 4 trai and quad in 4 quad " << endl;
411
412 cout << endl;
413
414 cout << " Remark: if no argument is given when the arguments are read in " << datafile
415 << " file" << endl;
416 exit(3);
417 }
418 // some verification
419 NoMeshReconstruction = fmeshr != 0;
420 if (!fmeshback) fmeshback = fmeshr;
421 fileout = fmeshout || fam || fnopo || fftq || fam || famdba || famfmt || wbb || wBB;
422 if (!fileout && !foM) {
423 cerr << " No Output file a given " << endl;
424 MeshError(1);
425 }
426 if (maxsubdiv > boundmaxsubdiv || maxsubdiv <= 1.0) {
427 cerr << " -maxsubdiv " << maxsubdiv << " is not in ]1," << boundmaxsubdiv << "]" << endl;
428 exit(3);
429 }
430 if (iso) anisomax = 1;
431 if (!(fmetrix || fMbb)) NbSmooth = 0; // no metric -> no smoothing
432 if (!rbb) // to set the rbb filename by default
433 rbb = fMbb;
434 if (!rBB) // to set the rbb filename by default
435 rBB = fMBB;
436
437 if (fMbb && rbb) rbbeqMbb = !strcmp(rbb, fMbb);
438 if (fMBB && rBB) rBBeqMBB = !strcmp(rBB, fMBB);
439
440 #ifdef DRAWING
441 if (initgraphique) {
442 initgraphique( );
443 initgraph = 1;
444 withrgraphique = 0;
445 }
446 #endif
447
448 if (verbosity)
449 cout << bamgversion << " : ";
450 if (fgeom && fileout)
451 cout << " Construction of a mesh from the geometry : " << fgeom << endl
452 << " with less than " << nbvx << " vertices " << endl;
453 else if (fmeshback && fileout)
454 if (NoMeshReconstruction)
455 cout << " Modification of a adpated mesh " << fmeshback << " with less than " << nbvx
456 << " vertices" << endl;
457 else
458 cout << " Construction of a adpated mesh from the background mesh: " << fmeshback
459 << " with less than " << nbvx << " vertices" << endl;
460 else if (fmeshback && foM)
461 cout << " Construction of the metric file on the background mesh: " << fmeshback << endl;
462
463 if (fgeom && fileout) {
464 double time1;
465 if (verbosity)
466 cout << " Construction of a mesh from the geometry : " << fgeom << endl
467 << " with less than " << nbvx << " vertices " << endl;
468
469 time0 = CPUtime( );
470 Geometry Gh(fgeom);
471 hmin = Max(hmin, Gh.MinimalHmin( ));
472 hmax = Min(hmax, Gh.MaximalHmax( ));
473 if (fmetrix)
474 Gh.ReadMetric(fmetrix, hmin, hmax, coef);
475 else {
476 Int4 iv;
477 for (iv = 0; iv < Gh.nbv; iv++) {
478 MetricAnIso M = Gh[iv];
479 MatVVP2x2 Vp(M / coef);
480 Vp.Maxh(hmin);
481 Vp.Minh(hmax);
482 Gh.vertices[iv].m = Vp;
483 }
484 }
485
486 time1 = CPUtime( );
487 if (verbosity > 3) cout << " Cpu for read the geometry " << time1 - time0 << "s" << endl;
488 Triangles Th(nbvx, Gh);
489 if (costheta <= 1) Th.MakeQuadrangles(costheta);
490 if (allquad) Th.SplitElement(allquad == 2);
491
492 if (SplitEdgeWith2Boundary) Th.SplitInternalEdgeWithBorderVertices( );
493 Th.ReNumberingTheTriangleBySubDomain( );
494 time2 = CPUtime( );
495 if (verbosity > 3) Th.ShowHistogram( );
496 if (NbSmooth > 0) Th.SmoothingVertex(NbSmooth, omega);
497 if (verbosity > 3 && NbSmooth > 0) Th.ShowHistogram( );
498 #ifdef DRAWING
499 if (initgraph) {
500 Th.InitDraw( );
501 Th.Draw( );
502 Th.inquire( );
503 }
504 #endif
505 time3 = CPUtime( );
506 if (verbosity > 1)
507 cout << " Cpu for meshing :" << setw(8) << time2 - time1 << "s, for Smoothing "
508 << time3 - time2 << "s Nb Vertices/s = " << (Th.nbv) / (time2 - time1) << setw(0)
509 << endl;
510 if (fmeshout) Th.Write(fmeshout, Triangles::BDMesh);
511 if (famfmt) Th.Write(famfmt, Triangles::am_fmtMesh);
512 if (fam) Th.Write(fam, Triangles::amMesh);
513 if (famdba) Th.Write(famdba, Triangles::amdbaMesh);
514 if (fftq) Th.Write(fftq, Triangles::ftqMesh);
515 if (fmsh) Th.Write(fmsh, Triangles::mshMesh);
516 if (fnopo) Th.Write(fnopo, Triangles::NOPOMesh);
517 #ifdef DRAWING
518 if (initgraph) {
519
520 reffecran( );
521 Th.InitDraw( );
522 Th.Draw( );
523 Th.inquire( );
524 }
525 #endif
526
527 time3 = CPUtime( );
528 if (verbosity > 0) {
529 cout << " Cpu for meshing with io :" << setw(8) << time3 - time0
530 << "s Nb Vertices/s = " << (Th.nbv) / (time3 - time0) << setw(0) << endl;
531 cout << " Nb vertices = " << Th.nbv;
532 if (Th.nbt - Th.NbOutT - Th.NbOfQuad * 2)
533 cout << " Nb Triangles = " << (Th.nbt - Th.NbOutT - Th.NbOfQuad * 2);
534 if (Th.NbOfQuad) cout << " Nb Quadrilaterals = " << Th.NbOfQuad;
535 cout << endl;
536 }
537
538 } else if (fmeshback && (fmetrix || fMbb || fMBB || NoMeshReconstruction) &&
539 (fileout || foM || (rbb && wbb) || (rBB && wBB))) {
540 // to be sure the value was not stupide
541
542 time0 = CPUtime( );
543 Triangles BTh(fmeshback, cutoffradian); // read the background mesh
544 hmin = Max(hmin, BTh.MinimalHmin( ));
545 hmax = Min(hmax, BTh.MaximalHmax( ));
546
547 BTh.MakeQuadTree( );
548 // time1 = CPUtime( );
549 if (fmetrix)
550 BTh.ReadMetric(fmetrix, hmin, hmax, coef);
551 else { // init with hmax
552 Metric Mhmax(hmax);
553 for (Int4 iv = 0; iv < BTh.nbv; iv++) BTh[iv].m = Mhmax;
554 }
555 if (fMbb) {
556 solMbb = ReadbbFile(fMbb, nbsolbb, lsolbb, 2, 2);
557 if (lsolbb != BTh.nbv) {
558 cerr << "fatal error nbsol " << nbsolbb << " " << lsolbb << " =! " << BTh.nbv << endl;
559 cerr << " size of sol incorrect " << endl;
560 MeshError(99);
561 }
562 assert(lsolbb == BTh.nbv);
563 BTh.IntersectConsMetric(solMbb, nbsolbb, 0, hmin, hmax, sqrt(err) * coef, 1e30,
564 AbsError ? 0.0 : CutOff, nbjacoby, Rescaling, power, ChoiseHessien);
565 if (!rbbeqMbb) delete[] solMbb, solMbb = 0;
566 }
567 if (fMBB) {
568 solMBB = ReadBBFile(fMBB, nbsolBB, lsolBB, typesolsBB, 2, 2);
569 if (lsolBB != BTh.nbv) {
570 cerr << "fatal error nbsol " << nbsolBB << " " << lsolBB << " =! " << BTh.nbv << endl;
571 cerr << " size of sol incorrect " << endl;
572 MeshError(99);
573 }
574 assert(lsolBB == BTh.nbv);
575 BTh.IntersectConsMetric(solMBB, nbsolBB, 0, hmin, hmax, sqrt(err) * coef, 1e30,
576 AbsError ? 0.0 : CutOff, nbjacoby, Rescaling, ChoiseHessien);
577 if (!rBBeqMBB) delete[] solMBB, solMBB = 0;
578 }
579
580 BTh.IntersectGeomMetric(errg, iso);
581 if (raison) BTh.SmoothMetric(raison);
582 BTh.MaxSubDivision(maxsubdiv);
583 //
584 if (iso) anisomax = 1;
585 BTh.BoundAnisotropy(anisomax, hminaniso);
586 time2 = CPUtime( );
587 if (foM) {
588 if (verbosity > 2) cout << " -- write Metric file " << foM << endl;
589
590 ofstream f(foM);
591 if (f) BTh.WriteMetric(f, iso);
592 }
593
594 if (fileout) {
595 if (NoMeshReconstruction)
596 if ((fmeshback == fmeshr) || (!strcmp(fmeshback, fmeshr)))
597 Thr = &BTh, Thb = 0; // back and r mesh are the same
598 else
599 Thr = new Triangles(fmeshr, cutoffradian), Thb = &BTh; // read the new
600
601 Triangles &Th(*(NoMeshReconstruction
602 ? new Triangles(*Thr, &Thr->Gh, Thb,
603 nbvx) // copy the mesh + free space to modification
604 : new Triangles(nbvx, BTh, KeepBackVertices) // construct a new mesh
605 ));
606 if (Thr != &BTh) delete Thr;
607
608 if (costheta <= 1.0) Th.MakeQuadrangles(costheta);
609 if (allquad) Th.SplitElement(allquad == 2);
610 if (SplitEdgeWith2Boundary) Th.SplitInternalEdgeWithBorderVertices( );
611 Th.ReNumberingTheTriangleBySubDomain( );
612 time3 = CPUtime( );
613
614 if (verbosity > 0)
615 cout << " Cpu time for meshing : " << time3 - time2
616 << "s Nb Vertices/s = " << (Th.nbv) / (time3 - time2) << endl;
617 if (verbosity > 3) Th.ShowHistogram( );
618 if (NbSmooth > 0) Th.SmoothingVertex(NbSmooth, omega);
619 if (verbosity > 3 && NbSmooth > 0) Th.ShowHistogram( );
620 #ifdef DRAWING
621 if (initgraph) {
622 Th.InitDraw( );
623 Th.Draw( );
624 Th.inquire( );
625 }
626 #endif
627 Th.BTh.ReMakeTriangleContainingTheVertex( );
628
629 if (fmeshout) Th.Write(fmeshout, Triangles::BDMesh);
630 if (famfmt) Th.Write(famfmt, Triangles::am_fmtMesh);
631 if (fam) Th.Write(fam, Triangles::amMesh);
632 if (famdba) Th.Write(famdba, Triangles::amdbaMesh);
633 if (fftq) Th.Write(fftq, Triangles::ftqMesh);
634 if (fmsh) Th.Write(fmsh, Triangles::mshMesh);
635 if (fnopo) Th.Write(fnopo, Triangles::NOPOMesh);
636
637 if ((rbb && wbb) || (rBB && wBB)) // the code for interpolation
638 {
639 if (verbosity > 1) {
640 if (rbb) cout << " -- interpolation P1 files : " << rbb << " -> " << wbb << endl;
641 if (rBB) cout << " -- interpolation P1 files : " << rBB << " -> " << wBB << endl;
642 }
643 double time00 = CPUtime( );
644 const int dim = 2;
645 // optimisation read only si rbb != fMbb
646
647 double *solbb = 0;
648 double *solBB = 0;
649
650 if (rbb) solbb = rbbeqMbb ? solMbb : ReadbbFile(rbb, nbsolbb, lsolbb, 2, 2);
651 if (rBB) solBB = rBBeqMBB ? solMBB : ReadBBFile(rBB, nbsolBB, lsolBB, typesolsBB, 2, 2);
652
653 // cout << " " << rbbeqMbb << " " << sol << " " << nbsol << " " << lsol << endl;
654 if (!solBB && !solbb) {
655 cerr << " Fatal Error "
656 << "solBB = " << solBB << " solbb= " << solbb << endl;
657 exit(2);
658 }
659
660 ofstream *fbb = wbb ? new ofstream(wbb) : 0;
661 ofstream *fBB = wBB ? new ofstream(wBB) : 0;
662 Int4 nbfieldBB = 0, nbfieldbb = nbsolbb;
663 if (fbb) *fbb << dim << " " << nbsolbb << " " << Th.nbv << " " << 2 << endl;
664 if (fBB) {
665 // int i;
666 *fBB << dim << " " << nbsolBB;
667 for (i = 0; i < nbsolBB; i++) *fBB << " " << (typesolsBB ? typesolsBB[i] + 1 : 1);
668 *fBB << " " << Th.nbv << " " << 2 << endl;
669 // assert(dim == 2); // dim is defined as const int = 2
670 for (i = 0; i < nbsolBB; i++) nbfieldBB += typesolsBB ? typesolsBB[i] + 1 : 1;
671 // this code is good only if dim == 2
672 }
673 cout << "nb of field BB " << nbfieldBB << endl;
674 // if(fBB) fBB->precision(15);
675 // if(fbb) fbb->precision(15);
676 for (i = 0; i < Th.nbv; i++) {
677 Int4 i0, i1, i2;
678 double a[3];
679 Icoor2 dete[3];
680 I2 I = Th.BTh.toI2(Th.vertices[i].r);
681 Triangle &tb = *Th.BTh.FindTriangleContening(I, dete);
682
683 if (tb.det > 0) { // internal point
684 a[0] = (Real8)dete[0] / tb.det;
685 a[1] = (Real8)dete[1] / tb.det;
686 a[2] = (Real8)dete[2] / tb.det;
687 i0 = Th.BTh.Number(tb[0]);
688 i1 = Th.BTh.Number(tb[1]);
689 i2 = Th.BTh.Number(tb[2]);
690 } else {
691 double aa, bb;
692
693 TriangleAdjacent ta = CloseBoundaryEdge(I, &tb, aa, bb).Adj( );
694
695 int k = ta;
696 Triangle &tc = *(Triangle *)ta;
697 i0 = Th.BTh.Number(tc[0]);
698 i1 = Th.BTh.Number(tc[1]);
699 i2 = Th.BTh.Number(tc[2]);
700 a[VerticesOfTriangularEdge[k][1]] = aa;
701 a[VerticesOfTriangularEdge[k][0]] = bb;
702 a[OppositeVertex[k]] = 1 - aa - bb;
703 }
704
705 Int4 ibb0 = nbfieldbb * i0;
706 Int4 ibb1 = nbfieldbb * i1;
707 Int4 ibb2 = nbfieldbb * i2;
708
709 Int4 iBB0 = nbfieldBB * i0;
710 Int4 iBB1 = nbfieldBB * i1;
711 Int4 iBB2 = nbfieldBB * i2;
712 Int4 j = 0;
713
714 for (j = 0; j < nbfieldbb; j++)
715 *fbb << " " << (a[0] * solbb[ibb0++] + a[1] * solbb[ibb1++] + a[2] * solbb[ibb2++]);
716
717 for (j = 0; j < nbfieldBB; j++)
718 *fBB << " " << (a[0] * solBB[iBB0++] + a[1] * solBB[iBB1++] + a[2] * solBB[iBB2++]);
719
720 if (fbb) *fbb << endl;
721 if (fBB) *fBB << endl;
722 }
723 if (fbb) delete fbb; // close
724 if (fBB) delete fBB; // close
725 if (solbb) delete[] solbb;
726 if (solBB) delete[] solBB;
727 double time04 = CPUtime( );
728 if (verbosity > 0)
729 cout << " Cpu time for interpolation " << time04 - time00 << " s" << endl;
730 }
731 time4 = CPUtime( );
732 if (verbosity > 0) {
733 cout << " Cpu time for meshing with io : " << time4 - time0
734 << "s Nb Vertices/s = " << Th.nbv / (time4 - time0) << endl
735 << " Nb vertices = " << Th.nbv;
736 if (Th.nbt - Th.NbOutT - Th.NbOfQuad * 2)
737 cout << " Nb Triangles = " << (Th.nbt - Th.NbOutT - Th.NbOfQuad * 2);
738 if (Th.NbOfQuad) cout << " Nb Quadrilaterals = " << Th.NbOfQuad;
739 cout << endl;
740 }
741
742 #ifdef DRAWING
743 if (initgraph) {
744 reffecran( );
745 Th.InitDraw( );
746 Th.Draw( );
747 Th.inquire( );
748 }
749 #endif
750
751 // cout << "delete &Th " ;
752 delete &Th;
753 // cout << " end " << endl;
754 }
755 } // if (fgeom && fmeshout)
756 // clean the
757 for (i = 1; i < datargc; i++) delete[] datargv[i];
758
759 #ifdef DRAWING
760 if (initgraph) closegraphique( );
761 initgraph = 0;
762 #endif
763 cout << flush;
764 return (0);
765 }
766