ParMooN
 All Classes Functions Variables Friends Pages
TeraHertzBreast.h
1 // ==========================================================================
2 // instationary problem
3 // ==========================================================================
4 
5 //===========================================================================
6 // example file
7 // =========================================================================
8 // exact solution in unit cube
9 void ExampleFile()
10 {
11 #define __ROBINBC__
12  OutPut("Example: TeraHertzBrain.h" << endl);
13 }
14 
15 void GridBoundCondition(double x, double y, double z, BoundCond &cond)
16 {
17  cond = DIRICHLET;
18 }
19 // exact solution
20 void Exact(double x, double y, double z, double *values)
21 {
22  values[0] = 0;
23  values[1] = 0;
24  values[2] = 0;
25  values[3] = 0;
26  values[4] = 0;
27 }
28 
29 // initial conditon
30 void InitialCondition(double x, double y, double z, double *values)
31 {
32  values[0] = TDatabase::ParamDB->P0;
33  values[0] = x;
34 }
35 
36 // kind of boundary condition
37 void BoundCondition(double x, double y, double z, BoundCond &cond)
38 {
39 
40 
41  if(TDatabase::ParamDB->P13==0)
42  { //Lside
43  if(x<225 && y>76 )
44  { cond = DIRICHLET; }
45  else
46  { cond = NEUMANN; }
47  }
48  else
49  {
50  if(x< 240 && y<-72. && z>191)
51  { cond = DIRICHLET; }
52  else
53  { cond = NEUMANN; }
54  }
55 
56 
57 
58 
59 }
60 
61 // value of boundary condition
62 void BoundValue(double x, double y, double z, double &value)
63 {
64  value = 0;
65 }
66 
67 void BilinearCoeffs(int n_points, double *X, double *Y, double *Z, double **parameters, double **coeffs)
68 {
69  double eps=0;
70  int i, Region_ID;
71  double *coeff;
72  double x, y, z, hK;
73 
74 
75  double rhsfact, alpha, char_L, beam_r, DimlessBeam_r;
76  double xp=0., yp=0., zp=0., SourceCoord, Sourceradius;
77 
78  if(TDatabase::ParamDB->P13==0)
79  { // center point is (291.977 -1.8733333333333333 159.89 ), old, wrong
80  // 2.922130000000000e+02 -1.8733333333333333e+00 1.583800000000000e+02 1
81  xp=2.922130000000000e+02;
82  zp=1.583800000000000e+02;
83  }
84  else
85  { // side 300.565 69.66666666666666 157.038
86  xp= 300.565;
87  zp=157.038;
88  }
89 
90 
91  beam_r = TDatabase::ParamDB->P11;
92  char_L = TDatabase::ParamDB->P12;
93 
94  hK = coeffs[0][2];
95  Region_ID = coeffs[0][1];
96  DimlessBeam_r = beam_r/char_L;
97 
98 
99  // cout << "Region_ID = " << Region_ID << endl;
100 
101  if(Region_ID==0) //0-Fat/adipost
102  {
103  eps = TDatabase::ParamDB->P1/(1100.*4483*char_L*char_L); ;
104  alpha = TDatabase::ParamDB->P7;
105  rhsfact = (alpha*TDatabase::ParamDB->P4)/(1100.0*4483*(22./7.)*beam_r*beam_r);
106  }
107  else
108  {
109  eps = TDatabase::ParamDB->P2/(1040.0*3500.*char_L*char_L);
110  alpha = TDatabase::ParamDB->P8;
111  rhsfact = alpha*TDatabase::ParamDB->P4/(1040.0*3500.*(22./7.)*beam_r*beam_r);
112  }
113 
114 // //2-CSF
115 // eps = TDatabase::ParamDB->P2/(997.0*3710*char_L*char_L);
116 // alpha = TDatabase::ParamDB->P8;
117 // rhsfact = alpha*TDatabase::ParamDB->P4/(997.0*3710*(22./7.)*beam_r*beam_r);
118 //
119  for(i=0;i<n_points;i++)
120  {
121  coeff = coeffs[i];
122 
123  x = X[i];
124  y = Y[i];
125  z = Z[i];
126 
127  // diffusion
128  coeff[0] = eps;
129  // convection in x direction
130  coeff[1] = 0;
131  // convection in y direction
132  coeff[2] = 0;
133  // convection in z direction
134  coeff[3] = 0;
135  // reaction
136  coeff[4] = 0;
137  // rhs
138 
139  if(TDatabase::ParamDB->P13==0)
140  { //Lside
141  Sourceradius = sqrt( (x-xp)*(x-xp) + (z-zp)*(z-zp));
142 
143  if(y<0.) y=0.;
144  SourceCoord = -y;
145  }
146  else
147  { // side
148  Sourceradius = sqrt( (x-xp)*(x-xp) + (z-zp)*(z-zp));
149 
150  if(y>0.) y=0.;
151  SourceCoord = y;
152  }
153 
154 
155  if(Sourceradius<=DimlessBeam_r)
156  {
157  coeff[5] = rhsfact*exp(alpha*SourceCoord*char_L);// f
158 
159 // cout << "Sourceradius : " << Sourceradius << endl;
160  if(TDatabase::ParamDB->P14< hK) TDatabase::ParamDB->P14=hK;
161  }
162  else
163  {coeff[5] = 0; }// f
164 
165  coeff[6] = 0;
166 
167  }
168 }
169 
170 void ReadAdditionalNModes(char *NODE, tetgenio &InAddPts)
171 {
172  int i, N_Vertices, index, N_Atribude, BDMarker;
173  char line[100];
174  std::ifstream dat(NODE);
175 
176  if (!dat)
177  {
178  cerr << "cannot open '" << NODE << "' for input" << endl;
179  exit(-1);
180  }
181 
182  dat.getline (line, 99); // first line is command line
183  dat >> N_Vertices >> N_Atribude >> BDMarker;
184  dat.getline (line, 99);
185 
186  InAddPts.numberofpoints = N_Vertices;
187  InAddPts.pointlist = new double[3*N_Vertices];
188 
189  //read from file
190  for(i=0;i<N_Vertices; i++)
191  {
192  dat.getline (line, 99);
193  dat >> index >> InAddPts.pointlist[3*i] >> InAddPts.pointlist[3*i+1] >> InAddPts.pointlist[3*i+2];
194  cout<< i << " vert X: " <<InAddPts.pointlist[3*i] << " vert Y: " <<InAddPts.pointlist[3*i+1] <<" vert Z: " <<InAddPts.pointlist[3*i+2] <<endl;
195 
196  }
197 
198  dat.close();
199 
200 // cout << "N_Vertices " <<N_Vertices <<endl;
201 //
202 // exit(0);
203 
204 } // ReadAdditionalNModes
205 
206 
207 
208 void ReadMeditMeshWithCells(char *SMESH, tetgenio &In)
209 {
210  int i, j, k, dimension, N_FVert, N_Faces, N_Vertices, N_CVert, N_Cells;
211  int BDComp_Min=10000000;
212  int *plist;
213 
214  double x;
215 
216  char line[100];
217  tetgenio::facet *F;
218  tetgenio::polygon *P;
219 
220  std::ifstream dat(SMESH);
221 
222  if (!dat)
223  {
224  cerr << "cannot open '" << SMESH << "' for input" << endl;
225  exit(-1);
226  }
227 
228  // check the dimension
229  while (!dat.eof())
230  {
231  dat >> line;
232 
233  if ( (!strcmp(line, "Dimension")) || (!strcmp(line, "dimension")) || (!strcmp(line, "DIMENSION")))
234  {
235  dat.getline (line, 99);
236  dat >> dimension;
237  break;
238  }
239  // read until end of line
240  dat.getline (line, 99);
241  }
242 
243  if(dimension!=3)
244  {
245  cerr << "dimension: " << dimension << endl;
246  exit(-1);
247  }
248 
249  // find the N_Vertices
250  while (!dat.eof())
251  {
252  dat >> line;
253 
254  if ( (!strcmp(line, "Vertices")) || (!strcmp(line, "vertices")) || (!strcmp(line, "VERTICES")) )
255  {
256  dat.getline (line, 99);
257  dat >> N_Vertices;
258  break;
259  }
260  // read until end of line
261  dat.getline (line, 99);
262  }
263 
264  In.numberofpoints = N_Vertices;
265  In.pointlist = new double[3*N_Vertices];
266 
267  //read from file
268  for(i=0;i<N_Vertices; i++)
269  {
270  dat.getline (line, 99);
271  dat >> x >> In.pointlist[3*i+1] >> In.pointlist[3*i+2];
272 
273  In.pointlist[3*i] = x - 3.606596999999999830777142; //center point
274 
275 // if(i==1063 || i==1061 || i==1062 || i==162 || i==1175 || i==517 || i==1064 || i==518 ) // i==1063 as center
276  if(i==1062 || i==916 || i==341 || i==914 || i==162 || i==1063 || i==1061) // i==1062 as center
277  cout<< i << " vert X: " <<In.pointlist[3*i] << " vert Y: " <<In.pointlist[3*i+1] <<" vert Z: " <<In.pointlist[3*i+2] <<endl;
278  }
279 
280  // find the N_Triangles
281  // find the N_Vertices
282  while (!dat.eof())
283  {
284  dat >> line;
285 
286  if ( (!strcmp(line, "Triangles")) || (!strcmp(line, "triangles")) || (!strcmp(line, "TRIANGLES")) )
287  {
288  N_FVert = 3;
289  dat.getline (line, 99);
290  dat >> N_Faces;
291  break;
292  }
293  else if ( (!strcmp(line, "Quadrilaterals")) || (!strcmp(line, "quadrilaterals")) || (!strcmp(line, "QUADRILATERALS")) )
294  {
295  N_FVert = 4;
296  dat.getline (line, 99);
297  dat >> N_Faces;
298  break;
299  }
300 
301  // read until end of line
302  dat.getline (line, 99);
303  }
304 
305  In.numberoffacets = N_Faces;
306  In.facetlist = new tetgenio::facet[In.numberoffacets];
307  In.facetmarkerlist = new int[In.numberoffacets];
308 
309  for(i=0;i<N_Faces; i++)
310  {
311  dat.getline (line, 99);
312 
313  F = &In.facetlist[i];
314  tetgenio::init(F);
315  F->numberofpolygons = 1;
316  F->polygonlist = new tetgenio::polygon[F->numberofpolygons];
317  F->numberofholes = 0;
318  F->holelist = NULL;
319  P = &F->polygonlist[0];
320  tetgenio::init(P);
321  P->numberofvertices = N_FVert;
322  P->vertexlist = new int[P->numberofvertices];
323 
324  for(j=0;j<N_FVert;j++)
325  {
326  dat >> k;
327  P->vertexlist[j] = k-1; // c numbering
328  }
329 
330  dat >> In.facetmarkerlist[i];
331 
332  if(BDComp_Min > In.facetmarkerlist[i])
333  BDComp_Min=In.facetmarkerlist[i];
334 
335  // cout << i<< " P->vertexlist[j]: " << In.facetmarkerlist[i] << endl;
336  } // for(i=0;i<N_Faces; i++)
337 
338  BDComp_Min--;
339 
340  for(i=0;i<N_Faces; i++)
341  In.facetmarkerlist[i] -= BDComp_Min;
342 
343 
345  while (!dat.eof())
346  {
347  dat >> line;
348 
349  if ( (!strcmp(line, "Tetrahedron")) || (!strcmp(line, "tetrahedron")) || (!strcmp(line, "TETRAHEDRON")) )
350  {
351  N_CVert = 4;
352  dat.getline (line, 99);
353  dat >> N_Cells;
354  break;
355  }
356  // read until end of line
357  dat.getline (line, 99);
358  }
359 
360  In.numberoftetrahedra = N_Cells;
361  In.numberofcorners = N_CVert;
362  In.numberoftetrahedronattributes = 1;
363  In.tetrahedronlist = new int[N_Cells * 4];
364  In.tetrahedronattributelist = new REAL[N_Cells];
365  In.tetrahedronvolumelist = new REAL[N_Cells];
366 
367  for(i=0; i<N_Cells; i++)
368  {
369  dat.getline (line, 99);
370  plist = &(In.tetrahedronlist[i * 4]);
371 
372  for(j=0;j<N_CVert;j++)
373  {
374  dat >> k;
375  plist[j] = k -1;
376  }// for(j=0;j<N_CV
377 
378 
379  dat >> In.tetrahedronattributelist[i];
380 // In.tetrahedronvolumelist [i] = 100.;
381  } // for(i=0; i<N_
382 
383 
384 
385 // cout << i<< " N_Cells : " << N_Cells << endl;
386 // exit(0);
387 
388  dat.close();
389 // exit(0);
390 } // ReadMeditMeshWithCells
391 
392 
393 void ReadMeditMesh(char *SMESH, tetgenio &In)
394 {
395  int i, j, k, dimension, N_FVert, N_Faces, N_Vertices;
396  int BDComp_Min=10000000;
397 
398  char line[100];
399  tetgenio::facet *F;
400  tetgenio::polygon *P;
401 
402  std::ifstream dat(SMESH);
403 
404  if (!dat)
405  {
406  cerr << "cannot open '" << SMESH << "' for input" << endl;
407  exit(-1);
408  }
409 
410  // check the dimension
411  while (!dat.eof())
412  {
413  dat >> line;
414 
415  if ( (!strcmp(line, "Dimension")) || (!strcmp(line, "dimension")) || (!strcmp(line, "DIMENSION")))
416  {
417  dat.getline (line, 99);
418  dat >> dimension;
419  break;
420  }
421  // read until end of line
422  dat.getline (line, 99);
423  }
424 
425  if(dimension!=3)
426  {
427  cerr << "dimension: " << dimension << endl;
428  exit(-1);
429  }
430 
431  // find the N_Vertices
432  while (!dat.eof())
433  {
434  dat >> line;
435 
436  if ( (!strcmp(line, "Vertices")) || (!strcmp(line, "vertices")) || (!strcmp(line, "VERTICES")) )
437  {
438  dat.getline (line, 99);
439  dat >> N_Vertices;
440  break;
441  }
442  // read until end of line
443  dat.getline (line, 99);
444  }
445 
446  In.numberofpoints = N_Vertices;
447  In.pointlist = new double[3*N_Vertices];
448 
449  //read from file
450  for(i=0;i<N_Vertices; i++)
451  {
452  dat.getline (line, 99);
453  dat >> In.pointlist[3*i] >> In.pointlist[3*i+1] >> In.pointlist[3*i+2];
454  //cout<< i << " vert X: " <<In.pointlist[3*i] << " vert Y: " <<In.pointlist[3*i+1] <<endl;
455  }
456 
457  // find the N_Triangles
458  // find the N_Vertices
459  while (!dat.eof())
460  {
461  dat >> line;
462 
463  if ( (!strcmp(line, "Triangles")) || (!strcmp(line, "triangles")) || (!strcmp(line, "TRIANGLES")) )
464  {
465  N_FVert = 3;
466  dat.getline (line, 99);
467  dat >> N_Faces;
468  break;
469  }
470  else if ( (!strcmp(line, "Quadrilaterals")) || (!strcmp(line, "quadrilaterals")) || (!strcmp(line, "QUADRILATERALS")) )
471  {
472  N_FVert = 4;
473  dat.getline (line, 99);
474  dat >> N_Faces;
475  break;
476  }
477 
478  // read until end of line
479  dat.getline (line, 99);
480  }
481 
482  In.numberoffacets = N_Faces;
483  In.facetlist = new tetgenio::facet[In.numberoffacets];
484  In.facetmarkerlist = new int[In.numberoffacets];
485 
486  for(i=0;i<N_Faces; i++)
487  {
488  dat.getline (line, 99);
489 
490  F = &In.facetlist[i];
491  tetgenio::init(F);
492  F->numberofpolygons = 1;
493  F->polygonlist = new tetgenio::polygon[F->numberofpolygons];
494  F->numberofholes = 0;
495  F->holelist = NULL;
496  P = &F->polygonlist[0];
497  tetgenio::init(P);
498  P->numberofvertices = N_FVert;
499  P->vertexlist = new int[P->numberofvertices];
500 
501  for(j=0;j<N_FVert;j++)
502  {
503  dat >> k;
504  P->vertexlist[j] = k-1; // c numbering
505  }
506 
507  dat >> In.facetmarkerlist[i];
508 
509  if(BDComp_Min > In.facetmarkerlist[i])
510  BDComp_Min=In.facetmarkerlist[i];
511 
512  // cout << i<< " P->vertexlist[j]: " << In.facetmarkerlist[i] << endl;
513  } // for(i=0;i<N_Faces; i++)
514 
515  BDComp_Min--;
516 
517  for(i=0;i<N_Faces; i++)
518  In.facetmarkerlist[i] -= BDComp_Min;
519 
520 // cout << i<< " P->vertexlist[j]: " << BDComp_Min << endl;
521 
522 
523  dat.close();
524 // exit(0);
525 } // ReadMeditMesh
526 
527 
528 void DeleteDomain(TDomain *&Domain)
529 {
530  int i, j, k, N_Cells, N_RootCells;
531  int CurrVertex, N_Faces, N_Vertices, ID;
532 
533  TBaseCell **CellTree, **SurfCellTree, *cell;
534  TGridCell **DelCell;
535  TVertex **VertexDel;
536  TCollection *coll;
537 
538 
539  Domain->GetTreeInfo(CellTree,N_RootCells);
540  coll = Domain->GetCollection(It_Finest, 0);
541  N_Cells = coll->GetN_Cells();
542 
543 
544  // cout<<"N_RootCells: "<<N_RootCells<<endl;
545  // remove all existing vertices and joints
546  VertexDel = new TVertex*[8*N_RootCells];
547  CurrVertex = 0;
548 
549  for(i=0;i<N_Cells;i++)
550  {
551  cell = coll->GetCell(i);
552  N_Faces = cell->GetN_Faces();
553  N_Vertices = cell->GetN_Vertices();
554  for(j=0;j<N_Faces;j++)
555  {
556  if(CurrVertex==0)
557  {
558  VertexDel[CurrVertex++] = cell->GetVertex(j);
559  }
560  else
561  {
562  ID = 0;
563  for(k=0;k<CurrVertex;k++)
564  if(VertexDel[k]==cell->GetVertex(j))
565  {
566  ID = 1; break;
567  }
568  if(ID!=1)
569  {
570  VertexDel[CurrVertex++] = cell->GetVertex(j);
571  }
572  } // else if(CurrVertex==0)
573 
574  ID = 0;
575  for(k=0;k<CurrVertex;k++)
576  if(VertexDel[k]==cell->GetVertex((j+1)%N_Vertices))
577  {
578  ID = 1; break;
579  }
580  if(ID!=1)
581  {
582  VertexDel[CurrVertex++] = cell->GetVertex((j+1)%N_Vertices);
583  }
584 
585  ID = 0;
586  for(k=0;k<CurrVertex;k++)
587  if(VertexDel[k]==cell->GetVertex((j+2)%N_Vertices))
588  {
589  ID = 1; break;
590  }
591  if(ID!=1)
592  {
593  VertexDel[CurrVertex++] = cell->GetVertex((j+2)%N_Vertices);
594  }
595  if(N_Faces==6) // If old cell is hexahedrol
596  {
597  ID = 0;
598  for(k=0;k<CurrVertex;k++)
599  if(VertexDel[k]==cell->GetVertex((j+4)%N_Vertices))
600  {
601  ID = 1; break;
602  }
603  if(ID!=1)
604  {
605  VertexDel[CurrVertex++] = cell->GetVertex((j+4)%N_Vertices);
606  }
607  }
608  } // for j
609  } // for i
610 
611  for(i=0;i<CurrVertex;i++)
612  delete VertexDel[i];
613  delete [] VertexDel;
614  OutPut(CurrVertex<<" vertices were deleted"<<endl);
615 
616  // remove all existing cells and joints
617 
618  for(i=0;i<N_RootCells;i++)
619  delete (TGridCell*)CellTree[i];
620  delete [] CellTree;
621  OutPut(N_RootCells<<" cells were deleted"<<endl);
622 
623 
624 /*
625  cout << "DeleteDomain" <<endl;
626  exit(0);
627  */
628 }
629 
630 
631 void TetrameshCreate(TDomain *&Domain)
632 {
633  int i, j, k, l, dimension, N_Vertices;
634  int N_FVert, N_Faces, *Facemarkerlist, *Facelist, N_RootCells;
635  int v1, v2, v3, v4, CellMarker, RefLevel=0, *Tetrahedrals, *PointNeighb, maxEpV, *Tetrahedrals_loc;
636  int a, b, c, Neib[2], Neighb_tmp, CurrNeib, len1, len2, len3, CurrComp, N_Points ;
637  int N_Cells, MaxLen, jj, N_SourcePts;
638  const int *EdgeVertex, *TmpFV, *TmpLen;
639  double X, Y, Z, DispX, DispY, DispZ;
640  double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
641  double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
642  double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
643  char *MESH, line[100];
644  MESH = TDatabase::ParamDB->SMESHFILE;
645  TVertex **NewVertices;
646  TBaseCell **CellTree, *cell;
647  TJoint *Joint;
648  TBoundComp3D *bdcomp;
649  TShapeDesc *ShapeDesc;
650  TCollection *coll;
651 
652 #ifdef _MPI
653  int rank, size;
654  MPI_Comm Comm = TDatabase::ParamDB->Comm;
655  MPI_Comm_rank(Comm, &rank);
656  MPI_Comm_size(Comm, &size);
657 #endif
658 
660  DeleteDomain(Domain);
662  std::ifstream dat(MESH);
663 
664  if (!dat)
665  {
666  cerr << "cannot open '" << MESH << "' for input" << endl;
667  exit(-1);
668  }
669 
670  // check the dimension
671  while (!dat.eof())
672  {
673  dat >> line;
674 
675  if ( (!strcmp(line, "Dimension")) || (!strcmp(line, "dimension")) || (!strcmp(line, "DIMENSION")))
676  {
677  dat.getline (line, 99);
678  dat >> dimension;
679  break;
680  }
681  // read until end of line
682  dat.getline (line, 99);
683  }
684 
685  if(dimension!=3)
686  {
687  cerr << "dimension: " << dimension << endl;
688  exit(-1);
689  }
690  // find the N_Vertices
691  while (!dat.eof())
692  {
693  dat >> line;
694 
695  if ( (!strcmp(line, "Vertices")) || (!strcmp(line, "vertices")) || (!strcmp(line, "VERTICES")) )
696  {
697  dat.getline (line, 99);
698  dat >> N_Vertices;
699  break;
700  }
701  // read until end of line
702  dat.getline (line, 99);
703  }
704 
705  // generate new vertices
706  NewVertices = new TVertex*[N_Vertices];
707 
708  //read from file
709  for(i=0;i<N_Vertices; i++)
710  {
711  dat.getline (line, 99);
712  dat >> X >> Y >> Z;
713  NewVertices[i] = new TVertex( (X - DispX), (Y - DispY), (Z- DispZ));
714  // set bounding box
715  if (X > Xmax) Xmax = X;
716  if (X < Xmin) Xmin = X;
717  if (Y > Ymax) Ymax = Y;
718  if (Y < Ymin) Ymin = Y;
719  if (Z > Zmax) Zmax = Z;
720  if (Z < Zmin) Zmin = Z;
721  }
722  // set bounding box
723  StartX = Xmin;
724  StartY = Ymin;
725  StartZ = Zmin;
726  BoundX = Xmax - Xmin;
727  BoundY = Ymax - Ymin;
728  BoundZ = Zmax - Zmin;
729 
730  Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
731  // find the N_RootCells
732  while (!dat.eof())
733  {
734  dat >> line;
735 
736  if ( (!strcmp(line, "Tetrahedra")) || (!strcmp(line, "tetrahedra")) || (!strcmp(line, "TETRAHEDRA")) )
737  {
738  dat.getline (line, 99);
739  dat >> N_RootCells;
740  break;
741  }
742  // read until end of line
743  dat.getline (line, 99);
744  }
745  // generate new cells
746  cout << "number of root cells "<<N_RootCells <<endl;
747  if (CellTree == NULL)
748  cout << "no celltree"<<endl;
749  CellTree = new TBaseCell*[N_RootCells];
750  Tetrahedrals = new int[4*N_RootCells];
751 
752  for (i=0;i<N_RootCells;i++)
753  {
754  dat.getline (line, 99);
755  dat >> v1 >> v2 >> v3 >> v4 >> CellMarker;
756  Tetrahedrals[4*i ] = v1 -1;
757  Tetrahedrals[4*i + 1] = v2 -1;
758  Tetrahedrals[4*i + 2] = v3 -1;
759  Tetrahedrals[4*i + 3] = v4 -1;
760 
761 // for(j=0;j<N_Domains;j++)
762 // {
763 // if (CellMarker==j+1)
764 // {
765 // ncellsindomain[j]+=1;
766 // continue;
767 // }
768 // }
769 
770  CellTree[i] = new TMacroCell(TDatabase::RefDescDB[Tetrahedron], RefLevel);
771  CellTree[i]->SetRegionID(CellMarker);
772  CellTree[i]->SetAsLayerCell(1);
773  CellTree[i]->SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
774  CellTree[i]->SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
775  CellTree[i]->SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
776  CellTree[i]->SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
777 
778  CellTree[i]->SetClipBoard(i);
779  ((TMacroCell *) CellTree[i])->SetSubGridID(0);
780  }
781 
782  dat.close();
783  Domain->SetTreeInfo(CellTree, N_RootCells);
784  // initialize iterators
785  TDatabase::IteratorDB[It_EQ]->SetParam(Domain);
786  TDatabase::IteratorDB[It_LE]->SetParam(Domain);
787  TDatabase::IteratorDB[It_Finest]->SetParam(Domain);
788  TDatabase::IteratorDB[It_Between]->SetParam(Domain);
789  TDatabase::IteratorDB[It_OCAF]->SetParam(Domain);
790 
791  // search neighbours
792  PointNeighb = new int[N_Vertices];
793 #ifdef _MPI
794  if(rank==0)
795 #endif
796  cout<<"Numberofpoints "<<N_Vertices<<endl;
797  memset(PointNeighb, 0, N_Vertices*SizeOfInt);
798 
799  for (i=0;i<4*N_RootCells;i++)
800  PointNeighb[Tetrahedrals[i]]++;
801 
802  maxEpV = 0;
803  for (i=0;i<N_Vertices;i++)
804  if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
805  delete [] PointNeighb;
806 
807 #ifdef _MPI
808  if(rank==0)
809 #endif
810  cout<<"maxEpV "<< maxEpV<<endl;
811  PointNeighb = new int[++maxEpV * N_Vertices];
812  memset(PointNeighb, 0, maxEpV*N_Vertices*SizeOfInt);
813  // every vertex contains "maxEpV" columns
814  // for every vertex at first colomn contains the number of cells containing this vertex
815  // at further columns we set the index of corresponding cells containing this vertex
816  // cout<<"maxEpV*N_Vertices "<<maxEpV*N_Vertices<<endl;
817  for(i=0;i<4*N_RootCells;i++)
818  {
819  j = Tetrahedrals[i]*maxEpV;
820  PointNeighb[j]++;
821  PointNeighb[j + PointNeighb[j]] = i / 4;
822  }
824  coll=Domain->GetCollection(It_Finest, 0);
825  N_Cells = coll->GetN_Cells();
826 
827  for(i=0;i<N_Cells;i++)
828  {
829  cell = coll->GetCell(i);
830  ShapeDesc= cell->GetShapeDesc();
831  ShapeDesc->GetFaceVertex(TmpFV, TmpLen, MaxLen);
832  N_Faces = cell->GetN_Faces();
833  Tetrahedrals_loc = Tetrahedrals+4*i;
834 
835  for(jj=0;jj<N_Faces;jj++)
836  if(cell->GetJoint(jj) == NULL)
837  {
838  N_Points = TmpLen[jj];
839  if(N_Points!=3)
840  {
841  cerr << "Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
842  exit(-1);
843  }
844 
845  //printf(" TetrameshGen Null joint \n" );
846  a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
847  b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
848  c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];
849  Neib[0] = -1;
850  Neib[1] = -1;
851  CurrNeib = 0;
852  len1 = PointNeighb[a*maxEpV];
853  len2 = PointNeighb[b*maxEpV];
854  len3 = PointNeighb[c*maxEpV];
855  // find the index of the cells containing current face with point indices a,b,c
856  for (j=1;j<=len1;j++)
857  {
858  Neighb_tmp = PointNeighb[a*maxEpV + j];
859  for (k=1;k<=len2;k++)
860  {
861  if(Neighb_tmp == PointNeighb[b*maxEpV + k])
862  {
863  for(l=1;l<=len3;l++)
864  if(Neighb_tmp == PointNeighb[c*maxEpV + l])
865  {
866  Neib[CurrNeib++] = Neighb_tmp;
867  break;
868  }
869  }
870  }
871  if (CurrNeib == 2) break;
872  } // for (j=1;j<=len1;j++)
873  if(CurrNeib==1)
874  {
875  bdcomp = Domain->GetBdPart(0)->GetBdComp(CurrComp);
876 
877  if(bdcomp->GetTSofXYZ(NewVertices[a]->GetX(), NewVertices[a]->GetY(),
878  NewVertices[a]->GetY(), T[1], S[1]) ||
879  bdcomp->GetTSofXYZ(NewVertices[b]->GetX(), NewVertices[b]->GetY(),
880  NewVertices[b]->GetY(), T[2], S[2]) ||
881  bdcomp->GetTSofXYZ(NewVertices[c]->GetX(), NewVertices[c]->GetY(),
882  NewVertices[c]->GetY(), T[3], S[3]) )
883  {
884  cerr<<"Error: could not set parameter values"<<endl;
885  OutPut(NewVertices[a]<<endl);
886  OutPut(NewVertices[b]<<endl);
887  OutPut(NewVertices[c]<<endl);
888  exit(0);
889  }
890  Joint = new TBoundFace(bdcomp, T, S);
891  }
892  else
893  {
894  if (CurrNeib != 2)
895  cerr << "Error !!!!!!!! not enough neighbours!" << endl;
896  Joint = new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
897  }
898  // First element containing the current face
899  // find the local index for the point 'a' on the cell
900  for (j=0;j<4;j++)
901  if (Tetrahedrals[4*Neib[0]+j] == a) break;
902  // find the local index for the point 'b' on the cell
903  for (k=0;k<4;k++)
904  if (Tetrahedrals[4*Neib[0]+k] == b) break;
905  // find the local index for the point 'c' on the cell
906  for (l=0;l<4;l++)
907  if (Tetrahedrals[4*Neib[0]+l] == c) break;
908  l = l*100 + k*10 + j;
909 
910  switch (l) // j will contain the local index for the current face
911  {
912  case 210: case 21: case 102:
913  case 120: case 12: case 201:
914  j = 0;
915  break;
916  case 310: case 31: case 103:
917  case 130: case 13: case 301:
918  j = 1;
919  break;
920  case 321: case 132: case 213:
921  case 231: case 123: case 312:
922  j = 2;
923  break;
924  case 230: case 23: case 302:
925  case 320: case 32: case 203:
926  j = 3;
927  break;
928 
929  default:
930  Error("Unable to set the face !!!!!!!!!!!!" << endl);
931  exit(0);
932  } // switch (l)
933  CellTree[Neib[0]]->SetJoint(j, Joint);
934 
935  // second cell containing the current face
936  if (Neib[1] != -1) // second element containing the current face
937  {
938  // find the local index for the point 'a' on the cell
939  for (j=0;j<4;j++)
940  if (Tetrahedrals[4*Neib[1]+j] == a) break;
941 
942  // find the local index for the point 'b' on the cell
943  for (k=0;k<4;k++)
944  if (Tetrahedrals[4*Neib[1]+k] == b) break;
945 
946  // find the local index for the point 'c' on the cell
947  for (l=0;l<4;l++)
948  if (Tetrahedrals[4*Neib[1]+l] == c) break;
949 
950  l = l*100 + k*10 + j;
951  switch (l) // j will contain the local index for the current face
952  {
953  case 210: case 21: case 102:
954  case 120: case 12: case 201:
955  j = 0;
956  break;
957  case 310: case 31: case 103:
958  case 130: case 13: case 301:
959  j = 1;
960  break;
961  case 321: case 132: case 213:
962  case 231: case 123: case 312:
963  j = 2;
964  break;
965  case 230: case 23: case 302:
966  case 320: case 32: case 203:
967  j = 3;
968  break;
969 
970  default:
971  Error("Unable to set the face !!!!!!!!!!!!" << endl);
972  exit(0);
973  }
974  CellTree[Neib[1]]->SetJoint(j, Joint);
975  } // if (Neib[1] != -1)
976 
977  if (Joint->GetType() == InterfaceJoint3D ||
978  Joint->GetType() == IsoInterfaceJoint3D)
979  {
980  ((TInterfaceJoint3D*)Joint)->SetMapType();
981  ((TInterfaceJoint3D*)(Joint))->CheckOrientation();
982  }
983  else
984  if (Joint->GetType() == JointEqN)
985  Joint->SetMapType();
986 
987  } // if(cell->GetJoint(j)
988  } // for(i=0;i<N_Cells;i++)
989 delete [] PointNeighb;
990  cout<< "Yestetrameshcreate " <<"\n";
991 }
992 
993 void TetrameshGen(TDomain *&Domain)
994 {
995  //======================================================================
996  // Tetgen for grid generation begin
997  //======================================================================
998  int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
999  int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
1000  int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
1001  int *Tetrahedrals, *PointNeighb, dimension;
1002  int *Facelist, *Facemarkerlist;
1003 
1004  double *Coordinates, N_x, N_y, N_z;
1005  double *Vertices;
1006  double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
1007  double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
1008  double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
1009 
1010  tetgenio In, Out;
1011 
1012  TBaseCell **CellTree, **SurfCellTree;
1013  TGridCell **DelCell;
1014  TVertex **VertexDel, **NewVertices, **NewSurfVertices;
1015  TBoundPart *BoundPart;
1016  TBdPlane **UpdateFaceParams;
1017  TJoint *Joint;
1018  TBoundComp3D *bdcomp;
1019  TCollection *coll, *SurfColl;
1020  TBaseCell *cell;
1021  TBdSphere *UpdateParam;
1022 
1023  char *SMESH, line[100];
1024  SMESH = TDatabase::ParamDB->SMESHFILE;
1025 
1026 #ifdef _MPI
1027  int rank, size;
1028  MPI_Comm Comm = TDatabase::ParamDB->Comm;
1029 
1030  MPI_Comm_rank(Comm, &rank);
1031  MPI_Comm_size(Comm, &size);
1032 
1033 #endif
1034 
1035 
1036 
1038 #ifdef _MPI
1039  if(rank==0)
1040 #endif
1041  {
1042 
1043 // In.initialize();
1044 // Out.initialize();
1045 
1046  std::ostringstream opts;
1047  opts << " ";
1048 
1049  opts.seekp(std::ios::beg);
1050  opts<<'p'; // Tetrahedralize the PLC. Switches are chosen to read a PLC (p)
1051 // opts<<'r'; // -r Reconstructs a previously generated mesh
1052 // opts<<"q"<<1.49; // quality mesh generation(q) with a specified quality bound
1053 // opts<<"a"<<0.1; // maximum volume constraint
1054 // opts<<'i'; // Inserts a list of additional points into mesh.
1055  opts<<'z'; // numbers all output items starting from zero
1056 // opts<<'d'; // Detect intersections of PLC facets.
1057  opts<<'f'; // Outputs all faces (including non-boundary)
1058  opts<<'e'; // Outputs a list of edges of the triangulation
1059 // opts<<'I'; // Suppresses mesh iteration numbers.
1060  opts<<'C'; // Checks the consistency of the final mesh.
1061 // opts<<'Q'; // Quiet: No terminal output except errors.
1062 // opts<<'g'; // Outputs mesh to .mesh file for viewing by Medit
1063  opts<<'Y'; // Suppresses boundary facets/segments splitting
1064 // opts<<'V'; //verbose mode
1065  opts<<ends;
1066 
1067 
1068 // cout << " SMESHFILE is " << SMESH << endl;
1070  ReadMeditMesh(SMESH, In);
1071 
1072 
1073 // for(i=0;i<In.numberofpoints;i++)
1074 // OutPut(i<<" (x, y, z) = "<<
1075 // In.pointlist[3*i]<<' '<<In.pointlist[3*i+1]<<' '<<In.pointlist[3*i+2]<<endl);
1076 
1077  // Calling tetrahedralize function of 3dtetgen mesh generator
1078  tetrahedralize((char*)opts.str().c_str(), &In, &Out);
1079 
1080  } // if(rank==0)
1081 
1082 
1083  // output: coordinates of all vertices
1084 // for(i=0;i<Out.numberofpoints;i++)
1085 // OutPut(" (x, y, z) = "<< Out.pointlist[3*i]<<' '<<Out.pointlist[3*i+1]<<' '<<Out.pointlist[3*i+2]<<endl);
1086 
1087  Domain->GetTreeInfo(CellTree,N_RootCells);
1088  coll = Domain->GetCollection(It_Finest, 0);
1089  N_Cells = coll->GetN_Cells();
1090 
1091 // cout<<"N_RootCells: "<<N_RootCells<<endl;
1092  // remove all existing vertices and joints
1093  VertexDel = new TVertex*[8*N_RootCells];
1094  CurrVertex = 0;
1095 
1096  for(i=0;i<N_Cells;i++)
1097  {
1098  cell = coll->GetCell(i);
1099  N_Faces = cell->GetN_Faces();
1100  N_Vertices = cell->GetN_Vertices();
1101  for(j=0;j<N_Faces;j++)
1102  {
1103  if(CurrVertex==0)
1104  {
1105  VertexDel[CurrVertex++] = cell->GetVertex(j);
1106  }
1107  else
1108  {
1109  ID = 0;
1110  for(k=0;k<CurrVertex;k++)
1111  if(VertexDel[k]==cell->GetVertex(j))
1112  {
1113  ID = 1; break;
1114  }
1115  if(ID!=1)
1116  {
1117  VertexDel[CurrVertex++] = cell->GetVertex(j);
1118  }
1119  } // else if(CurrVertex==0)
1120 
1121  ID = 0;
1122  for(k=0;k<CurrVertex;k++)
1123  if(VertexDel[k]==cell->GetVertex((j+1)%N_Vertices))
1124  {
1125  ID = 1; break;
1126  }
1127  if(ID!=1)
1128  {
1129  VertexDel[CurrVertex++] = cell->GetVertex((j+1)%N_Vertices);
1130  }
1131 
1132  ID = 0;
1133  for(k=0;k<CurrVertex;k++)
1134  if(VertexDel[k]==cell->GetVertex((j+2)%N_Vertices))
1135  {
1136  ID = 1; break;
1137  }
1138  if(ID!=1)
1139  {
1140  VertexDel[CurrVertex++] = cell->GetVertex((j+2)%N_Vertices);
1141  }
1142  if(N_Faces==6) // If old cell is hexahedrol
1143  {
1144  ID = 0;
1145  for(k=0;k<CurrVertex;k++)
1146  if(VertexDel[k]==cell->GetVertex((j+4)%N_Vertices))
1147  {
1148  ID = 1; break;
1149  }
1150  if(ID!=1)
1151  {
1152  VertexDel[CurrVertex++] = cell->GetVertex((j+4)%N_Vertices);
1153  }
1154  }
1155  } // for j
1156  } // for i
1157 
1158  for(i=0;i<CurrVertex;i++)
1159  delete VertexDel[i];
1160  delete [] VertexDel;
1161  OutPut(CurrVertex<<" vertices were deleted"<<endl);
1162 
1163  // remove all existing cells and joints
1164 
1165  for(i=0;i<N_RootCells;i++)
1166  delete (TGridCell*)CellTree[i];
1167  delete [] CellTree;
1168  OutPut(N_RootCells<<" cells were deleted"<<endl);
1169 
1170 
1171 #ifdef _MPI
1172  if(rank==0)
1173 #endif
1174  {
1175  N_RootCells = Out.numberoftetrahedra;
1176  N_G = Out.numberofpoints;
1177 
1178  // allocate auxillary fields
1179  Coordinates = Out.pointlist;
1180  Tetrahedrals = Out.tetrahedronlist;
1181  }
1182 
1183 
1184 
1185 #ifdef _MPI
1186  MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
1187  MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
1188 
1189 
1190  if(rank!=0)
1191  {
1192  Coordinates = new double [3*N_G];
1193  Tetrahedrals = new int[4*N_RootCells];
1194  }
1195 
1196  MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);
1197  MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);
1198 #endif
1199 
1200  // generate new vertices
1201  NewVertices = new TVertex*[N_G];
1202 
1203  for (i=0;i<N_G;i++)
1204  {
1205  NewVertices[i] = new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);
1206 
1207  // set bounding box
1208  if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
1209  if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
1210  if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
1211  if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
1212  if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
1213  if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
1214  }
1215 
1216  // set bounding box
1217  StartX = Xmin;
1218  StartY = Ymin;
1219  StartZ = Zmin;
1220  BoundX = Xmax - Xmin;
1221  BoundY = Ymax - Ymin;
1222  BoundZ = Zmax - Zmin;
1223 
1224 
1225  Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
1226 // cout<<Xmin <<" "<<Ymin <<" "<<Zmin<<endl;
1227 // cout<<Xmax <<" "<<Ymax <<" "<<Zmax<<endl;
1228 
1229  CellTree = new TBaseCell*[N_RootCells];
1230  // output of each tetraheron vertex indices (four vertices for each)
1231 // for (i=0;i<N_RootCells;i++)
1232 // cout<< Tetrahedrals[4*i]<<" "<<Tetrahedrals[4*i + 1]<<" "
1233 // <<Tetrahedrals[4*i + 2]<<" "<<Tetrahedrals[4*i + 3]<<endl;
1234 
1235 
1236  for (i=0;i<N_RootCells;i++)
1237  {
1238  CellTree[i] = new TMacroCell(TDatabase::RefDescDB[Tetrahedron],
1239  RefLevel);
1240 
1241  CellTree[i]->SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
1242  CellTree[i]->SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
1243  CellTree[i]->SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
1244  CellTree[i]->SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
1245 
1246  CellTree[i]->SetClipBoard(i);
1247  ((TMacroCell *) CellTree[i])->SetSubGridID(0);
1248  }
1249 
1250  Domain->SetTreeInfo(CellTree, N_RootCells);
1251 
1252 
1253  // initialize iterators
1254  TDatabase::IteratorDB[It_EQ]->SetParam(Domain);
1255  TDatabase::IteratorDB[It_LE]->SetParam(Domain);
1256  TDatabase::IteratorDB[It_Finest]->SetParam(Domain);
1257  TDatabase::IteratorDB[It_Between]->SetParam(Domain);
1258  TDatabase::IteratorDB[It_OCAF]->SetParam(Domain);
1259 
1260 
1261  // search neighbours
1262  PointNeighb = new int[N_G];
1263 #ifdef _MPI
1264  if(rank==0)
1265 #endif
1266  cout<<"numberofpoints "<<N_G<<endl;
1267  memset(PointNeighb, 0, N_G *SizeOfInt);
1268 
1269  for (i=0;i<4*N_RootCells;i++)
1270  PointNeighb[Tetrahedrals[i]]++;
1271 
1272  for (i=0;i<N_G;i++)
1273  if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
1274  delete [] PointNeighb;
1275 
1276 #ifdef _MPI
1277  if(rank==0)
1278 #endif
1279  cout<<"maxEpV "<< maxEpV<<endl;
1280 
1281  PointNeighb = new int[++maxEpV * N_G];
1282 
1283  memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);
1284 
1285  // every vertex contains "maxEpV" columns
1286  // for every vertex at first colomn contains the number of cells containing this vertex
1287  // at further columns we set the index of corresponding cells containing this vertex
1288  // cout<<"maxEpV*N_G "<<maxEpV*N_G<<endl;
1289 
1290  for(i=0;i<4*N_RootCells;i++)
1291  {
1292  j = Tetrahedrals[i]*maxEpV;
1293  PointNeighb[j]++;
1294  //cout<<"j + PointNeighb[j] " << j <<endl;
1295  PointNeighb[j + PointNeighb[j]] = i / 4;
1296  }
1297  // output of PointNeighb columns for each point
1298 // for (i=0;i<N_G;i++)
1299 // {
1300 // for (j=0;j<maxEpV;j++)
1301 // cout<<" "<< PointNeighb[i*maxEpV+j];
1302 // cout<<endl;
1303 // }
1304 
1305 
1306 
1307  // generate new faces
1308 
1309 #ifdef _MPI
1310  if(rank==0)
1311 #endif
1312  {
1313  N_G = Out.numberoftrifaces;
1314  Facelist = Out.trifacelist;
1315  Facemarkerlist = Out.trifacemarkerlist;
1316  }
1317 
1318 #ifdef _MPI
1319  MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
1320 
1321  if(rank!=0)
1322  {
1323  Facelist = new int [3*N_G];
1324  Facemarkerlist = new int[N_G];
1325  }
1326 
1327 
1328  MPI_Bcast(Facelist, 3*N_G, MPI_INT, 0, Comm);
1329  MPI_Bcast(Facemarkerlist, N_G, MPI_INT, 0, Comm);
1330 
1331  if(rank==0)
1332 #endif
1333  cout<<"numberoftrifaces "<<N_G<<endl;
1334 
1335  for (i=0;i<N_G;i++)
1336  {
1337  a = Facelist[3*i];
1338  b = Facelist[3*i+1];
1339  c = Facelist[3*i+2];
1340 
1341 // cout<<" "<< a<<" "<< b<<" "<< c<<endl;
1342 
1343  Neib[0] = -1;
1344  Neib[1] = -1;
1345  CurrNeib = 0;
1346 
1347  len1 = PointNeighb[a*maxEpV];
1348  len2 = PointNeighb[b*maxEpV];
1349  len3 = PointNeighb[c*maxEpV];
1350 
1351  // find the index of the cells containing current face with point indices a,b,c
1352  for (j=1;j<=len1;j++)
1353  {
1354  Neighb_tmp = PointNeighb[a*maxEpV + j];
1355  for (k=1;k<=len2;k++)
1356  {
1357  if (Neighb_tmp == PointNeighb[b*maxEpV + k])
1358  {
1359  for (l=1;l<=len3;l++)
1360  if (Neighb_tmp == PointNeighb[c*maxEpV + l])
1361  {
1362  Neib[CurrNeib++] = Neighb_tmp;
1363  break;
1364  }
1365  }
1366  }
1367  if (CurrNeib == 2) break;
1368  }
1369  // cout<<"CurrNeib " << CurrNeib <<endl;
1370 // cout<<"Facemarkerlist[i] : "<<Facemarkerlist[i]<<endl;
1371  if (Facemarkerlist[i]) // 0 for inner edges and Boundcomp+1 for Boundedge respect
1372  {
1373 
1374  CurrComp = Facemarkerlist[i] - 1;
1375 // cout<<"Boundary face CurrComp: "<<CurrComp<<endl;
1376 // exit(0);
1377  CurrComp = 0; // not yet implemented fully
1378 
1379 
1380 
1381  bdcomp = Domain->GetBdPart(0)->GetBdComp(CurrComp);
1382 
1383 
1384  if(bdcomp->GetTSofXYZ(NewVertices[a]->GetX(), NewVertices[a]->GetY(),
1385  NewVertices[a]->GetY(), T[1], S[1]) ||
1386  bdcomp->GetTSofXYZ(NewVertices[b]->GetX(), NewVertices[b]->GetY(),
1387  NewVertices[b]->GetY(), T[2], S[2]) ||
1388  bdcomp->GetTSofXYZ(NewVertices[c]->GetX(), NewVertices[c]->GetY(),
1389  NewVertices[c]->GetY(), T[3], S[3]) )
1390  {
1391  cerr<<"Error: could not set parameter values"<<endl;
1392  OutPut(NewVertices[a]<<endl);
1393  OutPut(NewVertices[b]<<endl);
1394  OutPut(NewVertices[c]<<endl);
1395  exit(0);
1396  }
1397 
1398  if (CurrNeib == 2)
1399  {
1400  if(bdcomp->IsFreeBoundary())
1401  Joint = new TIsoInterfaceJoint3D(bdcomp, T, S,
1402  CellTree[Neib[0]], CellTree[Neib[1]] );
1403  else
1404  Joint = new TInterfaceJoint3D(bdcomp, T, S,
1405  CellTree[Neib[0]], CellTree[Neib[1]] );
1406  }
1407  else
1408  {
1409  if(bdcomp->IsFreeBoundary())
1410  Joint = new TIsoBoundFace(bdcomp, T, S);
1411  else
1412  Joint = new TBoundFace(bdcomp, T, S);
1413  }
1414  }
1415  else
1416  {
1417 // // cout<<"Inner face"<<endl;
1418  if (CurrNeib != 2)
1419  cerr << "Error !!!!!!!! not enough neighbours!" << endl;
1420 
1421  Joint = new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
1422 
1423  }
1424 
1425  // First element containing the current face
1426  // find the local index for the point 'a' on the cell
1427  for (j=0;j<4;j++)
1428  if (Tetrahedrals[4*Neib[0]+j] == a) break;
1429 
1430  // find the local index for the point 'b' on the cell
1431  for (k=0;k<4;k++)
1432  if (Tetrahedrals[4*Neib[0]+k] == b) break;
1433 
1434  // find the local index for the point 'c' on the cell
1435  for (l=0;l<4;l++)
1436  if (Tetrahedrals[4*Neib[0]+l] == c) break;
1437 
1438  l = l*100 + k*10 + j;
1439 
1440 // cout<<""<< l <<endl;
1441 
1442  switch (l) // j will contain the local index for the current face
1443  {
1444  case 210: case 21: case 102:
1445  case 120: case 12: case 201:
1446  j = 0;
1447  break;
1448  case 310: case 31: case 103:
1449  case 130: case 13: case 301:
1450  j = 1;
1451  break;
1452  case 321: case 132: case 213:
1453  case 231: case 123: case 312:
1454  j = 2;
1455  break;
1456  case 230: case 23: case 302:
1457  case 320: case 32: case 203:
1458  j = 3;
1459  break;
1460 
1461  default:
1462  Error("Unable to set the face !!!!!!!!!!!!" << endl);
1463  exit(0);
1464  }
1465  CellTree[Neib[0]]->SetJoint(j, Joint);
1466 
1467  if (Neib[1] != -1) // second element containing the current face
1468  {
1469  // find the local index for the point 'a' on the cell
1470  for (j=0;j<4;j++)
1471  if (Tetrahedrals[4*Neib[1]+j] == a) break;
1472 
1473  // find the local index for the point 'b' on the cell
1474  for (k=0;k<4;k++)
1475  if (Tetrahedrals[4*Neib[1]+k] == b) break;
1476 
1477  // find the local index for the point 'c' on the cell
1478  for (l=0;l<4;l++)
1479  if (Tetrahedrals[4*Neib[1]+l] == c) break;
1480 
1481  l = l*100 + k*10 + j;
1482 
1483 // cout<<""<< l <<endl;
1484 
1485  switch (l) // j will contain the local index for the current face
1486  {
1487  case 210: case 21: case 102:
1488  case 120: case 12: case 201:
1489  j = 0;
1490  break;
1491  case 310: case 31: case 103:
1492  case 130: case 13: case 301:
1493  j = 1;
1494  break;
1495  case 321: case 132: case 213:
1496  case 231: case 123: case 312:
1497  j = 2;
1498  break;
1499  case 230: case 23: case 302:
1500  case 320: case 32: case 203:
1501  j = 3;
1502  break;
1503 
1504  default:
1505  Error("Unable to set the face !!!!!!!!!!!!" << endl);
1506  exit(0);
1507  }
1508  CellTree[Neib[1]]->SetJoint(j, Joint);
1509  }
1510 
1511  if (Joint->GetType() == InterfaceJoint3D ||
1512  Joint->GetType() == IsoInterfaceJoint3D)
1513  {
1514  ((TInterfaceJoint3D*)Joint)->SetMapType();
1515  ((TInterfaceJoint3D*)(Joint))->CheckOrientation();
1516  }
1517  else
1518  if (Joint->GetType() == JointEqN)
1519  Joint->SetMapType();
1520 
1521  }
1522 
1523 #ifdef _MPI
1524  if(rank==0)
1525 #endif
1526  {
1527 // In.deinitialize();
1528 // Out.deinitialize();
1529  }
1530 
1531 
1532 
1533 
1534 #ifdef _MPI
1535  if(rank!=0)
1536  {
1537  delete [] Tetrahedrals ;
1538  delete [] Coordinates;
1539  delete [] Facelist;
1540  delete [] Facemarkerlist;
1541  }
1542 #endif
1543 
1544  delete [] PointNeighb;
1545 
1546 // #ifdef _MPI
1547 // if(rank==3)
1548 // printf(" TetrameshGen %d \n", N_G );
1549 // MPI_Finalize();
1550 // #endif
1551 // exit(0);
1552 //
1553 
1554 }
1555 
1556 
1558 void TetraGen(TDomain *&Domain)
1559 {
1560  //======================================================================
1561  // Tetgen for grid generation begin
1562  //======================================================================
1563  int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
1564  int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
1565  int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
1566  int *Tetrahedrals, *PointNeighb, dimension;
1567  int jj, N_Points, MaxLen, *Tetrahedrals_loc;
1568  const int *EdgeVertex, *TmpFV, *TmpLen;
1569 
1570  double *Coordinates, N_x, N_y, N_z;
1571  double *Vertices;
1572  double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
1573  double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
1574  double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
1575 
1576  tetgenio In, Out;
1577 
1578  TBaseCell **CellTree, **SurfCellTree;
1579  TGridCell **DelCell;
1580  TVertex **VertexDel, **NewVertices, **NewSurfVertices;
1581  TBoundPart *BoundPart;
1582  TBdPlane **UpdateFaceParams;
1583  TJoint *Joint;
1584  TBoundComp3D *bdcomp;
1585  TCollection *coll, *SurfColl;
1586  TBaseCell *cell;
1587  TBdSphere *UpdateParam;
1588  TShapeDesc *ShapeDesc;
1589 
1590  char *SMESH, line[100];
1591  SMESH = TDatabase::ParamDB->SMESHFILE;
1592 
1593 #ifdef _MPI
1594  int rank, size;
1595  MPI_Comm Comm = TDatabase::ParamDB->Comm;
1596 
1597  MPI_Comm_rank(Comm, &rank);
1598  MPI_Comm_size(Comm, &size);
1599 
1600 #endif
1601 
1603 #ifdef _MPI
1604  if(rank==0)
1605 #endif
1606  {
1607 
1608 // In.initialize();
1609 // Out.initialize();
1610 
1611  std::ostringstream opts;
1612  opts << " ";
1613 
1614  opts.seekp(std::ios::beg);
1615  opts<<'p'; // Tetrahedralize the PLC. Switches are chosen to read a PLC (p)
1616 // opts<<'r'; // -r Reconstructs a previously generated mesh
1617  opts<<"q"<<1.20; // quality mesh generation(q) with a specified quality bound
1618 // opts<<"a"<<0.1; // maximum volume constraint
1619 // opts<<'i'; // Inserts a list of additional points into mesh.
1620  opts<<'z'; // numbers all output items starting from zero
1621 // opts<<'d'; // Detect intersections of PLC facets.
1622  opts<<'f'; // Outputs all faces (including non-boundary)
1623  opts<<'e'; // Outputs a list of edges of the triangulation
1624 // opts<<'I'; // Suppresses mesh iteration numbers.
1625  opts<<'C'; // Checks the consistency of the final mesh.
1626 // opts<<'Q'; // Quiet: No terminal output except errors.
1627 // opts<<'g'; // Outputs mesh to .mesh file for viewing by Medit
1628  opts<<'Y'; // Suppresses boundary facets/segments splitting
1629 // opts<<'V'; //verbose mode
1630  opts<<ends;
1631 
1632 
1633 // cout << " SMESHFILE is " << SMESH << endl;
1635  ReadMeditMesh(SMESH, In);
1636 
1637 
1638 // for(i=0;i<In.numberofpoints;i++)
1639 // OutPut(i<<" (x, y, z) = "<<
1640 // In.pointlist[3*i]<<' '<<In.pointlist[3*i+1]<<' '<<In.pointlist[3*i+2]<<endl);
1641 
1642  // Calling tetrahedralize function of 3dtetgen mesh generator
1643  tetrahedralize((char*)opts.str().c_str(), &In, &Out);
1644 
1645  // output: coordinates of all vertices
1646 // for(i=0;i<Out.numberofpoints;i++)
1647 // OutPut(" (x, y, z) = "<< Out.pointlist[3*i]<<' '<<Out.pointlist[3*i+1]<<' '<<Out.pointlist[3*i+2]<<endl);
1648 
1649  } // if(rank==0)
1650 
1652  DeleteDomain(Domain);
1653 
1654 
1655 #ifdef _MPI
1656  if(rank==0)
1657 #endif
1658  {
1659  N_RootCells = Out.numberoftetrahedra;
1660  N_G = Out.numberofpoints;
1661 
1662  // allocate auxillary fields
1663  Coordinates = Out.pointlist;
1664  Tetrahedrals = Out.tetrahedronlist;
1665  }
1666 
1667 #ifdef _MPI
1668  MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
1669  MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
1670 
1671 
1672  if(rank!=0)
1673  {
1674  Coordinates = new double [3*N_G];
1675  Tetrahedrals = new int[4*N_RootCells];
1676  }
1677 
1678  MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);
1679  MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);
1680 #endif
1681 
1682  // generate new vertices
1683  NewVertices = new TVertex*[N_G];
1684 
1685  for (i=0;i<N_G;i++)
1686  {
1687  NewVertices[i] = new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);
1688 
1689  // set bounding box
1690  if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
1691  if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
1692  if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
1693  if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
1694  if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
1695  if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
1696  }
1697 
1698  // set bounding box
1699  StartX = Xmin;
1700  StartY = Ymin;
1701  StartZ = Zmin;
1702  BoundX = Xmax - Xmin;
1703  BoundY = Ymax - Ymin;
1704  BoundZ = Zmax - Zmin;
1705 
1706 
1707  Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
1708 // cout<<Xmin <<" "<<Ymin <<" "<<Zmin<<endl;
1709 // cout<<Xmax <<" "<<Ymax <<" "<<Zmax<<endl;
1710 
1711  CellTree = new TBaseCell*[N_RootCells];
1712  // output of each tetraheron vertex indices (four vertices for each)
1713 // for (i=0;i<N_RootCells;i++)
1714 // cout<< Tetrahedrals[4*i]<<" "<<Tetrahedrals[4*i + 1]<<" "
1715 // <<Tetrahedrals[4*i + 2]<<" "<<Tetrahedrals[4*i + 3]<<endl;
1716 
1717 
1718  for (i=0;i<N_RootCells;i++)
1719  {
1720  CellTree[i] = new TMacroCell(TDatabase::RefDescDB[Tetrahedron],
1721  RefLevel);
1722 
1723  CellTree[i]->SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
1724  CellTree[i]->SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
1725  CellTree[i]->SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
1726  CellTree[i]->SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
1727 
1728  CellTree[i]->SetClipBoard(i);
1729  ((TMacroCell *) CellTree[i])->SetSubGridID(0);
1730  }
1731 
1732  Domain->SetTreeInfo(CellTree, N_RootCells);
1733 
1734 
1735  // initialize iterators
1736  TDatabase::IteratorDB[It_EQ]->SetParam(Domain);
1737  TDatabase::IteratorDB[It_LE]->SetParam(Domain);
1738  TDatabase::IteratorDB[It_Finest]->SetParam(Domain);
1739  TDatabase::IteratorDB[It_Between]->SetParam(Domain);
1740  TDatabase::IteratorDB[It_OCAF]->SetParam(Domain);
1741 
1742 
1743  // search neighbours
1744  PointNeighb = new int[N_G];
1745 #ifdef _MPI
1746  if(rank==0)
1747 #endif
1748  cout<<"numberofpoints "<<N_G<<endl;
1749  memset(PointNeighb, 0, N_G *SizeOfInt);
1750 
1751  for (i=0;i<4*N_RootCells;i++)
1752  PointNeighb[Tetrahedrals[i]]++;
1753 
1754  for (i=0;i<N_G;i++)
1755  if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
1756  delete [] PointNeighb;
1757 
1758 #ifdef _MPI
1759  if(rank==0)
1760 #endif
1761  cout<<"maxEpV "<< maxEpV<<endl;
1762 
1763  PointNeighb = new int[++maxEpV * N_G];
1764 
1765  memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);
1766 
1767  // every vertex contains "maxEpV" columns
1768  // for every vertex at first colomn contains the number of cells containing this vertex
1769  // at further columns we set the index of corresponding cells containing this vertex
1770  // cout<<"maxEpV*N_G "<<maxEpV*N_G<<endl;
1771 
1772  for(i=0;i<4*N_RootCells;i++)
1773  {
1774  j = Tetrahedrals[i]*maxEpV;
1775  PointNeighb[j]++;
1776  //cout<<"j + PointNeighb[j] " << j <<endl;
1777  PointNeighb[j + PointNeighb[j]] = i / 4;
1778  }
1779 
1780 
1781  // generate new faces
1782  coll=Domain->GetCollection(It_Finest, 0);
1783  N_Cells = coll->GetN_Cells();
1784  N_G = 0;
1785  for (i=0;i<N_Cells;i++)
1786  {
1787  cell = coll->GetCell(i);
1788  ShapeDesc= cell->GetShapeDesc();
1789  ShapeDesc->GetFaceVertex(TmpFV, TmpLen, MaxLen);
1790  N_Faces = cell->GetN_Faces();
1791  Tetrahedrals_loc = Tetrahedrals+4*i;
1792 
1793  for(jj=0;jj<N_Faces;jj++)
1794  if(cell->GetJoint(jj) == NULL)
1795  {
1796  N_G++;
1797  N_Points = TmpLen[jj];
1798 
1799  if(N_Points!=3)
1800  {
1801  cerr << "Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
1802  exit(-1);
1803  }
1804 
1805 
1806  //printf(" TetrameshGen Null joint \n" );
1807  a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
1808  b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
1809  c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];
1810 
1811 // cout<<" "<< a<<" "<< b<<" "<< c<<endl;
1812 
1813  Neib[0] = -1;
1814  Neib[1] = -1;
1815  CurrNeib = 0;
1816 
1817  len1 = PointNeighb[a*maxEpV];
1818  len2 = PointNeighb[b*maxEpV];
1819  len3 = PointNeighb[c*maxEpV];
1820 
1821  // find the index of the cells containing current face with point indices a,b,c
1822  for (j=1;j<=len1;j++)
1823  {
1824  Neighb_tmp = PointNeighb[a*maxEpV + j];
1825  for (k=1;k<=len2;k++)
1826  {
1827  if (Neighb_tmp == PointNeighb[b*maxEpV + k])
1828  {
1829  for (l=1;l<=len3;l++)
1830  if (Neighb_tmp == PointNeighb[c*maxEpV + l])
1831  {
1832  Neib[CurrNeib++] = Neighb_tmp;
1833  break;
1834  }
1835  }
1836  }
1837  if (CurrNeib == 2) break;
1838  }// for (j=1;j<=len1;j++)
1839 
1840  // cout<<"CurrNeib " << CurrNeib <<endl;
1841  // cout<<"Facemarkerlist[i] : "<<Facemarkerlist[i]<<endl;
1842  if (CurrNeib == 1) // 0 for inner edges and Boundcomp+1 for Boundedge respect
1843  {
1844  CurrComp = 0; // not yet implemented fully
1845 
1846 
1847 
1848  bdcomp = Domain->GetBdPart(0)->GetBdComp(CurrComp);
1849 
1850 
1851  if(bdcomp->GetTSofXYZ(NewVertices[a]->GetX(), NewVertices[a]->GetY(),
1852  NewVertices[a]->GetY(), T[1], S[1]) ||
1853  bdcomp->GetTSofXYZ(NewVertices[b]->GetX(), NewVertices[b]->GetY(),
1854  NewVertices[b]->GetY(), T[2], S[2]) ||
1855  bdcomp->GetTSofXYZ(NewVertices[c]->GetX(), NewVertices[c]->GetY(),
1856  NewVertices[c]->GetY(), T[3], S[3]) )
1857  {
1858  cerr<<"Error: could not set parameter values"<<endl;
1859  OutPut(NewVertices[a]<<endl);
1860  OutPut(NewVertices[b]<<endl);
1861  OutPut(NewVertices[c]<<endl);
1862  exit(0);
1863  }
1864 
1865  if (CurrNeib == 2)
1866  {
1867  if(bdcomp->IsFreeBoundary())
1868  Joint = new TIsoInterfaceJoint3D(bdcomp, T, S,
1869  CellTree[Neib[0]], CellTree[Neib[1]] );
1870  else
1871  Joint = new TInterfaceJoint3D(bdcomp, T, S,
1872  CellTree[Neib[0]], CellTree[Neib[1]] );
1873  }
1874  else
1875  {
1876  if(bdcomp->IsFreeBoundary())
1877  Joint = new TIsoBoundFace(bdcomp, T, S);
1878  else
1879  Joint = new TBoundFace(bdcomp, T, S);
1880  }
1881  }
1882  else
1883  {
1884  // cout<<"Inner face"<<endl;
1885  if (CurrNeib != 2)
1886  cerr << "Error !!!!!!!! not enough neighbours!" << endl;
1887 
1888  Joint = new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
1889  }
1890 
1891  // First element containing the current face
1892  // find the local index for the point 'a' on the cell
1893  for (j=0;j<4;j++)
1894  if (Tetrahedrals[4*Neib[0]+j] == a) break;
1895 
1896  // find the local index for the point 'b' on the cell
1897  for (k=0;k<4;k++)
1898  if (Tetrahedrals[4*Neib[0]+k] == b) break;
1899 
1900  // find the local index for the point 'c' on the cell
1901  for (l=0;l<4;l++)
1902  if (Tetrahedrals[4*Neib[0]+l] == c) break;
1903 
1904  l = l*100 + k*10 + j;
1905 
1906 // cout<<""<< l <<endl;
1907 
1908  switch (l) // j will contain the local index for the current face
1909  {
1910  case 210: case 21: case 102:
1911  case 120: case 12: case 201:
1912  j = 0;
1913  break;
1914  case 310: case 31: case 103:
1915  case 130: case 13: case 301:
1916  j = 1;
1917  break;
1918  case 321: case 132: case 213:
1919  case 231: case 123: case 312:
1920  j = 2;
1921  break;
1922  case 230: case 23: case 302:
1923  case 320: case 32: case 203:
1924  j = 3;
1925  break;
1926 
1927  default:
1928  Error("Unable to set the face !!!!!!!!!!!!" << endl);
1929  exit(0);
1930  }
1931  CellTree[Neib[0]]->SetJoint(j, Joint);
1932 
1933  if (Neib[1] != -1) // second element containing the current face
1934  {
1935  // find the local index for the point 'a' on the cell
1936  for (j=0;j<4;j++)
1937  if (Tetrahedrals[4*Neib[1]+j] == a) break;
1938 
1939  // find the local index for the point 'b' on the cell
1940  for (k=0;k<4;k++)
1941  if (Tetrahedrals[4*Neib[1]+k] == b) break;
1942 
1943  // find the local index for the point 'c' on the cell
1944  for (l=0;l<4;l++)
1945  if (Tetrahedrals[4*Neib[1]+l] == c) break;
1946 
1947  l = l*100 + k*10 + j;
1948 
1949 // cout<<""<< l <<endl;
1950 
1951  switch (l) // j will contain the local index for the current face
1952  {
1953  case 210: case 21: case 102:
1954  case 120: case 12: case 201:
1955  j = 0;
1956  break;
1957  case 310: case 31: case 103:
1958  case 130: case 13: case 301:
1959  j = 1;
1960  break;
1961  case 321: case 132: case 213:
1962  case 231: case 123: case 312:
1963  j = 2;
1964  break;
1965  case 230: case 23: case 302:
1966  case 320: case 32: case 203:
1967  j = 3;
1968  break;
1969 
1970  default:
1971  Error("Unable to set the face !!!!!!!!!!!!" << endl);
1972  exit(0);
1973  }
1974  CellTree[Neib[1]]->SetJoint(j, Joint);
1975  } // if (Neib[1] != -1)
1976 
1977  if (Joint->GetType() == InterfaceJoint3D ||
1978  Joint->GetType() == IsoInterfaceJoint3D)
1979  {
1980  ((TInterfaceJoint3D*)Joint)->SetMapType();
1981  ((TInterfaceJoint3D*)(Joint))->CheckOrientation();
1982  }
1983  else
1984  if (Joint->GetType() == JointEqN)
1985  Joint->SetMapType();
1986  } // if(cell->GetJoint(jj) == NULL)
1987  } // for (i=0;i<N_Cells;i++)
1988 
1989 
1990 #ifdef _MPI
1991  if(rank==0)
1992 #endif
1993  cout<<"numberoftrifaces after "<<N_G<<endl;
1994 
1995 #ifdef _MPI
1996  if(rank==0)
1997 #endif
1998  {
1999 // In.deinitialize();
2000 // Out.deinitialize();
2001  }
2002 
2003 
2004 
2005 
2006 #ifdef _MPI
2007  if(rank!=0)
2008  {
2009  delete [] Tetrahedrals ;
2010  delete [] Coordinates;
2011 // delete [] Facelist;
2012 // delete [] Facemarkerlist;
2013  }
2014 #endif
2015 
2016  delete [] PointNeighb;
2017 
2018 // #ifdef _MPI
2019 // if(rank==3)
2020 // printf(" TetrameshGen %d \n", N_G );
2021 // MPI_Finalize();
2022 // #endif
2023 // exit(0);
2024 //
2025 
2026 }
2027 
2028 
2030 void TetraGenWithInputCells(TDomain *&Domain)
2031 {
2032  //======================================================================
2033  // Tetgen for grid generation begin
2034  //======================================================================
2035  int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
2036  int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
2037  int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
2038  int *Tetrahedrals, *PointNeighb, dimension;
2039  int jj, N_Points, MaxLen, *Tetrahedrals_loc;
2040  const int *EdgeVertex, *TmpFV, *TmpLen;
2041  double *TetRegionlist;
2042 
2043  double *Coordinates, N_x, N_y, N_z;
2044  double *Vertices;
2045  double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
2046  double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
2047  double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
2048 
2049  tetgenio In, Out, InAddPts;
2050 
2051  TBaseCell **CellTree, **SurfCellTree;
2052  TGridCell **DelCell;
2053  TVertex **VertexDel, **NewVertices, **NewSurfVertices;
2054  TBoundPart *BoundPart;
2055  TBdPlane **UpdateFaceParams;
2056  TJoint *Joint;
2057  TBoundComp3D *bdcomp;
2058  TCollection *coll, *SurfColl;
2059  TBaseCell *cell;
2060  TBdSphere *UpdateParam;
2061  TShapeDesc *ShapeDesc;
2062 
2063  char *SMESH, *NODE, line[100];
2064  SMESH = TDatabase::ParamDB->SMESHFILE;
2065  NODE = TDatabase::ParamDB->GRAPEBASENAME;
2066 
2067 #ifdef _MPI
2068  int rank, size;
2069  MPI_Comm Comm = TDatabase::ParamDB->Comm;
2070 
2071  MPI_Comm_rank(Comm, &rank);
2072  MPI_Comm_size(Comm, &size);
2073 
2074 #endif
2075 
2077 #ifdef _MPI
2078  if(rank==0)
2079 #endif
2080  {
2081 
2082 // In.initialize();
2083 // Out.initialize();
2084 
2085  std::ostringstream opts;
2086  opts << " ";
2087 
2088  opts.seekp(std::ios::beg);
2089 // opts<<'p'; // Tetrahedralize the PLC. Switches are chosen to read a PLC (p)
2090  opts<<'r'; // -r Reconstructs a previously generated mesh
2091 // opts<<'q'<<1.2; // quality mesh generation(q) with a specified quality bound
2092 // opts<<'a'<<100; // maximum volume constraint
2093 // opts<<'a'; // maximum volume constraint
2094 // opts<<'i'; // Inserts a list of additional points into mesh.
2095  opts<<'z'; // numbers all output items starting from zero
2096 // opts<<'d'; // Detect intersections of PLC facets.
2097  opts<<'f'; // Outputs all faces (including non-boundary)
2098  opts<<'e'; // Outputs a list of edges of the triangulation
2099 // opts<<'I'; // Suppresses mesh iteration numbers.
2100  opts<<'C'; // Checks the consistency of the final mesh.
2101 // opts<<'Q'; // Quiet: No terminal output except errors.
2102 // opts<<'g'; // Outputs mesh to .mesh file for viewing by Medit
2103  opts<<'Y'; // Suppresses boundary facets/segments splitting
2104 // opts<<'A'; // output region attribute list
2105 // opts<<'V'; //verbose mode
2106 // opts<<'V';
2107 // opts<<'V';
2108 // opts<<'V';
2109 // opts<<'V';
2110  opts<<ends;
2111 
2112  In.initialize();
2113  Out.initialize();
2114 // InAddPts.initialize();
2115 // cout << " SMESHFILE is " << SMESH << endl;
2117  ReadMeditMeshWithCells(SMESH, In);
2118 // ReadAdditionalNModes(NODE, InAddPts);
2119 // In.load_medit(SMESH, 1);
2120 
2121 
2122 // InAddPts.load_node(NODE);
2123 // cout<< " InAddPts.numberofpoints " << InAddPts.numberofpoints <<endl;
2124 // exit(0);
2125 // for(i=0;i<In.numberofpoints;i++)
2126 // OutPut(i<<" (x, y, z) = "<<
2127 // In.pointlist[3*i]<<' '<<In.pointlist[3*i+1]<<' '<<In.pointlist[3*i+2]<<endl);
2128 // exit(0);
2129  // Calling tetrahedralize function of 3dtetgen mesh generator
2130 
2131 
2132  tetrahedralize((char*)opts.str().c_str(), &In, &Out);
2133 // tetrahedralize((char*)opts.str().c_str(), &In, &Out, &InAddPts);
2134 // exit(0);
2135  // output: coordinates of all vertices
2136 // for(i=0;i<Out.numberofpoints;i++)
2137 // OutPut(" (x, y, z) = "<< Out.pointlist[3*i]<<' '<<Out.pointlist[3*i+1]<<' '<<Out.pointlist[3*i+2]<<endl);
2138 
2139  } // if(rank==0)
2140 
2142  DeleteDomain(Domain);
2143 
2144 
2145 #ifdef _MPI
2146  if(rank==0)
2147 #endif
2148  {
2149  N_RootCells = Out.numberoftetrahedra;
2150  N_G = Out.numberofpoints;
2151 
2152  // allocate auxillary fields
2153  Coordinates = Out.pointlist;
2154  Tetrahedrals = Out.tetrahedronlist;
2155 
2156  TetRegionlist = Out.tetrahedronattributelist;
2157  }
2158 
2159 #ifdef _MPI
2160  MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
2161  MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
2162 
2163 
2164  if(rank!=0)
2165  {
2166  Coordinates = new double [3*N_G];
2167  Tetrahedrals = new int[4*N_RootCells];
2168  TetRegionlist = new double [N_RootCells];
2169  }
2170 
2171  MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);
2172  MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);
2173  MPI_Bcast(TetRegionlist, N_RootCells, MPI_DOUBLE, 0, Comm);
2174 #endif
2175 
2176  // generate new vertices
2177  NewVertices = new TVertex*[N_G];
2178 
2179  for (i=0;i<N_G;i++)
2180  {
2181  NewVertices[i] = new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);
2182 
2183  // set bounding box
2184  if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
2185  if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
2186  if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
2187  if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
2188  if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
2189  if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
2190  }
2191 
2192  // set bounding box
2193  StartX = Xmin;
2194  StartY = Ymin;
2195  StartZ = Zmin;
2196  BoundX = Xmax - Xmin;
2197  BoundY = Ymax - Ymin;
2198  BoundZ = Zmax - Zmin;
2199 
2200 
2201  Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
2202 // cout<<Xmin <<" "<<Ymin <<" "<<Zmin<<endl;
2203 // cout<<Xmax <<" "<<Ymax <<" "<<Zmax<<endl;
2204 
2205  CellTree = new TBaseCell*[N_RootCells];
2206  // output of each tetraheron vertex indices (four vertices for each)
2207 // for (i=0;i<N_RootCells;i++)
2208 // cout<< Tetrahedrals[4*i]<<" "<<Tetrahedrals[4*i + 1]<<" "
2209 // <<Tetrahedrals[4*i + 2]<<" "<<Tetrahedrals[4*i + 3]<<endl;
2210 
2211 
2212  for (i=0;i<N_RootCells;i++)
2213  {
2214  CellTree[i] = new TMacroCell(TDatabase::RefDescDB[Tetrahedron],
2215  RefLevel);
2216 
2217  CellTree[i]->SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
2218  CellTree[i]->SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
2219  CellTree[i]->SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
2220  CellTree[i]->SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
2221 
2222  CellTree[i]->SetClipBoard(i);
2223  ((TMacroCell *) CellTree[i])->SetSubGridID(0);
2224 
2225 // cout << i<< " TetRegionlist : " << TetRegionlist[i] << endl;
2226  CellTree[i]->SetRegionID((int)TetRegionlist[i]);
2227  }
2228 // exit(0);
2229  Domain->SetTreeInfo(CellTree, N_RootCells);
2230 
2231 
2232  // initialize iterators
2233  TDatabase::IteratorDB[It_EQ]->SetParam(Domain);
2234  TDatabase::IteratorDB[It_LE]->SetParam(Domain);
2235  TDatabase::IteratorDB[It_Finest]->SetParam(Domain);
2236  TDatabase::IteratorDB[It_Between]->SetParam(Domain);
2237  TDatabase::IteratorDB[It_OCAF]->SetParam(Domain);
2238 
2239 
2240  // search neighbours
2241  PointNeighb = new int[N_G];
2242 #ifdef _MPI
2243  if(rank==0)
2244 #endif
2245  cout<<"numberofpoints "<<N_G<<endl;
2246  memset(PointNeighb, 0, N_G *SizeOfInt);
2247 
2248  for (i=0;i<4*N_RootCells;i++)
2249  PointNeighb[Tetrahedrals[i]]++;
2250 
2251  for (i=0;i<N_G;i++)
2252  if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
2253  delete [] PointNeighb;
2254 
2255 #ifdef _MPI
2256  if(rank==0)
2257 #endif
2258  cout<<"maxEpV "<< maxEpV<<endl;
2259 
2260  PointNeighb = new int[++maxEpV * N_G];
2261 
2262  memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);
2263 
2264  // every vertex contains "maxEpV" columns
2265  // for every vertex at first colomn contains the number of cells containing this vertex
2266  // at further columns we set the index of corresponding cells containing this vertex
2267  // cout<<"maxEpV*N_G "<<maxEpV*N_G<<endl;
2268 
2269  for(i=0;i<4*N_RootCells;i++)
2270  {
2271  j = Tetrahedrals[i]*maxEpV;
2272  PointNeighb[j]++;
2273  //cout<<"j + PointNeighb[j] " << j <<endl;
2274  PointNeighb[j + PointNeighb[j]] = i / 4;
2275  }
2276 
2277 
2278  // generate new faces
2279  coll=Domain->GetCollection(It_Finest, 0);
2280  N_Cells = coll->GetN_Cells();
2281  N_G = 0;
2282  for (i=0;i<N_Cells;i++)
2283  {
2284  cell = coll->GetCell(i);
2285  ShapeDesc= cell->GetShapeDesc();
2286  ShapeDesc->GetFaceVertex(TmpFV, TmpLen, MaxLen);
2287  N_Faces = cell->GetN_Faces();
2288  Tetrahedrals_loc = Tetrahedrals+4*i;
2289 
2290  for(jj=0;jj<N_Faces;jj++)
2291  if(cell->GetJoint(jj) == NULL)
2292  {
2293  N_G++;
2294  N_Points = TmpLen[jj];
2295 
2296  if(N_Points!=3)
2297  {
2298  cerr << "Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
2299  exit(-1);
2300  }
2301 
2302 
2303  //printf(" TetrameshGen Null joint \n" );
2304  a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
2305  b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
2306  c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];
2307 
2308 // cout<<" "<< a<<" "<< b<<" "<< c<<endl;
2309 
2310  Neib[0] = -1;
2311  Neib[1] = -1;
2312  CurrNeib = 0;
2313 
2314  len1 = PointNeighb[a*maxEpV];
2315  len2 = PointNeighb[b*maxEpV];
2316  len3 = PointNeighb[c*maxEpV];
2317 
2318  // find the index of the cells containing current face with point indices a,b,c
2319  for (j=1;j<=len1;j++)
2320  {
2321  Neighb_tmp = PointNeighb[a*maxEpV + j];
2322  for (k=1;k<=len2;k++)
2323  {
2324  if (Neighb_tmp == PointNeighb[b*maxEpV + k])
2325  {
2326  for (l=1;l<=len3;l++)
2327  if (Neighb_tmp == PointNeighb[c*maxEpV + l])
2328  {
2329  Neib[CurrNeib++] = Neighb_tmp;
2330  break;
2331  }
2332  }
2333  }
2334  if (CurrNeib == 2) break;
2335  }// for (j=1;j<=len1;j++)
2336 
2337  // cout<<"CurrNeib " << CurrNeib <<endl;
2338  // cout<<"Facemarkerlist[i] : "<<Facemarkerlist[i]<<endl;
2339  if (CurrNeib == 1) // 0 for inner edges and Boundcomp+1 for Boundedge respect
2340  {
2341  CurrComp = 0; // not yet implemented fully
2342 
2343 
2344 
2345  bdcomp = Domain->GetBdPart(0)->GetBdComp(CurrComp);
2346 
2347 
2348  if(bdcomp->GetTSofXYZ(NewVertices[a]->GetX(), NewVertices[a]->GetY(),
2349  NewVertices[a]->GetY(), T[1], S[1]) ||
2350  bdcomp->GetTSofXYZ(NewVertices[b]->GetX(), NewVertices[b]->GetY(),
2351  NewVertices[b]->GetY(), T[2], S[2]) ||
2352  bdcomp->GetTSofXYZ(NewVertices[c]->GetX(), NewVertices[c]->GetY(),
2353  NewVertices[c]->GetY(), T[3], S[3]) )
2354  {
2355  cerr<<"Error: could not set parameter values"<<endl;
2356  OutPut(NewVertices[a]<<endl);
2357  OutPut(NewVertices[b]<<endl);
2358  OutPut(NewVertices[c]<<endl);
2359  exit(0);
2360  }
2361 
2362  if (CurrNeib == 2)
2363  {
2364  if(bdcomp->IsFreeBoundary())
2365  Joint = new TIsoInterfaceJoint3D(bdcomp, T, S,
2366  CellTree[Neib[0]], CellTree[Neib[1]] );
2367  else
2368  Joint = new TInterfaceJoint3D(bdcomp, T, S,
2369  CellTree[Neib[0]], CellTree[Neib[1]] );
2370  }
2371  else
2372  {
2373  if(bdcomp->IsFreeBoundary())
2374  Joint = new TIsoBoundFace(bdcomp, T, S);
2375  else
2376  Joint = new TBoundFace(bdcomp, T, S);
2377  }
2378  }
2379  else
2380  {
2381  // cout<<"Inner face"<<endl;
2382  if (CurrNeib != 2)
2383  cerr << "Error !!!!!!!! not enough neighbours!" << endl;
2384 
2385  Joint = new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
2386  }
2387 
2388  // First element containing the current face
2389  // find the local index for the point 'a' on the cell
2390  for (j=0;j<4;j++)
2391  if (Tetrahedrals[4*Neib[0]+j] == a) break;
2392 
2393  // find the local index for the point 'b' on the cell
2394  for (k=0;k<4;k++)
2395  if (Tetrahedrals[4*Neib[0]+k] == b) break;
2396 
2397  // find the local index for the point 'c' on the cell
2398  for (l=0;l<4;l++)
2399  if (Tetrahedrals[4*Neib[0]+l] == c) break;
2400 
2401  l = l*100 + k*10 + j;
2402 
2403 // cout<<""<< l <<endl;
2404 
2405  switch (l) // j will contain the local index for the current face
2406  {
2407  case 210: case 21: case 102:
2408  case 120: case 12: case 201:
2409  j = 0;
2410  break;
2411  case 310: case 31: case 103:
2412  case 130: case 13: case 301:
2413  j = 1;
2414  break;
2415  case 321: case 132: case 213:
2416  case 231: case 123: case 312:
2417  j = 2;
2418  break;
2419  case 230: case 23: case 302:
2420  case 320: case 32: case 203:
2421  j = 3;
2422  break;
2423 
2424  default:
2425  Error("Unable to set the face !!!!!!!!!!!!" << endl);
2426  exit(0);
2427  }
2428  CellTree[Neib[0]]->SetJoint(j, Joint);
2429 
2430  if (Neib[1] != -1) // second element containing the current face
2431  {
2432  // find the local index for the point 'a' on the cell
2433  for (j=0;j<4;j++)
2434  if (Tetrahedrals[4*Neib[1]+j] == a) break;
2435 
2436  // find the local index for the point 'b' on the cell
2437  for (k=0;k<4;k++)
2438  if (Tetrahedrals[4*Neib[1]+k] == b) break;
2439 
2440  // find the local index for the point 'c' on the cell
2441  for (l=0;l<4;l++)
2442  if (Tetrahedrals[4*Neib[1]+l] == c) break;
2443 
2444  l = l*100 + k*10 + j;
2445 
2446 // cout<<""<< l <<endl;
2447 
2448  switch (l) // j will contain the local index for the current face
2449  {
2450  case 210: case 21: case 102:
2451  case 120: case 12: case 201:
2452  j = 0;
2453  break;
2454  case 310: case 31: case 103:
2455  case 130: case 13: case 301:
2456  j = 1;
2457  break;
2458  case 321: case 132: case 213:
2459  case 231: case 123: case 312:
2460  j = 2;
2461  break;
2462  case 230: case 23: case 302:
2463  case 320: case 32: case 203:
2464  j = 3;
2465  break;
2466 
2467  default:
2468  Error("Unable to set the face !!!!!!!!!!!!" << endl);
2469  exit(0);
2470  }
2471  CellTree[Neib[1]]->SetJoint(j, Joint);
2472  } // if (Neib[1] != -1)
2473 
2474  if (Joint->GetType() == InterfaceJoint3D ||
2475  Joint->GetType() == IsoInterfaceJoint3D)
2476  {
2477  ((TInterfaceJoint3D*)Joint)->SetMapType();
2478  ((TInterfaceJoint3D*)(Joint))->CheckOrientation();
2479  }
2480  else
2481  if (Joint->GetType() == JointEqN)
2482  Joint->SetMapType();
2483  } // if(cell->GetJoint(jj) == NULL)
2484  } // for (i=0;i<N_Cells;i++)
2485 
2486 
2487 #ifdef _MPI
2488  if(rank==0)
2489 #endif
2490  cout<<"numberoftrifaces after "<<N_G<<endl;
2491 
2492 #ifdef _MPI
2493  if(rank==0)
2494 #endif
2495  {
2496 // In.deinitialize();
2497 // Out.deinitialize();
2498  }
2499 
2500 
2501 
2502 
2503 #ifdef _MPI
2504  if(rank!=0)
2505  {
2506  delete [] Tetrahedrals ;
2507  delete [] Coordinates;
2508  delete [] TetRegionlist;
2509 // delete [] Facemarkerlist;
2510  }
2511 #endif
2512 
2513  delete [] PointNeighb;
2514 
2515 // #ifdef _MPI
2516 // if(rank==3)
2517 // printf(" TetrameshGen %d \n", N_G );
2518 // MPI_Finalize();
2519 // #endif
2520 // exit(0);
2521 //
2522 
2523 }
2524 
2525 
2526 
2527 
TJoint::SetMapType
void SetMapType()
Definition: Joint.C:107
TBdPlane
Definition: BdPlane.h:18
TDomain::GetTreeInfo
void GetTreeInfo(TBaseCell **&celltree, int &N_rootcells)
get tree of cells
Definition: Domain.h:176
TBoundComp3D::GetTSofXYZ
virtual int GetTSofXYZ(double X, double Y, double Z, double &T, double &S)=0
TJoint::GetType
JointType GetType()
Definition: Joint.h:75
TCollection::GetCell
TBaseCell * GetCell(int i) const
return Cell with index i in Cells-array
Definition: Collection.h:50
TBaseCell::GetN_Faces
int GetN_Faces()
return the number of faces of the cell
Definition: BaseCell.h:190
TBoundPart::GetBdComp
TBoundComp2D * GetBdComp(int i)
Definition: BoundPart.h:49
TBaseCell::SetJoint
int SetJoint(int J_i, TJoint *J)
set the pointer to face J_i to J
Definition: BaseCell.h:168
TParaDB::P0
double P0
Definition: Database.h:533
TShapeDesc::GetFaceVertex
int GetFaceVertex(const int *&TmpFV, const int *&TmpLen, int &MaxLen)
Definition: ShapeDesc.h:124
TDomain
contains the boundary description, the virtual cell tree and macro grid
Definition: Domain.h:36
TBaseCell::SetClipBoard
void SetClipBoard(int value)
set value in ClipBoard
Definition: BaseCell.h:263
TIterator::SetParam
int SetParam(TDomain *domain)
Definition: Iterator.C:17
TDatabase::IteratorDB
static TIterator ** IteratorDB
Definition: Database.h:1131
TIsoBoundFace
Definition: IsoBoundFace.h:18
TBaseCell::GetShapeDesc
TShapeDesc * GetShapeDesc() const
return shape descriptor of refinement descriptor
Definition: BaseCell.h:134
TCollection
store cells in an array, used by cell iterators
Definition: Collection.h:18
TBaseCell::GetVertex
virtual TVertex * GetVertex(int Vert_i)=0
return the pointer to vertex with number i
TMacroCell
represent a unit of the macro grid
Definition: MacroCell.h:15
TJointEqN
Definition: JointEqN.h:20
TCollection::GetN_Cells
int GetN_Cells() const
return number of cells
Definition: Collection.h:46
TShapeDesc
Definition: ShapeDesc.h:29
TBaseCell::GetJoint
TJoint * GetJoint(int J_i)
return the pointer to face with number i
Definition: BaseCell.h:175
TParaDB::Comm
MPI_Comm Comm
Definition: Database.h:978
TJoint
Definition: Joint.h:48
TDomain::GetBdPart
TBoundPart * GetBdPart(int i)
get i-th boundary part
Definition: Domain.h:172
TBoundComp3D
Definition: BoundComp3D.h:17
TVertex::GetY
double GetY() const
Definition: Vertex.h:99
TBaseCell::GetN_Vertices
int GetN_Vertices()
return the number of vertices of the cell
Definition: BaseCell.h:179
TBoundPart
Definition: BoundPart.h:21
TDomain::SetTreeInfo
void SetTreeInfo(TBaseCell **celltree, int N_rootcells)
set tree of cells
Definition: Domain.h:183
TBoundComp::IsFreeBoundary
bool IsFreeBoundary() const
Definition: BoundComp.h:57
TDomain::GetCollection
TCollection * GetCollection(Iterators it, int level)
produce a collection with all cells returned by iterator it
Definition: Domain.C:1982
TIsoInterfaceJoint3D
Definition: IsoInterfaceJoint3D.h:19
TBaseCell
information for finite element data structure
Definition: BaseCell.h:25
TVertex
Definition: Vertex.h:19
TDatabase::RefDescDB
static TRefDesc ** RefDescDB
Definition: Database.h:1125
TBaseCell::SetRegionID
void SetRegionID(int val)
set region number to this cell
Definition: BaseCell.h:347
TBoundFace
Definition: BoundFace.h:18
TVertex::GetX
double GetX() const
Definition: Vertex.h:96
TBaseCell::SetVertex
virtual int SetVertex(int Vert_i, TVertex *Vert)=0
set the pointer of vertex Vert_i to Vert
TGridCell
represent geometric information of the cell
Definition: GridCell.h:15
TBaseCell::SetAsLayerCell
void SetAsLayerCell(int val)
set as LayerCell cell
Definition: BaseCell.h:355
TBdSphere
Definition: BdSphere.h:18
TDatabase::ParamDB
static TParamDB * ParamDB
Definition: Database.h:1134
TInterfaceJoint3D
Definition: InterfaceJoint3D.h:18
Generated on Tue May 10 2016 09:07:56 for ParMooN by   doxygen 1.8.5