msbGrid  1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gridfactory.hh
Go to the documentation of this file.
1 /*****************************************************************************
2 * This program is part of the msbGrid software *
3 * *
4 * msbGrid stands for multi-structured block Grid generator *
5 * *
6 * msbGrid is a free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 2 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * msbGrid is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * See the file COPYING for full copying permissions. *
17 *****************************************************************************/
23 // TODO implement run time evaluation
24 // TODO implement overall structured superposed grid
25 
26 #ifndef _GRID_FACTORY_HH_
27 #define _GRID_FACTORY_HH_
28 
29 #include <sstream>
30 
31 #include <block.hh>
32 #include <io.hh>
33 #include <stepper.hh>
34 
35 namespace msbGrid
36 {
37 template< typename CellT >
38 class GridFactory;
39 
40 const std :: string newCl(const std :: string &clName, const Scalar &clValue);
41 }
42 
43 template< typename CellT >
44 class msbGrid :: GridFactory
45 {
46 public:
47  typedef CellT Cell;
48  typedef msbGrid :: GridFactory< Cell > ThisT;
49 
50  enum { dim = Cell :: dim };
51  typedef msbGrid :: Color Color;
52 
53  typedef msbGrid :: Scalar Scalar;
54  typedef msbGrid :: InputData InputData;
55  typedef msbGrid :: Point Point;
56  typedef msbGrid :: Edge Edge;
57  typedef msbGrid :: Triangle Triangle;
58  typedef msbGrid :: Quadrangle Quadrangle;
59  typedef msbGrid :: Stepper< ThisT > Stepper;
60  typedef msbGrid :: PointsMultiBlock PointsMultiBlock;
61  typedef msbGrid :: EdgesBlock EdgesBlock;
62  typedef msbGrid :: CellsBlock< Triangle > TrianglesBlock;
63  typedef msbGrid :: CellsBlock< Quadrangle > QuadranglesBlock;
64 
65  typedef msbGrid :: MFilesWriter MFilesWriter;
66  typedef msbGrid :: GeoFilesWriter GeoFilesWriter;
67 
68 private:
69  const InputData *inData_;
70 
71  int totCelNumber_ /*, totQdNumber_*/;
72  std :: string outputDirName_;
73 
74  msbGrid :: PointsBlock blPtExtBound_;
75  msbGrid :: EdgesBlock blEdExtBound_;
76 
77  PointsMultiBlock vBlPtI_, vBlPtJ_;
78  std :: vector< EdgesBlock > vBlEd_;
79  std :: vector< TrianglesBlock > vBlTr_;
80  std :: vector< QuadranglesBlock > vBlQd_;
81 
82  std :: vector< std :: vector< Point * > > vvBoundPtPtr_; /*vector of Point instances on block boundaries*/
83  std :: vector< std :: vector< Edge * > > vvBlBoundEdPtr_; /*vector of Edge instances on block boundaries*/
84 
85  std :: vector< std :: vector< Triangle * > > vvBoundTrPtr_; /* OnBlockBoundaryTriangles */
86  std :: vector< std :: vector< Quadrangle * > > vvBoundQdPtr_; /* OnBlockBoundaryQuadrangles */
87 
88  std :: vector< Edge * > vEdExtBound_; /*vector of Edge instances on the external boundary*/
89 
90  /* stepper */
91  Stepper stepper_;
92 
93 public:
94 
95  GridFactory(const InputData * inData0)
96  {
97  inData_ = inData0;
98 
99  init_();
100  doWork_();
101  }
102 
103  ~GridFactory() {}
104 
105 private:
106  void init_()
107  {
108  totCelNumber_ = /*totQdNumber_ =*/ 0;
109  outputDirName_ = outputDirectoryName_( inData_->outputDirName() );
110 
111  vBlPtI_.init(inData_, msbGrid :: Color :: Red);
112  vBlPtJ_.init(inData_, msbGrid :: Color :: Green);
113  stepper_.init( inData_->blockNumber() );
114  // stepper_.iStepMax() = 11;
115  }
116 
117  const std :: string outputDirectoryName_(const std :: string outpuDirName0) const
118  {
119 #if VERBOSE_LEVEL >= 1
120  std :: cout << "-- output directory name\n";
121 #endif
122 
123  if ( outpuDirName0 != "" ) {
124 #if VERBOSE_LEVEL >= 1
125  std :: cout << "The output directory name is explicitly given : " << outpuDirName0 << "\n";
126 #endif
127  return outpuDirName0;
128  }
129 
130  std :: string xMinString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->xMin() ) )->str();
131  std :: string xMaxString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->xMax() ) )->str();
132 
133  std :: string bnCoeffxString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->blockNumberController() [ 0 ] ) )->str();
134  std :: string bnCoeffyString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->blockNumberController() [ 1 ] ) )->str();
135 
136  std :: string drString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->dr() ) )->str();
137  std :: string rafCoeffString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->rafCoeff() ) )->str();
138 
139  std :: string blockCenterRandSeqIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->blockCenterRandSeqId() ) )->str();
140  std :: string blockAngleRandSeqIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->blockAngleRandSeqId() ) )->str();
141 
142  std :: string blockCenterPertCoeffString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->blockCenterPertCoeff() ) )->str();
143 
144  std :: string dist2BoundCoeffString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->dist2BoundCoeff() ) )->str();
145  std :: string interBlockDistCoeffString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->interBlockDistCoeff() ) )->str();
146 
147  std :: string outpuDirName1 = "xMin" + xMinString + "_"
148  + "xMax" + xMaxString + "_"
149  + "bnCx" + bnCoeffxString + "_"
150  + "bnCy" + bnCoeffyString + "_"
151  + "dr" + drString + "_"
152  + "rafC" + rafCoeffString + "_"
153  + "bCId" + blockCenterRandSeqIdString + "_"
154  + "bAId" + blockAngleRandSeqIdString + "_"
155  + "bCPC" + blockCenterPertCoeffString + "_"
156  + "d2bC" + dist2BoundCoeffString + "_"
157  + "intbC" + interBlockDistCoeffString;
158 
159 #if VERBOSE_LEVEL >= 1
160  std :: cout << "The output directory name is computed : " << outpuDirName1 << "\n";
161 #endif
162 
163  return outpuDirName1;
164  }
165 
166  /* fill the domain "Gamma" by geometric entities: points, edges and cells */
167  void doWork_()
168  {
169 #if OUTPUT_MFILES_DIR_CREATED == 0
170  const std :: string createDirectoryCommandString = "mkdir -p " + inData_->outputPath() + "m/" + outputDirName_;
171  char *createDirectoryCommandChar = ( char * ) createDirectoryCommandString.c_str();
172  system(createDirectoryCommandChar);
173  #define OUTPUT_MFILES_DIR_CREATED 1
174 #endif
175  /*
176  * initialize points, edges and cells development mode identifiers (DEV_MOD),
177  * noting that real (or final) identifiers of geometric entities are set by setIdentifiers_()
178  */
179 #if DEV_MOD == 1
180  unsigned int pointsMaxId = 1, edgesMaxId = 1, cellsQdMaxId = 1;
181 #endif
182 
183  /* points on the external boundary (bounding box) */ // TODO change this to handle 3D
184 #if CLOCKWISE_ROTATION_SENSE == 0
185  blPtExtBound_.add( Point( inData_->xMin(), inData_->yMin( inData_->xMin() ) ) );
186  blPtExtBound_.add( Point( inData_->xMin(), inData_->yMax( inData_->xMin() ) ) );
187  blPtExtBound_.add( Point( inData_->xMax(), inData_->yMax( inData_->xMax() ) ) );
188  blPtExtBound_.add( Point( inData_->xMax(), inData_->yMin( inData_->xMax() ) ) );
189 #endif
190 #if CLOCKWISE_ROTATION_SENSE == 1
191  blPtExtBound_.add( Point( inData_->xMin(), inData_->yMin( inData_->xMin() ) ) );
192  blPtExtBound_.add( Point( inData_->xMax(), inData_->yMin( inData_->xMax() ) ) );
193  blPtExtBound_.add( Point( inData_->xMax(), inData_->yMax( inData_->xMax() ) ) );
194  blPtExtBound_.add( Point( inData_->xMin(), inData_->yMax( inData_->xMin() ) ) );
195 #endif
196 
197 #if DEV_MOD == 1
198  pointsMaxId = blPtExtBound_.putId(pointsMaxId);
199  blPtExtBound_.putMarkerProperties(/*symbol=*/ "'s'"
200  , /*size=*/ 14
201  , /*edgeColor=*/ Color :: Black
202  , /*faceColor=*/ Color :: Magenta);
203 #endif
204 
205  //#define ENABLE_DEV_MODE 1
206 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 1 ) )
207  MFilesWriter mEnities(inData_);
208  mEnities.openFile(outputDirName_ + "/mEnities");
209  mEnities.write(blPtExtBound_, true);
210  mEnities.closeFile();
211  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
212  std :: cerr << "exiting the program by calling : std::exit(1);\n";
213  std :: exit(1);
214 #endif
215 
216  /* edges on the external boundary (bounding box) */
217  blEdExtBound_.add( Edge(& blPtExtBound_.vEnt() [ 0 ], & blPtExtBound_.vEnt() [ 1 ]) );
218  blEdExtBound_.add( Edge(& blPtExtBound_.vEnt() [ 1 ], & blPtExtBound_.vEnt() [ 2 ]) );
219  blEdExtBound_.add( Edge(& blPtExtBound_.vEnt() [ 2 ], & blPtExtBound_.vEnt() [ 3 ]) );
220  blEdExtBound_.add( Edge(& blPtExtBound_.vEnt() [ 3 ], & blPtExtBound_.vEnt() [ 0 ]) );
221 
222  vEdExtBound_.clear(); /* vector of edge pointer needed for *.geo files output, TODO get rid of this confusing design */
223  for ( unsigned int ieExt = 0; ieExt < blEdExtBound_.size(); ++ieExt ) {
224  vEdExtBound_.push_back(& blEdExtBound_.vEnt() [ ieExt ]);
225  vEdExtBound_.back()->onBlockBoundary() = true;
226  }
227 
228 #if DEV_MOD == 1
229  edgesMaxId = blEdExtBound_.putId(edgesMaxId);
230  blEdExtBound_.putColor(Color :: Blue);
231 #endif
232 
233  //#define ENABLE_DEV_MODE 2
234 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 2 ) )
235  MFilesWriter mEnities(inData_);
236  mEnities.openFile(outputDirName_ + "/mEnities");
237  mEnities.write(blPtExtBound_, true);
238  mEnities.write(blEdExtBound_, true);
239  mEnities.closeFile();
240  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
241  std :: cerr << "exiting the program by calling : std::exit(1);\n";
242  std :: exit(1);
243 #endif
244 
245  vBlPtI_.fill(stepper_);
246  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
247 #if DEV_MOD == 1
248  pointsMaxId = vBlPtI_ [ iBl ].putId(pointsMaxId);
249  vBlPtI_ [ iBl ].putMarkerProperties(/*symbol=*/ "'s'"
250  , /*size=*/ 14
251  , /*edgeColor=*/ Color :: Black
252  , /*faceColor=*/ Color :: Red);
253 #endif
254  }
255 
256  //#define ENABLE_DEV_MODE 3
257 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 3 ) )
258  MFilesWriter mEnities(inData_);
259  mEnities.openFile(outputDirName_ + "/mEnities");
260  mEnities.write(blPtExtBound_, true);
261  mEnities.write(blEdExtBound_, true);
262  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
263  mEnities.write(vBlPtI_ [ iBl ], true);
264  }
265 
266  mEnities.closeFile();
267  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
268  std :: cerr << "exiting the program by calling : std::exit(1);\n";
269  std :: exit(1);
270 #endif
271 
272  stepper_.reset();
273  vBlPtJ_.fill(stepper_);
274  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
275 #if DEV_MOD == 1
276  pointsMaxId = vBlPtJ_ [ iBl ].putId(pointsMaxId);
277  vBlPtJ_ [ iBl ].putMarkerProperties(/*symbol=*/ "'s'"
278  , /*size=*/ 10
279  , /*edgeColor=*/ Color :: Black
280  , /*faceColor=*/ Color :: Green);
281 #endif
282  }
283 
284  //#define ENABLE_DEV_MODE 4
285 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 4 ) )
286  MFilesWriter mEnities(inData_);
287  mEnities.openFile(outputDirName_ + "/mEnities");
288  mEnities.write(blPtExtBound_, true);
289  mEnities.write(blEdExtBound_, true);
290  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
291  mEnities.write(vBlPtI_ [ iBl ], true);
292  mEnities.write(vBlPtJ_ [ iBl ], false);
293  }
294 
295  mEnities.closeFile();
296  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
297  std :: cerr << "exiting the program by calling : std::exit(1);\n";
298  std :: exit(1);
299 #endif
300 
301  /* fill vpNeighbor0_ and vpNeighbor1_ of each point of vBlPtI_ and vBlPtJ_ */
302  fillvpNeighbor_();
303 
304  setOnBlockBoundaryPtQd1_();
305 
306  //#define ENABLE_DEV_MODE 5
307 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 5 ) )
308  MFilesWriter mEnities(inData_);
309  mEnities.openFile(outputDirName_ + "/mEnities");
310  mEnities.write(blPtExtBound_, true);
311  mEnities.write(blEdExtBound_, true);
312  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
313  mEnities.write(vBlPtI_ [ iBl ], true);
314  mEnities.write(vBlPtJ_ [ iBl ], false);
315  }
316 
317  mEnities.closeFile();
318  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
319  std :: cerr << "exiting the program by calling : std::exit(1);\n";
320  std :: exit(1);
321 #endif
322 
323  markPtQd1_();
324 
325  //#define ENABLE_DEV_MODE 6
326 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 6 ) )
327  MFilesWriter mEnities(inData_);
328  mEnities.openFile(outputDirName_ + "/mEnities");
329  mEnities.write(blPtExtBound_, true);
330  mEnities.write(blEdExtBound_, true);
331  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
332  mEnities.write(vBlPtI_ [ iBl ], true);
333  mEnities.write(vBlPtJ_ [ iBl ], true);
334  }
335 
336  mEnities.closeFile();
337  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
338  std :: cerr << "exiting the program by calling : std::exit(1);\n";
339  std :: exit(1);
340 #endif
341 
342  /* fill vBlEd_ */
343  fillvBlEd_();
344 
345  /* clean vBlEd_ : mark then remove duplicated edges from vBlEd_ */
346  cleanvBlEd_();
347 
348  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
349 #if DEV_MOD == 1
350  edgesMaxId = vBlEd_ [ iBl ].putId(edgesMaxId);
351  vBlEd_ [ iBl ].putColor(Color :: Black);
352 #endif
353  }
354 
355  //#define ENABLE_DEV_MODE 7
356 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 7 ) )
357  MFilesWriter mEnities(inData_);
358  mEnities.openFile(outputDirName_ + "/mEnities");
359  mEnities.write(blPtExtBound_, true);
360  mEnities.write(blEdExtBound_, true);
361  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
362  mEnities.write(vBlPtI_ [ iBl ], true);
363  mEnities.write(vBlPtJ_ [ iBl ], false);
364  mEnities.write(vBlEd_ [ iBl ], true);
365  }
366 
367  mEnities.closeFile();
368  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
369  std :: cerr << "exiting the program by calling : std::exit(1);\n";
370  std :: exit(1);
371 #endif
372 
373  /* fill vEd_ (vector of edge instance addresses) of each point */
374  fillPointsVEd_();
375 
376  /* fill vBlTr_ */
377  fillvBlQd_();
378 
379  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
380 #if DEV_MOD == 1
381  cellsQdMaxId = vBlQd_ [ iBl ].putId(cellsQdMaxId);
382  vBlQd_ [ iBl ].putColor(Color :: Black);
383 #endif
384  }
385 
386  //#define ENABLE_DEV_MODE 8
387 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 8 ) )
388  MFilesWriter mEnities(inData_);
389  mEnities.openFile(outputDirName_ + "/mEnities");
390  mEnities.write(blPtExtBound_, true);
391  mEnities.write(blEdExtBound_, true);
392  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
393  mEnities.write(vBlPtI_ [ iBl ], true);
394  mEnities.write(vBlPtJ_ [ iBl ], false);
395  mEnities.write(vBlEd_ [ iBl ], false);
396  mEnities.write(vBlQd_ [ iBl ], true);
397  }
398 
399  mEnities.closeFile();
400  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
401  std :: cerr << "exiting the program by calling : std::exit(1);\n";
402  std :: exit(1);
403 #endif
404 
405  /* fill vCel_ (vector of cell instance addresses) of each point */
406  fillPointsVCel_();
407 
408  /* fill vCel_ (vector of cell instance addresses) of each edge */
409  fillEdgesVCel_();
410 
411  setOnBlockBoundaryCelQd1_(Black);
412 
413  //#define ENABLE_DEV_MODE 9
414 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 9 ) )
415  MFilesWriter mEnities(inData_);
416  mEnities.openFile(outputDirName_ + "/mEnities");
417  mEnities.write(blPtExtBound_, true);
418  mEnities.write(blEdExtBound_, true);
419  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
420  mEnities.write(vBlPtI_ [ iBl ], true);
421  mEnities.write(vBlPtJ_ [ iBl ], false);
422  mEnities.write(vBlEd_ [ iBl ], false);
423  mEnities.write(vBlQd_ [ iBl ], true);
424  }
425 
426  mEnities.closeFile();
427  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
428  std :: cerr << "exiting the program by calling : std::exit(1);\n";
429  std :: exit(1);
430 #endif
431 
432  markPtQd2_();
433 
434  setOnBlockBoundEdQd_();
435 
436  //#define ENABLE_DEV_MODE 10
437 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 10 ) )
438  MFilesWriter mEnities(inData_);
439  mEnities.openFile(outputDirName_ + "/mEnities");
440  mEnities.write(blPtExtBound_, true);
441  mEnities.write(blEdExtBound_, true);
442  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
443  mEnities.write(vBlPtI_ [ iBl ], true);
444  mEnities.write(vBlPtJ_ [ iBl ], false);
445  mEnities.write(vBlEd_ [ iBl ], false);
446  mEnities.write(vBlQd_ [ iBl ], true);
447  }
448 
449  mEnities.closeFile();
450  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
451  std :: cerr << "exiting the program by calling : std::exit(1);\n";
452  std :: exit(1);
453 #endif
454 
455  setOnBlockBoundaryPtQd2_();
456 
457  //#define ENABLE_DEV_MODE 11
458 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 11 ) )
459  MFilesWriter mEnities(inData_);
460  mEnities.openFile(outputDirName_ + "/mEnities");
461  mEnities.write(blPtExtBound_, true);
462  mEnities.write(blEdExtBound_, true);
463  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
464  mEnities.write(vBlPtI_ [ iBl ], true);
465  mEnities.write(vBlPtJ_ [ iBl ], false);
466  mEnities.write(vBlEd_ [ iBl ], false);
467  mEnities.write(vBlQd_ [ iBl ], true);
468  }
469 
470  mEnities.closeFile();
471  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
472  std :: cerr << "exiting the program by calling : std::exit(1);\n";
473  std :: exit(1);
474 #endif
475 
476  setOnBlockBoundaryCelQd2_();
477 
478  //#define ENABLE_DEV_MODE 12
479 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 12 ) )
480  MFilesWriter mEnities(inData_);
481  mEnities.openFile(outputDirName_ + "/mEnities");
482  mEnities.write(blPtExtBound_, true);
483  mEnities.write(blEdExtBound_, true);
484  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
485  mEnities.write(vBlPtI_ [ iBl ], true);
486  mEnities.write(vBlPtJ_ [ iBl ], false);
487  mEnities.write(vBlEd_ [ iBl ], false);
488  mEnities.write(vBlQd_ [ iBl ], true);
489  }
490 
491  mEnities.closeFile();
492  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
493  std :: cerr << "exiting the program by calling : std::exit(1);\n";
494  std :: exit(1);
495 #endif
496 
497  markCelQd1_(Magenta);
498 
499  //#define ENABLE_DEV_MODE 13
500 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 13 ) )
501  MFilesWriter mEnities(inData_);
502  mEnities.openFile(outputDirName_ + "/mEnities");
503  mEnities.write(blPtExtBound_, true);
504  mEnities.write(blEdExtBound_, true);
505  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
506  if ( iBl == 13 ) {
507  mEnities.write(vBlPtI_ [ iBl ], true);
508  mEnities.write(vBlPtJ_ [ iBl ], false);
509  mEnities.write(vBlEd_ [ iBl ], false);
510  mEnities.write(vBlQd_ [ iBl ], true);
511  } else {
512  mEnities.write(vBlPtI_ [ iBl ], false);
513  mEnities.write(vBlPtJ_ [ iBl ], false);
514  mEnities.write(vBlEd_ [ iBl ], false);
515  mEnities.write(vBlQd_ [ iBl ], false);
516  }
517  }
518 
519  for ( int iQd1 = 0; iQd1 < vBlQd_.size(); ++iQd1 ) {
520  if ( vBlQd_ [ 13 ].vEnt() [ iQd1 ].Id() == 7102 ) {
521  bool isolatedCell = true;
522  for ( unsigned int iEd1 = 0; iEd1 < 4; ++iEd1 ) {
523  std :: cout << "vBlQd_[13].vEnt()[iQd1].vEd()[iEd1]->vCel().size() = "
524  << vBlQd_ [ 13 ].vEnt() [ iQd1 ].vEd() [ iEd1 ]->vCel().size() << "\n";
525 
526  if ( vBlQd_ [ 13 ].vEnt() [ iQd1 ].vEd() [ iEd1 ]->vCel().size() != 1 ) {
527  isolatedCell = false;
528  continue;
529  }
530  }
531  }
532  }
533 
534  mEnities.closeFile();
535  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
536  std :: cerr << "exiting the program by calling : std::exit(1);\n";
537  std :: exit(1);
538 #endif
539 
540  /* fill block vvBlBoundEdPtr_ */
541  fillvvBlBoundEdPtr_();
542 
543  //#define ENABLE_DEV_MODE 14
544 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 14 ) )
545  MFilesWriter mEnities(inData_);
546  mEnities.openFile(outputDirName_ + "/mEnities");
547  mEnities.write(blPtExtBound_, true);
548  mEnities.write(blEdExtBound_, true);
549  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
550  mEnities.write(vBlPtI_ [ iBl ], true);
551  mEnities.write(vBlPtJ_ [ iBl ], false);
552  mEnities.write(vBlEd_ [ iBl ], false);
553  mEnities.write(vBlQd_ [ iBl ], true);
554  }
555 
556  mEnities.closeFile();
557  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
558  std :: cerr << "exiting the program by calling : std::exit(1);\n";
559  std :: exit(1);
560 #endif
561 
562  setIdentifiers_();
563 
564  //#define ENABLE_DEV_MODE 15
565 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 15 ) )
566  MFilesWriter mEnities(inData_);
567  mEnities.openFile(outputDirName_ + "/mEnities");
568  mEnities.write(blPtExtBound_, true);
569  mEnities.write(blEdExtBound_, true);
570  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
571  mEnities.write(vBlPtI_ [ iBl ], true);
572  mEnities.write(vBlPtJ_ [ iBl ], false);
573  mEnities.write(vBlEd_ [ iBl ], false);
574  mEnities.write(vBlQd_ [ iBl ], true);
575  }
576 
577  mEnities.closeFile();
578  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
579  std :: cerr << "exiting the program by calling : std::exit(1);\n";
580  std :: exit(1);
581 #endif
582 
583 #if ( ( ENABLE_MFILES_OUTPUT == 1 ) && ( WRITE_MFILES == 1 ) )
584  /* write *.m output files */
585  // outputMFiles_();
586 #endif
587 
588 #if ENABLE_GEOFILES_OUTPUT == 1
589  /* write *.geo output files */
590  outputGeoFiles_();
591 #endif
592 
593  //#define ENABLE_DEV_MODE 16
594 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 16 ) )
595  MFilesWriter mEnities(inData_);
596  mEnities.openFile(outputDirName_ + "/mEnities");
597  mEnities.write(blPtExtBound_, true);
598  mEnities.write(blEdExtBound_, true);
599  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
600  mEnities.write(vBlPtI_ [ iBl ], true);
601  // mEnities.write(vBlPtJ_ [ iBl ], false);
602  mEnities.write(vBlEd_ [ iBl ], false);
603  mEnities.write(vBlQd_ [ iBl ], true);
604  }
605 
606  mEnities.closeFile();
607  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
608  std :: cerr << "exiting the program by calling : std::exit(1);\n";
609  std :: exit(1);
610 #endif
611 
612  //#define ENABLE_DEV_MODE 17
613 #if ( ( DEV_MOD == 1 ) && ( ENABLE_DEV_MODE == 17 ) )
614  MFilesWriter mEnities(inData_);
615  mEnities.openFile(outputDirName_ + "/mEnities");
616  mEnities.write(blPtExtBound_, true);
617  mEnities.write(blEdExtBound_, true);
618  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
619  mEnities.write(vBlPtI_ [ iBl ], true);
620  // mEnities.write(vBlPtJ_ [ iBl ], false);
621  mEnities.write(vBlEd_ [ iBl ], false);
622  mEnities.write(vBlQd_ [ iBl ], true);
623  }
624 
625  mEnities.closeFile();
626  std :: cerr << "ENABLE_DEV_MODE = " << ENABLE_DEV_MODE << "\n";
627  std :: cerr << "exiting the program by calling : std::exit(1);\n";
628  std :: exit(1);
629 #endif
630  } /* void doWork_() */
631 
632  void fillvpNeighbor_()
633  {
634  /*
635  * handle vBlPtI_ : for each point from vBlPtI_ save all points from vBlPtJ_ which satisfy the neighborhood condition
636  * handle vBlPtJ_ : for each point from vBlPtJ_ save all points from vBlPtI_ which satisfy the neighborhood condition
637  */
638  fillvpNeighbor0_();
639 
640  /*
641  * handle vBlPtI_ : for each point from vBlPtI_ save all points from vBlPtI_ which satisfy the neighborhood condition
642  * handle vBlPtJ_ : for each point from vBlPtJ_ save all points from vBlPtJ_ which satisfy the neighborhood condition
643  */
644  fillvpNeighbor1_();
645  }
646 
647  void fillvpNeighbor0_()
648  {
649 #if VERBOSE_LEVEL >= 2
650  std :: cout << "-- fill vpNeighbor0_\n";
651 #endif
652  Scalar diag = inData_->dr() / sqrt(2.0);
653 
654  /* handle vBlPtI_ */
655  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
656  for ( unsigned int iPt = 0; iPt < vBlPtI_ [ iBl ].size(); ++iPt ) {
657  vBlPtI_ [ iBl ].vEnt() [ iPt ].vpNeighbor0().clear();
658  for ( unsigned int jPt = 0; jPt < vBlPtJ_ [ iBl ].size(); ++jPt ) {
659  Scalar diff = eucDistance(vBlPtI_ [ iBl ].vEnt() [ iPt ], vBlPtJ_ [ iBl ].vEnt() [ jPt ]);
660 
661  if ( std :: abs(diff - diag) <= ( +Constants :: TOLERANCE ) ) {
662  vBlPtI_ [ iBl ].vEnt() [ iPt ].vpNeighbor0().push_back(& vBlPtJ_ [ iBl ].vEnt() [ jPt ]);
663  }
664 
665  if ( vBlPtI_ [ iBl ].vEnt() [ iPt ].vpNeighbor0().size() == 4 ) {
666  break;
667  }
668  }
669  }
670  }
671 
672  /* handle vBlPtJ_ */
673  for ( unsigned int jBl = 0; jBl < inData_->blockNumber(); ++jBl ) {
674  for ( unsigned int jPt = 0; jPt < vBlPtJ_ [ jBl ].size(); ++jPt ) {
675  vBlPtJ_ [ jBl ].vEnt() [ jPt ].vpNeighbor0().clear();
676  for ( unsigned int iPt = 0; iPt < vBlPtI_ [ jBl ].size(); ++iPt ) {
677  Scalar diff = eucDistance(vBlPtJ_ [ jBl ].vEnt() [ jPt ], vBlPtI_ [ jBl ].vEnt() [ iPt ]);
678 
679  if ( std :: abs(diff - diag) <= ( +Constants :: TOLERANCE ) ) {
680  vBlPtJ_ [ jBl ].vEnt() [ jPt ].vpNeighbor0().push_back(& vBlPtI_ [ jBl ].vEnt() [ iPt ]);
681  }
682 
683  if ( vBlPtJ_ [ jBl ].vEnt() [ jPt ].vpNeighbor0().size() == 4 ) {
684  break;
685  }
686  }
687  }
688  }
689  }
690 
691  void fillvpNeighbor1_()
692  {
693 #if VERBOSE_LEVEL >= 2
694  std :: cout << "-- fill vpNeighbor1_\n";
695 #endif
696  Scalar diag = inData_->dr() / sqrt(2.0);
697 
698  /* handle vBlPtI_*/
699  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
700  for ( unsigned int iPt1 = 0; iPt1 < vBlPtI_ [ iBl ].size(); ++iPt1 ) {
701  vBlPtI_ [ iBl ].vEnt() [ iPt1 ].vpNeighbor1().clear();
702  for ( unsigned int iPt2 = 0; iPt2 < vBlPtI_ [ iBl ].size(); ++iPt2 ) {
703  Scalar diff = eucDistance(vBlPtI_ [ iBl ].vEnt() [ iPt1 ], vBlPtI_ [ iBl ].vEnt() [ iPt2 ]);
704 
705  if ( ( diff - diag ) >= ( -Constants :: TOLERANCE ) && ( diff - inData_->dr() ) <= ( +Constants :: TOLERANCE ) ) {
706  vBlPtI_ [ iBl ].vEnt() [ iPt1 ].vpNeighbor1().push_back(& vBlPtI_ [ iBl ].vEnt() [ iPt2 ]);
707  }
708 
709  if ( vBlPtI_ [ iBl ].vEnt() [ iPt1 ].vpNeighbor1().size() == 4 ) {
710  break;
711  }
712  }
713  }
714  }
715 
716  /* handle vBlPtJ_*/
717  for ( unsigned int jBl = 0; jBl < inData_->blockNumber(); ++jBl ) {
718  for ( unsigned int jPt1 = 0; jPt1 < vBlPtJ_ [ jBl ].size(); ++jPt1 ) {
719  vBlPtJ_ [ jBl ].vEnt() [ jPt1 ].vpNeighbor1().clear();
720  for ( unsigned int jPt2 = 0; jPt2 < vBlPtJ_ [ jBl ].size(); ++jPt2 ) {
721  Scalar diff = eucDistance(vBlPtJ_ [ jBl ].vEnt() [ jPt1 ], vBlPtJ_ [ jBl ].vEnt() [ jPt2 ]);
722 
723  if ( ( diff - diag ) >= ( -Constants :: TOLERANCE ) && ( diff - inData_->dr() ) <= ( +Constants :: TOLERANCE ) ) {
724  vBlPtJ_ [ jBl ].vEnt() [ jPt1 ].vpNeighbor1().push_back(& vBlPtJ_ [ jBl ].vEnt() [ jPt2 ]);
725  }
726 
727  if ( vBlPtJ_ [ jBl ].vEnt() [ jPt1 ].vpNeighbor1().size() == 4 ) {
728  break;
729  }
730  }
731  }
732  }
733  }
734 
735  void setOnBlockBoundaryPtQd1_()
736  {
737  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
738  for ( unsigned int iPt = 0; iPt < vBlPtI_ [ iBl ].size(); ++iPt ) {
739  if ( vBlPtI_ [ iBl ].vEnt() [ iPt ].vpNeighbor0().size() < 4 ) {
740  if ( !vBlPtI_ [ iBl ].vEnt() [ iPt ].onBlockBoundary() ) {
741  vBlPtI_ [ iBl ].vEnt() [ iPt ].onBlockBoundary() = true;
742 #if ENABLE_MFILES_OUTPUT == 1
743  vBlPtI_ [ iBl ].vEnt() [ iPt ].marker().size += 4;
744  vBlPtI_ [ iBl ].vEnt() [ iPt ].marker().edgeColor = Blue;
745  vBlPtI_ [ iBl ].vEnt() [ iPt ].marker().faceColor = Red;
746 #endif
747  }
748  }
749  }
750 
751  for ( unsigned int iPt = 0; iPt < vBlPtI_ [ iBl ].size(); ++iPt ) {
752  if ( vBlPtI_ [ iBl ].vEnt() [ iPt ].vpNeighbor1().size() < 4 ) {
753  if ( !vBlPtI_ [ iBl ].vEnt() [ iPt ].onBlockBoundary() ) {
754  vBlPtI_ [ iBl ].vEnt() [ iPt ].onBlockBoundary() = true;
755 #if ENABLE_MFILES_OUTPUT == 1
756  vBlPtI_ [ iBl ].vEnt() [ iPt ].marker().size += 4;
757  vBlPtI_ [ iBl ].vEnt() [ iPt ].marker().edgeColor = Magenta;
758  vBlPtI_ [ iBl ].vEnt() [ iPt ].marker().faceColor = Red;
759 #endif
760  }
761  }
762  }
763 
764  std :: vector< Point * > vPt1;
765  vPt1.clear();
766  for ( unsigned int iPt1 = 0; iPt1 < vBlPtI_ [ iBl ].size(); ++iPt1 ) {
767  if ( vBlPtI_ [ iBl ].vEnt() [ iPt1 ].vpNeighbor1().size() != 4 ) {
768  continue;
769  }
770 
771  /*
772  * Note: this criteria identify undesired points as "Yellow edge points"
773  * we take them for a moment until we use another criteria (cell based criteria)
774  * to identify them and take them out of the boundary "Black edge points"
775  */
776  if ( vBlPtI_ [ iBl ].vEnt() [ iPt1 ].onBlockBoundvpNeighbor1().size() >= 2 ) {
777  if ( !vBlPtI_ [ iBl ].vEnt() [ iPt1 ].onBlockBoundary() ) {
778  vPt1.push_back(& vBlPtI_ [ iBl ].vEnt() [ iPt1 ]);
779  }
780  }
781  }
782 
783  for ( unsigned int iPt1 = 0; iPt1 < vPt1.size(); ++iPt1 ) {
784  if ( !vPt1 [ iPt1 ]->onBlockBoundary() ) {
785  vPt1 [ iPt1 ]->onBlockBoundary() = true;
786 #if ENABLE_MFILES_OUTPUT == 1
787  vPt1 [ iPt1 ]->marker().size += 4;
788  vPt1 [ iPt1 ]->marker().edgeColor = Yellow;
789  vPt1 [ iPt1 ]->marker().faceColor = Red;
790 #endif
791  }
792  }
793  }
794  }
795 
796  void markPtQd1_()
797  {
798  for ( unsigned int jBl = 0; jBl < inData_->blockNumber(); ++jBl ) {
799  for ( unsigned int jPt = 0; jPt < vBlPtJ_ [ jBl ].size(); ++jPt ) {
800  /* check if the considered point has to be marked */
801  if ( vBlPtJ_ [ jBl ].vEnt() [ jPt ].vpNeighbor0().size() < 4 ) {
802  if ( !vBlPtJ_ [ jBl ].vEnt() [ jPt ].mark() ) {
803  vBlPtJ_ [ jBl ].vEnt() [ jPt ].mark() = true;
804 #if ENABLE_MFILES_OUTPUT == 1
805  vBlPtJ_ [ jBl ].vEnt() [ jPt ].marker().size += 4;
806  vBlPtJ_ [ jBl ].vEnt() [ jPt ].marker().faceColor = Cyan;
807 #endif
808  }
809  }
810  }
811  }
812  }
813 
814  void fillvBlEd_()
815  {
816 #if VERBOSE_LEVEL >= 2
817  std :: cout << "-- fill vBlEd_\n";
818 #endif
819  vBlEd_.resize( inData_->blockNumber() );
820  for ( unsigned int jBl = 0; jBl < inData_->blockNumber(); ++jBl ) {
821  vBlEd_ [ jBl ].clear();
822  for ( unsigned int jPt = 0; jPt < vBlPtJ_ [ jBl ].size(); ++jPt ) {
823  if ( vBlPtJ_ [ jBl ].vEnt() [ jPt ].vpNeighbor0().size() != 4 ) {
824  continue;
825  }
826 
827  /* save the considered point and its neighbors in a new vector */
828  std :: vector< Point * > vp0( vBlPtJ_ [ jBl ].vEnt() [ jPt ].vpNeighbor0() );
829  vp0.insert(vp0.begin(), & vBlPtJ_ [ jBl ].vEnt() [ jPt ]); /* insert in the beginning of the vector */
830  const int &nbrNeighbor = vBlPtJ_ [ jBl ].vEnt() [ jPt ].vpNeighbor0().size();
831 
832  for ( unsigned int iPt = 1; iPt < nbrNeighbor; ++iPt ) {
833  for ( unsigned int kPt = iPt + 1; kPt <= nbrNeighbor; ++kPt ) {
834  if ( std :: abs( eucDistance(* vp0 [ iPt ], * vp0 [ kPt ]) - inData_->dr() ) <= ( +Constants :: TOLERANCE ) ) {
835  vBlEd_ [ jBl ].add( Edge(vp0 [ iPt ], vp0 [ kPt ]) );
836  }
837  }
838  }
839  }
840 
841 #if VERBOSE_LEVEL >= 2
842  std :: cout << "vBlEd_[jBl].size() = " << vBlEd_ [ jBl ].size() << "\n";
843 #endif
844  }
845  }
846 
847  void cleanvBlEd_()
848  {
849 #if VERBOSE_LEVEL >= 2
850  std :: cout << "-- clean vBlEd_ : mark then remove duplicated edges from vBlEd_\n";
851  std :: cout << "---- mark duplicated edges in vBlEd_\n";
852 #endif
853  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
854  for ( unsigned int iEd1 = 0; iEd1 < vBlEd_ [ iBl ].size(); ++iEd1 ) {
855  for ( unsigned int iEd2 = ( iEd1 + 1 ); iEd2 < vBlEd_ [ iBl ].size(); ++iEd2 ) {
856  if ( superposed(vBlEd_ [ iBl ].vEnt() [ iEd1 ], vBlEd_ [ iBl ].vEnt() [ iEd2 ]) ) {
857  if ( !vBlEd_ [ iBl ].vEnt() [ iEd2 ].mark() ) {
858  vBlEd_ [ iBl ].vEnt() [ iEd2 ].mark() = true;
859  }
860  }
861  }
862  }
863  }
864 
865 #if VERBOSE_LEVEL >= 2
866  std :: cout << "---- remove duplicated edges from vBlEd_\n";
867 #endif
868  std :: vector< Edge > vec1;
869  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
870  vec1.clear();
871  for ( unsigned int iEd1 = 0; iEd1 < vBlEd_ [ iBl ].size(); ++iEd1 ) {
872  if ( !vBlEd_ [ iBl ].vEnt() [ iEd1 ].mark() ) {
873  vec1.push_back(vBlEd_ [ iBl ].vEnt() [ iEd1 ]);
874  }
875  }
876 
877  vBlEd_ [ iBl ].vEnt().clear();
878  for ( unsigned int iEd1 = 0; iEd1 < vec1.size(); ++iEd1 ) {
879  vBlEd_ [ iBl ].vEnt().push_back(vec1 [ iEd1 ]);
880  }
881 
882 #if VERBOSE_LEVEL >= 2
883  std :: cout << "vBlEd_[iBl].size() = " << vBlEd_ [ iBl ].size() << "\n";
884 #endif
885  }
886  }
887 
888  void fillPointsVEd_()
889  {
890 #if VERBOSE_LEVEL >= 2
891  std :: cout << "-- fill vEd_ (vector of edge instance addresses) of each point\n";
892 #endif
893  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
894  for ( unsigned int iEd1 = 0; iEd1 < vBlEd_ [ iBl ].size(); ++iEd1 ) {
895  vBlEd_ [ iBl ].vEnt() [ iEd1 ].vPt() [ 0 ]->vEd().push_back(& vBlEd_ [ iBl ].vEnt() [ iEd1 ]); // first point vPt()[0]
896  vBlEd_ [ iBl ].vEnt() [ iEd1 ].vPt() [ 1 ]->vEd().push_back(& vBlEd_ [ iBl ].vEnt() [ iEd1 ]); // second point vPt()[1]
897  }
898  }
899  }
900 
901  void fillvBlQd_()
902  {
903 #if VERBOSE_LEVEL >= 2
904  std :: cout << "-- fill vBlQd_\n";
905 #endif
906  vBlQd_.clear();
907  vBlQd_.resize( inData_->blockNumber() );
908  for ( unsigned int jBl1 = 0; jBl1 < inData_->blockNumber(); ++jBl1 ) {
909  vBlQd_ [ jBl1 ].clear();
910  for ( unsigned int jPt1 = 0; jPt1 < vBlPtJ_ [ jBl1 ].size(); ++jPt1 ) {
911  std :: vector< msbGrid :: Edge * > vEd1( neighborsFromvBlEd_( jBl1, inData_->dr() / 2.0, & vBlPtJ_ [ jBl1 ].vEnt() [ jPt1 ] ) );
912  if ( vEd1.size() == 4 ) {
913  Quadrangle qd1(vEd1);
914 
915  /* add this quadrangle to vBlQd_ */
916  vBlQd_ [ jBl1 ].add(qd1);
917  }
918  }
919 
920 #if VERBOSE_LEVEL >= 2
921  std :: cout << "vBlQd_[jBl1].size() = " << vBlQd_ [ jBl1 ].size() << "\n";
922 #endif
923  }
924  }
925 
926  const std :: vector< msbGrid :: Edge * > neighborsFromvBlEd_(const int &iBl0, const Scalar &dist0, const Point *pt0) /*const*/
927  {
928  std :: vector< msbGrid :: Edge * > vEd1;
929  vEd1.clear();
930  for ( unsigned int iEd1 = 0; iEd1 < vBlEd_ [ iBl0 ].size(); ++iEd1 ) { /* TODO look only in the neighborhood */
931  if ( std :: abs(eucDistance(Point( vBlEd_ [ iBl0 ].vEnt() [ iEd1 ].pos() ), * pt0) - dist0) <= Constants :: TOLERANCE ) {
932  vEd1.push_back(& vBlEd_ [ iBl0 ].vEnt() [ iEd1 ]);
933  }
934 
935  if ( vEd1.size() == 4 ) {
936  continue;
937  }
938  }
939 
940  return vEd1;
941  }
942 
943  void fillPointsVCel_()
944  {
945 #if VERBOSE_LEVEL >= 2
946  std :: cout << "-- filling vCel_ (vector of cell instance addresses) of each point\n";
947 #endif
948  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
949  fillPointsVCelQd_(iBl);
950  }
951  }
952 
953  void fillEdgesVCel_()
954  {
955 #if VERBOSE_LEVEL >= 2
956  std :: cout << "-- filling vCel_ (vector of cell instance addresses) of each edge\n";
957 #endif
958  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
959  fillEdgesVCelQd_(iBl);
960  }
961  }
962 
963  void markPtQd2_()
964  {
965  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
966  for ( unsigned int iPt = 0; iPt < vBlPtI_ [ iBl ].size(); ++iPt ) {
967  /* check if the considered point has to be marked */
968  if ( vBlPtI_ [ iBl ].vEnt() [ iPt ].vCel().size() == 0 ) {
969  if ( !vBlPtI_ [ iBl ].vEnt() [ iPt ].mark() ) {
970  vBlPtI_ [ iBl ].vEnt() [ iPt ].mark() = true;
971 #if ENABLE_MFILES_OUTPUT == 1
972  vBlPtI_ [ iBl ].vEnt() [ iPt ].marker().size += 4;
973  vBlPtI_ [ iBl ].vEnt() [ iPt ].marker().edgeColor = Blue;
974  vBlPtI_ [ iBl ].vEnt() [ iPt ].marker().faceColor = Green;
975 #endif
976  }
977  }
978  }
979  }
980  }
981 
982  void fillPointsVCelQd_(const int &iBl0)
983  {
984  for ( unsigned int iQd1 = 0; iQd1 < vBlQd_ [ iBl0 ].size(); ++iQd1 ) {
985  vBlQd_ [ iBl0 ].vEnt() [ iQd1 ].vPt() [ 0 ]->vCel().push_back( & ( vBlQd_ [ iBl0 ].vEnt() [ iQd1 ] ) ); // first point vPt()[0]
986  vBlQd_ [ iBl0 ].vEnt() [ iQd1 ].vPt() [ 1 ]->vCel().push_back( & ( vBlQd_ [ iBl0 ].vEnt() [ iQd1 ] ) ); // second point vPt()[1]
987  vBlQd_ [ iBl0 ].vEnt() [ iQd1 ].vPt() [ 2 ]->vCel().push_back( & ( vBlQd_ [ iBl0 ].vEnt() [ iQd1 ] ) ); // third point vPt()[2]
988  vBlQd_ [ iBl0 ].vEnt() [ iQd1 ].vPt() [ 3 ]->vCel().push_back( & ( vBlQd_ [ iBl0 ].vEnt() [ iQd1 ] ) ); // forth point vPt()[3]
989  }
990  }
991 
992  void fillEdgesVCelQd_(const int &iBl0)
993  {
994  for ( unsigned int iQd1 = 0; iQd1 < vBlQd_ [ iBl0 ].size(); ++iQd1 ) {
995  vBlQd_ [ iBl0 ].vEnt() [ iQd1 ].vEd() [ 0 ]->vCel().push_back( & ( vBlQd_ [ iBl0 ].vEnt() [ iQd1 ] ) ); // first edge vEd()[0]
996  vBlQd_ [ iBl0 ].vEnt() [ iQd1 ].vEd() [ 1 ]->vCel().push_back( & ( vBlQd_ [ iBl0 ].vEnt() [ iQd1 ] ) ); // second edge vEd()[1]
997  vBlQd_ [ iBl0 ].vEnt() [ iQd1 ].vEd() [ 2 ]->vCel().push_back( & ( vBlQd_ [ iBl0 ].vEnt() [ iQd1 ] ) ); // third edge vEd()[2]
998  vBlQd_ [ iBl0 ].vEnt() [ iQd1 ].vEd() [ 3 ]->vCel().push_back( & ( vBlQd_ [ iBl0 ].vEnt() [ iQd1 ] ) ); // forth edge vEd()[3]
999  }
1000  }
1001 
1002  void setOnBlockBoundaryCelQd1_(const Color &color0 = Red)
1003  {
1004  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1005  for ( unsigned int iQd1 = 0; iQd1 < vBlQd_ [ iBl ].size(); ++iQd1 ) {
1006  int blBoundPtNumber = 0;
1007  for ( unsigned int iPt1 = 0; iPt1 < 4; ++iPt1 ) {
1008  if ( vBlQd_ [ iBl ].vEnt() [ iQd1 ].vPt() [ iPt1 ]->onBlockBoundary() ) {
1009  ++blBoundPtNumber;
1010  }
1011  }
1012 
1013  if ( blBoundPtNumber < 2 ) {
1014  continue;
1015  }
1016 
1017  vBlQd_ [ iBl ].vEnt() [ iQd1 ].onBlockBoundary() = true;
1018 #if ENABLE_MFILES_OUTPUT == 1
1019  vBlQd_ [ iBl ].vEnt() [ iQd1 ].color() = color0;
1020 #endif
1021  }
1022 
1023  for ( unsigned int iQd1 = 0; iQd1 < vBlQd_ [ iBl ].size(); ++iQd1 ) {
1024  if ( !vBlQd_ [ iBl ].vEnt() [ iQd1 ].onBlockBoundary() ) {
1025  continue;
1026  }
1027 
1028  for ( unsigned int iQd2 = 0; iQd2 < vBlQd_ [ iBl ].size(); ++iQd2 ) {
1029  if ( iQd2 == iQd1 ) {
1030  continue;
1031  }
1032 
1033  if ( !vBlQd_ [ iBl ].vEnt() [ iQd2 ].onBlockBoundary() ) {
1034  continue;
1035  }
1036 
1037  int commonBoundPtNumber = 0;
1038  for ( unsigned int iPt1 = 0; iPt1 < 4; ++iPt1 ) {
1039  if ( !vBlQd_ [ iBl ].vEnt() [ iQd1 ].vPt() [ iPt1 ]->onBlockBoundary() ) {
1040  continue;
1041  }
1042 
1043  for ( unsigned int iPt2 = 0; iPt2 < 4; ++iPt2 ) {
1044  if ( !vBlQd_ [ iBl ].vEnt() [ iQd2 ].vPt() [ iPt2 ]->onBlockBoundary() ) {
1045  continue;
1046  }
1047 
1048  if ( vBlQd_ [ iBl ].vEnt() [ iQd1 ].vPt() [ iPt1 ] == vBlQd_ [ iBl ].vEnt() [ iQd2 ].vPt() [ iPt2 ] ) {
1049  ++commonBoundPtNumber;
1050  }
1051  }
1052  }
1053 
1054  if ( commonBoundPtNumber != 2 ) {
1055  continue;
1056  }
1057 
1058  int qd1BoundPtNumber = 0;
1059  for ( unsigned int iPt1 = 0; iPt1 < 4; ++iPt1 ) {
1060  if ( vBlQd_ [ iBl ].vEnt() [ iQd1 ].vPt() [ iPt1 ]->onBlockBoundary() ) {
1061  ++qd1BoundPtNumber;
1062  }
1063  }
1064 
1065  if ( qd1BoundPtNumber == 2 ) {
1066  vBlQd_ [ iBl ].vEnt() [ iQd1 ].onBlockBoundary() = false;
1067 #if ENABLE_MFILES_OUTPUT == 1
1068  vBlQd_ [ iBl ].vEnt() [ iQd1 ].color() = Black;
1069 #endif
1070  }
1071 
1072  int qd2BoundPtNumber = 0;
1073  for ( unsigned int iPt2 = 0; iPt2 < 4; ++iPt2 ) {
1074  if ( vBlQd_ [ iBl ].vEnt() [ iQd2 ].vPt() [ iPt2 ]->onBlockBoundary() ) {
1075  ++qd2BoundPtNumber;
1076  }
1077  }
1078 
1079  if ( qd2BoundPtNumber == 2 ) {
1080  vBlQd_ [ iBl ].vEnt() [ iQd2 ].onBlockBoundary() = false;
1081 #if ENABLE_MFILES_OUTPUT == 1
1082  vBlQd_ [ iBl ].vEnt() [ iQd2 ].color() = Black;
1083 #endif
1084  }
1085  }
1086  }
1087 
1088  vvBoundQdPtr_.clear();
1089  vvBoundQdPtr_.resize( inData_->blockNumber() );
1090  for ( unsigned int iQd1 = 0; iQd1 < vBlQd_ [ iBl ].size(); ++iQd1 ) {
1091  if ( vBlQd_ [ iBl ].vEnt() [ iQd1 ].onBlockBoundary() ) {
1092  vvBoundQdPtr_ [ iBl ].push_back(& vBlQd_ [ iBl ].vEnt() [ iQd1 ]);
1093  }
1094  }
1095 
1096 #if VERBOSE_LEVEL >= 2
1097  std :: cout << "vvBoundQdPtr_[iBl].size() = " << vvBoundQdPtr_ [ iBl ].size() << "\n";
1098 #endif
1099  }
1100  }
1101 
1102  void setOnBlockBoundEdQd_(const Color &color0 = Blue)
1103  {
1104  for ( unsigned int iBl1 = 0; iBl1 < inData_->blockNumber(); ++iBl1 ) {
1105  for ( unsigned int iEd1 = 0; iEd1 < vBlEd_ [ iBl1 ].size(); ++iEd1 ) {
1106  if ( vBlEd_ [ iBl1 ].vEnt() [ iEd1 ].vPt() [ 0 ]->onBlockBoundary() &&
1107  vBlEd_ [ iBl1 ].vEnt() [ iEd1 ].vPt() [ 1 ]->onBlockBoundary() ) {
1108  if ( !vBlEd_ [ iBl1 ].vEnt() [ iEd1 ].onBlockBoundary() ) {
1109  vBlEd_ [ iBl1 ].vEnt() [ iEd1 ].onBlockBoundary() = true;
1110 #if ENABLE_MFILES_OUTPUT == 1
1111  vBlEd_ [ iBl1 ].vEnt() [ iEd1 ].color() = color0;
1112 #endif
1113  }
1114  }
1115  }
1116  }
1117 
1118  for ( unsigned int iBl1 = 0; iBl1 < inData_->blockNumber(); ++iBl1 ) {
1119  for ( unsigned int iQd1 = 0; iQd1 < vBlQd_ [ iBl1 ].size(); ++iQd1 ) {
1120  if ( !vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].onBlockBoundary() ) {
1121  continue;
1122  }
1123 
1124  for ( unsigned int iQd2 = 0; iQd2 < vBlQd_ [ iBl1 ].size(); ++iQd2 ) {
1125  if ( iQd2 == iQd1 ) {
1126  continue;
1127  }
1128 
1129  for ( unsigned int iEd1 = 0; iEd1 < 4; ++iEd1 ) {
1130  if ( !vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].vEd() [ iEd1 ]->onBlockBoundary() ) {
1131  continue;
1132  }
1133 
1134  for ( unsigned int iEd2 = 0; iEd2 < 4; ++iEd2 ) {
1135  if ( !vBlQd_ [ iBl1 ].vEnt() [ iQd2 ].vEd() [ iEd2 ]->onBlockBoundary() ) {
1136  continue;
1137  }
1138 
1139  if ( vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].vEd() [ iEd1 ] == vBlQd_ [ iBl1 ].vEnt() [ iQd2 ].vEd() [ iEd2 ] ) {
1140  vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].vEd() [ iEd1 ]->onBlockBoundary() = false;
1141 #if ENABLE_MFILES_OUTPUT == 1
1142  vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].vEd() [ iEd1 ]->color() = Black;
1143 #endif
1144  }
1145  }
1146  }
1147  }
1148  }
1149  }
1150  }
1151 
1152  void fillvvBlBoundEdPtr_()
1153  {
1154 #if VERBOSE_LEVEL >= 2
1155  std :: cout << "-- fill block boundary : vvBlBoundEdPtr_\n";
1156 #endif
1157 
1158  vvBlBoundEdPtr_.clear();
1159  vvBlBoundEdPtr_.resize( inData_->blockNumber() );
1160 #if 0
1161  for ( unsigned int iBl1 = 0; iBl1 < inData_->blockNumber(); ++iBl1 ) {
1162  vvBlBoundEdPtr_ [ iBl1 ].clear();
1163  for ( unsigned int iQd1 = 0; iQd1 < vBlQd_ [ iBl1 ].size(); ++iQd1 ) {
1164  if ( !vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].onBlockBoundary() ) {
1165  continue;
1166  }
1167 
1168  short int onBlockBoundEdNumber = 0;
1169  for ( unsigned int iEd1 = 0; iEd1 < 4; ++iEd1 ) {
1170  if ( vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].vEd() [ iEd1 ]->onBlockBoundary() ) {
1171  ++onBlockBoundEdNumber;
1172  }
1173  }
1174 
1175  if ( onBlockBoundEdNumber != 4 ) {
1176  for ( unsigned int iEd1 = 0; iEd1 < 4; ++iEd1 ) {
1177  if ( vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].vEd() [ iEd1 ]->onBlockBoundary() ) {
1178  vvBlBoundEdPtr_ [ iBl1 ].push_back(vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].vEd() [ iEd1 ]);
1179  }
1180  }
1181  }
1182 
1183  if ( onBlockBoundEdNumber == 4 ) {
1184  for ( unsigned int iEd1 = 0; iEd1 < 4; ++iEd1 ) {
1185  vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].vEd() [ iEd1 ]->mark() = true;
1186  }
1187 
1188  for ( unsigned int iPt1 = 0; iPt1 < 4; ++iPt1 ) {
1189  if ( vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].vPt() [ iPt1 ]->vCel().size() == 1 ) {
1190  vBlQd_ [ iBl1 ].vEnt() [ iQd1 ].vPt() [ iPt1 ]->mark() = true;
1191  }
1192  }
1193  }
1194  }
1195 
1196  #if VERBOSE_LEVEL >= 2
1197  std :: cout << "vvBlBoundEdPtr_[iBl1].size() = " << vvBlBoundEdPtr_ [ iBl1 ].size() << "\n";
1198  #endif
1199  }
1200 
1201 #endif
1202 
1203 #if 1
1204  for ( unsigned int iBl1 = 0; iBl1 < inData_->blockNumber(); ++iBl1 ) {
1205  vvBlBoundEdPtr_ [ iBl1 ].clear();
1206  for ( unsigned int iEd1 = 0; iEd1 < vBlEd_ [ iBl1 ].size(); ++iEd1 ) {
1207  if ( vBlEd_ [ iBl1 ].vEnt() [ iEd1 ].onBlockBoundary() ) {
1208  vvBlBoundEdPtr_ [ iBl1 ].push_back(& vBlEd_ [ iBl1 ].vEnt() [ iEd1 ]);
1209  }
1210  }
1211 
1212  #if VERBOSE_LEVEL >= 2
1213  std :: cout << "vvBlBoundEdPtr_[iBl1].size() = " << vvBlBoundEdPtr_ [ iBl1 ].size() << "\n";
1214  #endif
1215  }
1216 
1217 #endif
1218  }
1219 
1220  void setOnBlockBoundaryPtQd2_()
1221  {
1222  for ( unsigned int iBl1 = 0; iBl1 < inData_->blockNumber(); ++iBl1 ) {
1223  for ( unsigned int iPt1 = 0; iPt1 < vBlPtI_ [ iBl1 ].size(); ++iPt1 ) {
1224  if ( !vBlPtI_ [ iBl1 ].vEnt() [ iPt1 ].onBlockBoundary() ) {
1225  continue;
1226  }
1227 
1228  short int onBlockBoundEdNumber = 0;
1229  for ( unsigned int iEd1 = 0; iEd1 < vBlPtI_ [ iBl1 ].vEnt() [ iPt1 ].vEd().size(); ++iEd1 ) {
1230  if ( vBlPtI_ [ iBl1 ].vEnt() [ iPt1 ].vEd() [ iEd1 ]->onBlockBoundary() ) {
1231  ++onBlockBoundEdNumber;
1232  }
1233  }
1234 
1235  if ( onBlockBoundEdNumber == 0 ) {
1236  vBlPtI_ [ iBl1 ].vEnt() [ iPt1 ].onBlockBoundary() = false;
1237 #if ENABLE_MFILES_OUTPUT == 1
1238  vBlPtI_ [ iBl1 ].vEnt() [ iPt1 ].marker().size -= 4;
1239  vBlPtI_ [ iBl1 ].vEnt() [ iPt1 ].marker().edgeColor = Black;
1240  vBlPtI_ [ iBl1 ].vEnt() [ iPt1 ].marker().faceColor = Red;
1241 #endif
1242  }
1243  }
1244  }
1245 
1246  fillvvBoundPtPtr_();
1247  }
1248 
1249  void fillvvBoundPtPtr_()
1250  {
1251  vvBoundPtPtr_.clear();
1252  vvBoundPtPtr_.resize( inData_->blockNumber() );
1253  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1254  for ( unsigned int iPt = 0; iPt < vBlPtI_ [ iBl ].size(); ++iPt ) {
1255  if ( vBlPtI_ [ iBl ].vEnt() [ iPt ].onBlockBoundary() ) {
1256  vvBoundPtPtr_ [ iBl ].push_back(& vBlPtI_ [ iBl ].vEnt() [ iPt ]);
1257  }
1258  }
1259  }
1260  }
1261 
1262  void setOnBlockBoundaryCelQd2_(const Color &color0 = Black)
1263  {
1264  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1265  for ( unsigned int iQd1 = 0; iQd1 < vBlQd_ [ iBl ].size(); ++iQd1 ) {
1266  short int onBlockBoundEdNumber = 0;
1267  for ( unsigned int iEd1 = 0; iEd1 < 4; ++iEd1 ) {
1268  if ( vBlQd_ [ iBl ].vEnt() [ iQd1 ].vEd() [ iEd1 ]->onBlockBoundary() ) {
1269  ++onBlockBoundEdNumber;
1270  }
1271  }
1272 
1273  if ( onBlockBoundEdNumber == 0 ) {
1274  vBlQd_ [ iBl ].vEnt() [ iQd1 ].onBlockBoundary() = false;
1275 #if ENABLE_MFILES_OUTPUT == 1
1276  vBlQd_ [ iBl ].vEnt() [ iQd1 ].color() = color0;
1277 #endif
1278  }
1279  }
1280  }
1281  }
1282 
1283  void markCelQd1_(const Color &color0 = Black)
1284  {
1285  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1286  if ( vBlQd_ [ iBl ].size() == 0 ) {
1287  std :: cerr << "warning : block number " << iBl << " is empty!!\n";
1288  continue;
1289  }
1290 
1291  if ( vBlQd_ [ iBl ].size() == 1 ) {
1292  std :: cerr << "warning : block number " << iBl << " has only one cell!!\n";
1293  continue;
1294  }
1295 
1296  for ( unsigned int iQd1 = 0; iQd1 < vBlQd_ [ iBl ].size(); ++iQd1 ) {
1297  if ( !vBlQd_ [ iBl ].vEnt() [ iQd1 ].onBlockBoundary() ) {
1298  continue;
1299  }
1300 
1301  bool isolatedCell = true;
1302  for ( unsigned int iEd1 = 0; iEd1 < 4; ++iEd1 ) {
1303  if ( vBlQd_ [ iBl ].vEnt() [ iQd1 ].vEd() [ iEd1 ]->vCel().size() != 1 ) {
1304  isolatedCell = false;
1305  continue;
1306  }
1307  }
1308 
1309  if ( isolatedCell ) {
1310  std :: cerr << "isolated cell in block number " << iBl << "\n";
1311  if ( !vBlQd_ [ iBl ].vEnt() [ iQd1 ].mark() ) {
1312  vBlQd_ [ iBl ].vEnt() [ iQd1 ].mark() = true;
1313 #if ENABLE_MFILES_OUTPUT == 1
1314  vBlQd_ [ iBl ].vEnt() [ iQd1 ].color() = color0;
1315 #endif
1316  }
1317  }
1318  }
1319  }
1320  }
1321 
1322  void setIdentifiers_()
1323  {
1324 #if VERBOSE_LEVEL >= 2
1325  std :: cout << "-- set identifiers\n";
1326  std :: cout << "---- Points identification : blPtExtBound_, vBlPtI_ and vBlPtJ_\n";
1327 #endif
1328  int long idPt1 = 1;
1329  idPt1 = blPtExtBound_.putId(idPt1);
1330  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1331  vBlPtI_ [ iBl ].Id() = iBl + 1;
1332  idPt1 = vBlPtI_ [ iBl ].putId(idPt1);
1333  }
1334 
1335 #if VERBOSE_LEVEL >= 1
1336  std :: cout << "totPtNumber_ = " << ( idPt1 - 1 ) << "\n";
1337 #endif
1338 
1339 #if VERBOSE_LEVEL >= 2
1340  std :: cout << "---- Edges identification : blEdExtBound_ and vBlEd_\n";
1341 #endif
1342  int long idEd1 = 1;
1343  idEd1 = blEdExtBound_.putId(idEd1);
1344  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1345  vBlEd_ [ iBl ].Id() = iBl + 1;
1346  idEd1 = vBlEd_ [ iBl ].putId(idEd1);
1347  }
1348 
1349 #if VERBOSE_LEVEL >= 1
1350  std :: cout << "totEdNumber_ = " << ( idEd1 - 1 ) << "\n";
1351 #endif
1352 
1353 #if VERBOSE_LEVEL >= 2
1354  std :: cout << "---- Quadrangles identification : vBlQd_\n";
1355 #endif
1356 
1357  const int firstQdId = inData_->blockNumber() + 1 + 2 /*one for IG and 1 for EHM*/;
1358  int long idQd1 = firstQdId;
1359  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1360  vBlQd_ [ iBl ].Id() = iBl + 1;
1361  idQd1 = vBlQd_ [ iBl ].putId(idQd1);
1362  }
1363 
1364  // totQdNumber_ = idQd1 - 1;
1365 
1366 #if VERBOSE_LEVEL >= 1
1367  std :: cout << "totQdNumber_ = " << totQdNumber_ << "\n";
1368 #endif
1369  }
1370 
1371  void outputMFiles_() // TODO revise this method.
1372  {
1373 #if OUTPUT_MFILES_DIR_CREATED != 1
1374  const std :: string createDirectoryCommandString = "mkdir -p " + inData_->outputPath() + "m/" + outputDirName_;
1375  char *createDirectoryCommandChar = ( char * ) createDirectoryCommandString.c_str();
1376  system(createDirectoryCommandChar);
1377 #endif
1378 
1379  std :: cout << "-- write *.m output files\n";
1380  /* matlab writer instances */
1381  MFilesWriter mwAllEntities(inData_),
1382  mwvBlPtI_(inData_), mwvBlPtJ_(inData_), mwblPtExtBound_(inData_),
1383  mwvvBoundPtPtr_(inData_), mwvvBlBoundEdPtr_(inData_);
1384  MFilesWriter mwvvBoundQdPtr_(inData_), mwvBlQd_(inData_);
1385 
1386  /* open matlab output files */
1387  mwAllEntities.openFile(outputDirName_ + "/allEntities");
1388  mwvBlPtI_.openFile(outputDirName_ + "/vBlPtI_");
1389  mwvBlPtJ_.openFile(outputDirName_ + "/vBlPtJ_");
1390  mwblPtExtBound_.openFile(outputDirName_ + "/blPtExtBound_");
1391  mwvvBoundPtPtr_.openFile(outputDirName_ + "/vvBoundPtPtr_");
1392  mwvvBlBoundEdPtr_.openFile(outputDirName_ + "/vvBlBoundEdPtr_");
1393 
1394  mwvvBoundQdPtr_.openFile(outputDirName_ + "/vvBoundQdPtr_");
1395  mwvBlQd_.openFile(outputDirName_ + "/vBlQd_");
1396 
1397  /* handle points */
1398  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1399  vBlPtI_ [ iBl ].putMarkerSymbol("'o'");
1400  mwvBlPtI_.write(vBlPtI_ [ iBl ], false);
1401  mwAllEntities.write(vBlPtI_ [ iBl ], false);
1402 
1403  vBlPtJ_ [ iBl ].putMarkerSymbol("'s'");
1404  mwvBlPtJ_.write(vBlPtJ_ [ iBl ], false);
1405  mwAllEntities.write(vBlPtJ_ [ iBl ], false);
1406  }
1407 
1408  blPtExtBound_.putMarkerSymbol("'o'");
1409  mwblPtExtBound_.write(blPtExtBound_, false);
1410  mwAllEntities.write(blPtExtBound_, false);
1411  mwvBlPtI_.closeFile();
1412  mwvBlPtJ_.closeFile();
1413  mwblPtExtBound_.closeFile();
1414 
1415  /* handle quadrangles */
1416  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1417  mwvBlQd_.write(vBlQd_ [ iBl ], false);
1418  mwAllEntities.write(vBlQd_ [ iBl ], false);
1419  }
1420 
1421  mwvBlQd_.closeFile();
1422  mwAllEntities.closeFile();
1423 
1424  /* handle boundaries */
1425  for ( unsigned int iBl = 0; iBl < vvBoundPtPtr_.size(); ++iBl ) {
1426  for ( unsigned int i = 0; i < vvBoundPtPtr_ [ iBl ].size(); ++i ) {
1427  mwvvBoundPtPtr_.write(* vvBoundPtPtr_ [ iBl ] [ i ], true);
1428  }
1429  }
1430 
1431  mwvvBoundPtPtr_.closeFile();
1432 
1433  /* output points in the same file as Edge instances */
1434  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1435  mwvvBlBoundEdPtr_.write(vBlPtI_ [ iBl ], false);
1436  }
1437 
1438  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1439  mwvvBlBoundEdPtr_.write(vBlPtJ_ [ iBl ], false);
1440  }
1441 
1442  /* output Edge instances */
1443  for ( unsigned int iBl = 0; iBl < ( inData_->blockNumber() - 1 ); ++iBl ) {
1444  for ( unsigned int i = 0; i < vvBlBoundEdPtr_ [ iBl ].size(); ++i ) {
1445  mwvvBlBoundEdPtr_.write(* vvBlBoundEdPtr_ [ iBl ] [ i ], false);
1446  }
1447  }
1448 
1449  mwvvBlBoundEdPtr_.closeFile();
1450 
1451  /* output points in the same file as Quadrangle instances */
1452  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1453  mwvvBoundQdPtr_.write(vBlPtI_ [ iBl ], false);
1454  }
1455 
1456  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1457  mwvvBoundQdPtr_.write(vBlPtJ_ [ iBl ], false);
1458  }
1459 
1460  /* output Quadrangle instances */
1461  for ( unsigned int iBl = 0; iBl < vvBoundQdPtr_.size(); ++iBl ) {
1462  for ( unsigned int i = 0; i < vvBoundQdPtr_ [ iBl ].size(); ++i ) {
1463  mwvvBoundQdPtr_.write(* vvBoundQdPtr_ [ iBl ] [ i ], false);
1464  }
1465  }
1466 
1467  mwvvBoundQdPtr_.closeFile();
1468  }
1469 
1470  void outputGeoFiles_()
1471  {
1472  const std :: string createDirectoryCommandString = "mkdir -p " + inData_->outputPath() + "geo/" + outputDirName_;
1473  char *createDirectoryCommandChar = ( char * ) createDirectoryCommandString.c_str();
1474  system(createDirectoryCommandChar);
1475  #define OUTPUT_GEOFILES_DIR_CREATED 1
1476 
1477 #if WRITE_MSBGRID_GEOFILES == 1
1478  /* msbGrid output */
1479  outputMsbGridGeoFiles_();
1480 #endif
1481 
1482 #if WRITE_FULLY_UNSTRUCTURED_GEOFILE == 1
1483  /* fullyUnstructured output */
1484  outputFullyUnstructuredGeoFile_();
1485 #endif
1486  }
1487 
1488  void outputMsbGridGeoFiles_()
1489  {
1490 #if VERBOSE_LEVEL >= 1
1491  std :: cout << "-- write msbGrid \"*.geo\" files\n";
1492 #endif
1493 
1494  GeoFilesWriter gwblPtExtBound_(inData_), gwvBlPtI_(inData_);
1495  GeoFilesWriter gwvEdExtBound_(inData_), gwvBlEd_(inData_);
1496  GeoFilesWriter gwvBlQd_(inData_);
1497  GeoFilesWriter gwBlockBoundary(inData_);
1498 
1499  gwblPtExtBound_.openFile(outputDirName_ + "/blPtExtBound_");
1500  gwvBlPtI_.openFile(outputDirName_ + "/vBlPtI_");
1501  gwvEdExtBound_.openFile(outputDirName_ + "/vEdExtBound_");
1502  gwvBlEd_.openFile(outputDirName_ + "/vBlEd_");
1503  gwvBlQd_.openFile(outputDirName_ + "/vBlQd_");
1504  gwBlockBoundary.openFile(outputDirName_ + "/blockBoundary_");
1505 
1506  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1507  gwvBlPtI_.write(vBlPtI_ [ iBl ]);
1508  gwvBlEd_.write(vBlEd_ [ iBl ]);
1509  gwvBlQd_.write(vBlQd_ [ iBl ], /*physical=*/ true);
1510  }
1511 
1512  std :: string firstSurfaceIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << vBlQd_.front().vEnt().front().Id() /*1*/ ) )->str();
1513  std :: string lastSurfaceIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << vBlQd_.back().vEnt().back().Id() /*totQdNumber_*/ ) )->str();
1514  gwvBlQd_.write("\nColor White {Surface {" + firstSurfaceIdString + ":" + lastSurfaceIdString + "};}\n");
1515 
1516  gwvEdExtBound_.write("\n//////////////////////////////\n");
1517  gwvEdExtBound_.write("// external boundary lines //\n");
1518  gwvEdExtBound_.write("//////////////////////////////\n");
1519  for ( unsigned int iEd = 0; iEd < vEdExtBound_.size(); ++iEd ) {
1520  gwvEdExtBound_.write(* vEdExtBound_ [ iEd ]);
1521  }
1522 
1523  gwblPtExtBound_.write("\n//////////////////////////////\n");
1524  gwblPtExtBound_.write("// external boundary points //\n");
1525  gwblPtExtBound_.write("//////////////////////////////\n");
1526  blPtExtBound_.putCl( inData_->dr() );
1527  gwblPtExtBound_.write(blPtExtBound_);
1528 
1529  gwBlockBoundary.write("\n///////////////////////\n");
1530  gwBlockBoundary.write("// block line Loops //\n");
1531  gwBlockBoundary.write("///////////////////////\n");
1532  const int blockBoundaryIdMin = vBlQd_.back().vEnt().back().Id() + 1;
1533  int blockBoundaryId = blockBoundaryIdMin;
1534  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1535  std :: string blockIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << vBlQd_ [ iBl ].Id() ) )->str();
1536  std :: string blockBoundaryIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << blockBoundaryId ) )->str();
1537  gwBlockBoundary.write("\n/// Boundary of Block " + blockIdString + ", Line Loop " + blockBoundaryIdString + "\n");
1538  gwBlockBoundary.writeLineLoop(vvBlBoundEdPtr_ [ iBl ], blockBoundaryId);
1539  ++blockBoundaryId;
1540  }
1541 
1542  const int blockBoundaryIdMax = blockBoundaryId - 1;
1543 
1544  gwBlockBoundary.write("\n//////////////////////////////////\n");
1545  gwBlockBoundary.write("// external boundary line Loop //\n");
1546  gwBlockBoundary.write("//////////////////////////////////\n");
1547  const int extBoundaryId = 1;
1548  std :: string extBoundaryIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << extBoundaryId ) )->str();
1549  gwBlockBoundary.write("\n/// Line Loop " + extBoundaryIdString + "\n");
1550  gwBlockBoundary.writeLineLoop(vEdExtBound_, extBoundaryId);
1551 
1552  gwBlockBoundary.write("\n/////////////////////////////////\n");
1553  gwBlockBoundary.write("// inter-blocks plane surface //\n");
1554  gwBlockBoundary.write("/////////////////////////////////\n");
1555  const int ig_Physical_Id = vBlQd_.back().Id() + 1;
1556  const int ig_Plane_Id = vBlQd_.back().vEnt().back().Id() + 1;
1557  std :: string ig_Physical_IdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << ig_Physical_Id ) )->str();
1558  std :: string ig_Plane_IdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << ig_Plane_Id ) )->str();
1559  std :: string blockBoundaryIdMinString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << blockBoundaryIdMin ) )->str();
1560  std :: string blockBoundaryIdMaxString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << blockBoundaryIdMax ) )->str();
1561  gwBlockBoundary.write("Plane Surface(" + ig_Plane_IdString + ") = {" + extBoundaryIdString + "," + blockBoundaryIdMinString + ":" + blockBoundaryIdMaxString + "};\n");
1562  gwBlockBoundary.write("Physical Surface(" + ig_Physical_IdString + ") = {" + ig_Plane_IdString + "};\n");
1563  gwBlockBoundary.write("\nColor Red {Surface {" + ig_Plane_IdString + "};}\n");
1564 
1565  gwvBlPtI_.closeFile();
1566  gwvBlEd_.closeFile();
1567  gwvEdExtBound_.closeFile();
1568  gwblPtExtBound_.closeFile();
1569  gwBlockBoundary.closeFile();
1570  gwvBlQd_.closeFile();
1571 
1572  GeoFilesWriter gwmsbGrid_(inData_);
1573  const std :: string msbGridGeoFileName = outputDirName_ + "/msbGrid_" + outputDirName_;
1574  gwmsbGrid_.openFile(msbGridGeoFileName);
1575 
1576  std :: string drString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->dr() ) )->str();
1577  std :: string xminString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->xMin() ) )->str();
1578  std :: string xmaxString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->xMax() ) )->str();
1579  std :: string yminString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->yMin( inData_->xMin() ) ) )->str();
1580  std :: string ymaxString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->yMax( inData_->xMax() ) ) )->str();
1581  std :: string zString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << 0. ) )->str();
1582  std :: string blockNumberString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->blockNumber() ) )->str();
1583  std :: string ehmCString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << 4. ) )->str();
1584 
1585  gwmsbGrid_.write("dr = " + drString + ";\n");
1586  gwmsbGrid_.write("xmin = " + xminString + ";\n");
1587  gwmsbGrid_.write("xmax = " + xmaxString + ";\n");
1588  gwmsbGrid_.write("ymin = " + yminString + ";\n");
1589  gwmsbGrid_.write("ymax = " + ymaxString + ";\n");
1590  gwmsbGrid_.write("z = " + zString + ";\n");
1591  gwmsbGrid_.write("blocks_number = " + blockNumberString + ";\n");
1592  gwmsbGrid_.write("ehm_c = " + ehmCString + ";\n\n");
1593 
1594  gwmsbGrid_.write("Include \"" + gwmsbGrid_.geoFile("blPtExtBound_") + "\";\n");
1595  gwmsbGrid_.write("Include \"" + gwmsbGrid_.geoFile("vBlPtI_") + "\";\n");
1596 
1597  gwmsbGrid_.write("Include \"" + gwmsbGrid_.geoFile("vEdExtBound_") + "\";\n");
1598  gwmsbGrid_.write("Include \"" + gwmsbGrid_.geoFile("vBlEd_") + "\";\n");
1599 
1600  gwmsbGrid_.write("Include \"" + gwmsbGrid_.geoFile("vBlQd_") + "\";\n");
1601 
1602  gwmsbGrid_.write("Include \"" + gwmsbGrid_.geoFile("blockBoundary_") + "\";\n");
1603 #if ADD_RECOMBINEALL == 1
1604  gwmsbGrid_.write("\nInclude \"../../../ehm.geo\";\n");
1605  gwmsbGrid_.write("\nMesh.RecombineAll = 1;\n");
1606 #endif
1607  gwmsbGrid_.closeFile();
1608  }
1609 
1610  void outputFullyUnstructuredGeoFile_() // TODO revise this method.
1611  {
1612 #if VERBOSE_LEVEL >= 1
1613  std :: cout << "-- write fullyUnstructured.geo\n";
1614 #endif
1615  GeoFilesWriter fullyUnstructured(inData_);
1616  const std :: string fullyUnstructuredGeoFileName = outputDirName_ + "/fullyUnstructured_" + outputDirName_;
1617  fullyUnstructured.openFile(fullyUnstructuredGeoFileName);
1618  const std :: string cl0( newCl( "cl0", ( inData_->dr() / inData_->rafCoeff() ) ) );
1619  fullyUnstructured.write(cl0);
1620  fullyUnstructured.write("//////////////////\n");
1621  fullyUnstructured.write("// block points //\n");
1622  fullyUnstructured.write("//////////////////\n");
1623  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1624  std :: string blockIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << ( iBl + 1 ) ) )->str();
1625  fullyUnstructured.write("\n/// Block " + blockIdString + "\n");
1626  for ( unsigned int i = 0; i < vvBoundPtPtr_ [ iBl ].size(); ++i ) {
1627  fullyUnstructured.write(* vvBoundPtPtr_ [ iBl ] [ i ], "cl0");
1628  }
1629  }
1630 
1631  fullyUnstructured.write("\n//////////////////////////////\n");
1632  fullyUnstructured.write("// external boundary points //\n");
1633  fullyUnstructured.write("//////////////////////////////\n");
1634  fullyUnstructured.write(blPtExtBound_, "cl0");
1635 
1636  fullyUnstructured.write("\n//////////////////\n");
1637  fullyUnstructured.write("// block lines //\n");
1638  fullyUnstructured.write("//////////////////\n");
1639 
1640  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1641  std :: string blockIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << ( iBl + 1 ) ) )->str();
1642  fullyUnstructured.write("\n/// Block " + blockIdString + "\n");
1643  for ( unsigned int i = 0; i < vvBlBoundEdPtr_ [ iBl ].size(); ++i ) {
1644  fullyUnstructured.write(* vvBlBoundEdPtr_ [ iBl ] [ i ]);
1645  }
1646  }
1647 
1648  fullyUnstructured.write("\n//////////////////////////////\n");
1649  fullyUnstructured.write("// external boundary lines //\n");
1650  fullyUnstructured.write("//////////////////////////////\n");
1651 
1652  for ( unsigned int iEd = 0; iEd < vEdExtBound_.size(); ++iEd ) {
1653  fullyUnstructured.write(* vEdExtBound_ [ iEd ]);
1654  }
1655 
1656  fullyUnstructured.write("\n///////////////////////\n");
1657  fullyUnstructured.write("// block line Loops //\n");
1658  fullyUnstructured.write("///////////////////////\n");
1659 
1660  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1661  std :: string blockIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << ( iBl + 1 ) ) )->str();
1662  fullyUnstructured.write("\n/// Block " + blockIdString + "\n");
1663  fullyUnstructured.writeLineLoop(vvBlBoundEdPtr_ [ iBl ], iBl + 1);
1664  }
1665 
1666  fullyUnstructured.write("\n//////////////////////////////////\n");
1667  fullyUnstructured.write("// external boundary line Loop //\n");
1668  fullyUnstructured.write("//////////////////////////////////\n");
1669  fullyUnstructured.writeLineLoop(vEdExtBound_, inData_->blockNumber() + 1);
1670 
1671  fullyUnstructured.write("\n///////////////////////////\n");
1672  fullyUnstructured.write("// block plane surfaces //\n");
1673  fullyUnstructured.write("///////////////////////////\n");
1674  const std :: string cl1( newCl( "cl1", ( inData_->dr() / inData_->rafCoeff() ) /*Scalar(1.0)*/ ) );
1675  fullyUnstructured.write(cl1);
1676 
1677  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1678  if ( vvBoundPtPtr_ [ iBl ].size() < 1 ) {
1679  continue;
1680  }
1681 
1682  Point blockBarycenter;
1683  for ( unsigned int iPt1 = 0; iPt1 < vvBoundPtPtr_ [ iBl ].size(); ++iPt1 ) {
1684  blockBarycenter.x() [ 0 ] += vvBoundPtPtr_ [ iBl ] [ iPt1 ]->x() [ 0 ];
1685  blockBarycenter.x() [ 1 ] += vvBoundPtPtr_ [ iBl ] [ iPt1 ]->x() [ 1 ];
1686  }
1687 
1688  blockBarycenter.x() [ 0 ] /= vvBoundPtPtr_ [ iBl ].size();
1689  blockBarycenter.x() [ 1 ] /= vvBoundPtPtr_ [ iBl ].size();
1690  // pBlockId = newp; Point(pBlockId) = {xbc, ybc, 0.0, cl1};
1691  std :: string blockIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << ( iBl + 1 ) ) )->str();
1692  std :: stringstream newPoint;
1693  newPoint << "p" + blockIdString + " = newp; Point(p" + blockIdString + ") = {" << +blockBarycenter.x() [ 0 ] << ", " << blockBarycenter.x() [ 1 ] << ", " << Scalar(0.0) << ", cl1};\n";
1694  fullyUnstructured.write( newPoint.str() );
1695 
1696  std :: stringstream ss;
1697  ss << "Plane Surface(" << blockIdString << ") = {" << blockIdString << "};\n";
1698  // Point{newp} In Surface{BlockId};
1699  ss << "Point{" << "p" + blockIdString << "} In Surface{" << blockIdString << "};\n";
1700 
1701  // ss << "Physical Surface(" + blockIdString + ") = {" + blockIdString + "};\n\n";
1702 
1703  fullyUnstructured.write( ss.str() );
1704  }
1705 
1706  std :: string firstSurfaceIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << 1 ) )->str();
1707  std :: string lastSurfaceIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << inData_->blockNumber() ) )->str();
1708  fullyUnstructured.write("\nPhysical Surface(1) = {" + firstSurfaceIdString + ":" + lastSurfaceIdString + "};\n");
1709  fullyUnstructured.write("\nColor White {Surface {" + firstSurfaceIdString + ":" + lastSurfaceIdString + "};}\n");
1710 
1711  fullyUnstructured.write("/////////////////////////////////\n");
1712  fullyUnstructured.write("// inter-blocks plane surface //\n");
1713  fullyUnstructured.write("/////////////////////////////////\n");
1714  const int extBoundaryId = inData_->blockNumber() + 1;
1715  std :: string extBoundaryIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << extBoundaryId ) )->str();
1716  std :: stringstream ss;
1717  ss << "Plane Surface(" << extBoundaryId << ") = {";
1718  fullyUnstructured.write( ss.str() );
1719  for ( unsigned int iBl = 0; iBl < inData_->blockNumber(); ++iBl ) {
1720  std :: string blockIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << ( iBl + 1 ) ) )->str();
1721  fullyUnstructured.write(blockIdString);
1722  fullyUnstructured.write(",");
1723  }
1724 
1725  std :: string blockIdString = static_cast< std :: ostringstream * >( & ( std :: ostringstream() << ( inData_->blockNumber() + 1 ) ) )->str();
1726  fullyUnstructured.write(blockIdString);
1727  fullyUnstructured.write("};\n");
1728 
1729  std :: stringstream ss0;
1730  ss0 << "Physical Surface(" << blockIdString << ") = {" << blockIdString << "};\n";
1731  ss0 << "Color Red {Surface {" + blockIdString + "};}\n\n";
1732  fullyUnstructured.write( ss0.str() );
1733 #if ADD_RECOMBINEALL == 1
1734  fullyUnstructured.write("\nMesh.RecombineAll = 1;\n");
1735 #endif
1736 
1737  fullyUnstructured.write("// ## END OF FILE ##\n");
1738  fullyUnstructured.closeFile();
1739  }
1740 };
1741 
1742 const std :: string msbGrid :: newCl(const std :: string &clName, const Scalar &clValue)
1743 {
1744  std :: stringstream scl;
1745  scl << "\n// Characteristic length, " << clName << std :: endl;
1746  scl << std :: setw(1) << std :: fixed << clName << " = " << clValue << ";\n\n";
1747  return scl.str();
1748 }
1749 
1750 #endif /* #ifndef _GRID_FACTORY_HH_ */