12 OutPut(
"Example: import3Dmesh.h" << endl);
15 void GridBoundCondition(
double x,
double y,
double z, BoundCond &cond)
21 void SolidBoundCondition(
double x,
double y,
double z, BoundCond &cond)
26 void FluidBoundCondition(
double x,
double y,
double z, BoundCond &cond)
31 void Exact(
double x,
double y,
double z,
double *values)
41 void InitialCondition(
double x,
double y,
double z,
double *values)
48 void BoundCondition(
double x,
double y,
double z, BoundCond &cond)
61 if(x< 240 && y<-72. && z>191)
69 void BoundValue(
double x,
double y,
double z,
double &value)
74 void BilinearCoeffs(
int n_points,
double *X,
double *Y,
double *Z,
double **parameters,
double **coeffs)
82 double rhsfact, alpha, char_L, beam_r, DimlessBeam_r;
83 double xp=0., yp=0., zp=0., SourceCoord, Sourceradius;
88 xp=2.922130000000000e+02;
89 zp=1.583800000000000e+02;
102 Region_ID = coeffs[0][1];
103 DimlessBeam_r = beam_r/char_L;
126 for(i=0;i<n_points;i++)
148 Sourceradius = sqrt( (x-xp)*(x-xp) + (z-zp)*(z-zp));
155 Sourceradius = sqrt( (x-xp)*(x-xp) + (z-zp)*(z-zp));
162 if(Sourceradius<=DimlessBeam_r)
164 coeff[5] = rhsfact*exp(alpha*SourceCoord*char_L);
177 void ReadAdditionalNModes(
char *NODE, tetgenio &InAddPts)
179 int i, N_Vertices, index, N_Atribude, BDMarker;
181 std::ifstream dat(NODE);
185 cerr <<
"cannot open '" << NODE <<
"' for input" << endl;
189 dat.getline (line, 99);
190 dat >> N_Vertices >> N_Atribude >> BDMarker;
191 dat.getline (line, 99);
193 InAddPts.numberofpoints = N_Vertices;
194 InAddPts.pointlist =
new double[3*N_Vertices];
197 for(i=0;i<N_Vertices; i++)
199 dat.getline (line, 99);
200 dat >> index >> InAddPts.pointlist[3*i] >> InAddPts.pointlist[3*i+1] >> InAddPts.pointlist[3*i+2];
201 cout<< i <<
" vert X: " <<InAddPts.pointlist[3*i] <<
" vert Y: " <<InAddPts.pointlist[3*i+1] <<
" vert Z: " <<InAddPts.pointlist[3*i+2] <<endl;
213 void ReadMeditMeshWithCells(
char *SMESH, tetgenio &In)
215 int i, j, k, dimension, N_FVert, N_Faces, N_Vertices, N_CVert, N_Cells;
216 int BDComp_Min=10000000;
223 tetgenio::polygon *P;
225 std::ifstream dat(SMESH);
229 cerr <<
"cannot open '" << SMESH <<
"' for input" << endl;
238 if ( (!strcmp(line,
"Dimension")) || (!strcmp(line,
"dimension")) || (!strcmp(line,
"DIMENSION")))
240 dat.getline (line, 99);
245 dat.getline (line, 99);
250 cerr <<
"dimension: " << dimension << endl;
259 if ( (!strcmp(line,
"Vertices")) || (!strcmp(line,
"vertices")) || (!strcmp(line,
"VERTICES")) )
261 dat.getline (line, 99);
266 dat.getline (line, 99);
269 In.numberofpoints = N_Vertices;
270 In.pointlist =
new double[3*N_Vertices];
273 for(i=0;i<N_Vertices; i++)
275 dat.getline (line, 99);
276 dat >> x >> In.pointlist[3*i+1] >> In.pointlist[3*i+2];
278 In.pointlist[3*i] = x - 3.606596999999999830777142;
281 if(i==1062 || i==916 || i==341 || i==914 || i==162 || i==1063 || i==1061)
282 cout<< i <<
" vert X: " <<In.pointlist[3*i] <<
" vert Y: " <<In.pointlist[3*i+1] <<
" vert Z: " <<In.pointlist[3*i+2] <<endl;
291 if ( (!strcmp(line,
"Triangles")) || (!strcmp(line,
"triangles")) || (!strcmp(line,
"TRIANGLES")) )
294 dat.getline (line, 99);
298 else if ( (!strcmp(line,
"Quadrilaterals")) || (!strcmp(line,
"quadrilaterals")) || (!strcmp(line,
"QUADRILATERALS")) )
301 dat.getline (line, 99);
307 dat.getline (line, 99);
310 In.numberoffacets = N_Faces;
311 In.facetlist =
new tetgenio::facet[In.numberoffacets];
312 In.facetmarkerlist =
new int[In.numberoffacets];
314 for(i=0;i<N_Faces; i++)
316 dat.getline (line, 99);
318 F = &In.facetlist[i];
320 F->numberofpolygons = 1;
321 F->polygonlist =
new tetgenio::polygon[F->numberofpolygons];
322 F->numberofholes = 0;
324 P = &F->polygonlist[0];
326 P->numberofvertices = N_FVert;
327 P->vertexlist =
new int[P->numberofvertices];
329 for(j=0;j<N_FVert;j++)
332 P->vertexlist[j] = k-1;
335 dat >> In.facetmarkerlist[i];
337 if(BDComp_Min > In.facetmarkerlist[i])
338 BDComp_Min=In.facetmarkerlist[i];
345 for(i=0;i<N_Faces; i++)
346 In.facetmarkerlist[i] -= BDComp_Min;
354 if ( (!strcmp(line,
"Tetrahedron")) || (!strcmp(line,
"tetrahedron")) || (!strcmp(line,
"TETRAHEDRON")) )
357 dat.getline (line, 99);
362 dat.getline (line, 99);
365 In.numberoftetrahedra = N_Cells;
366 In.numberofcorners = N_CVert;
367 In.numberoftetrahedronattributes = 1;
368 In.tetrahedronlist =
new int[N_Cells * 4];
369 In.tetrahedronattributelist =
new REAL[N_Cells];
370 In.tetrahedronvolumelist =
new REAL[N_Cells];
372 for(i=0; i<N_Cells; i++)
374 dat.getline (line, 99);
375 plist = &(In.tetrahedronlist[i * 4]);
377 for(j=0;j<N_CVert;j++)
384 dat >> In.tetrahedronattributelist[i];
397 void ReadMeditMesh(
char *SMESH, tetgenio &In)
399 int i, j, k, dimension, N_FVert, N_Faces, N_Vertices;
400 int BDComp_Min=10000000;
404 tetgenio::polygon *P;
406 std::ifstream dat(SMESH);
410 cerr <<
"cannot open '" << SMESH <<
"' for input" << endl;
419 if ( (!strcmp(line,
"Dimension")) || (!strcmp(line,
"dimension")) || (!strcmp(line,
"DIMENSION")))
421 dat.getline (line, 99);
426 dat.getline (line, 99);
431 cerr <<
"dimension: " << dimension << endl;
440 if ( (!strcmp(line,
"Vertices")) || (!strcmp(line,
"vertices")) || (!strcmp(line,
"VERTICES")) )
442 dat.getline (line, 99);
447 dat.getline (line, 99);
450 In.numberofpoints = N_Vertices;
451 In.pointlist =
new double[3*N_Vertices];
454 for(i=0;i<N_Vertices; i++)
456 dat.getline (line, 99);
457 dat >> In.pointlist[3*i] >> In.pointlist[3*i+1] >> In.pointlist[3*i+2];
467 if ( (!strcmp(line,
"Triangles")) || (!strcmp(line,
"triangles")) || (!strcmp(line,
"TRIANGLES")) )
470 dat.getline (line, 99);
474 else if ( (!strcmp(line,
"Quadrilaterals")) || (!strcmp(line,
"quadrilaterals")) || (!strcmp(line,
"QUADRILATERALS")) )
477 dat.getline (line, 99);
483 dat.getline (line, 99);
486 In.numberoffacets = N_Faces;
487 In.facetlist =
new tetgenio::facet[In.numberoffacets];
488 In.facetmarkerlist =
new int[In.numberoffacets];
490 for(i=0;i<N_Faces; i++)
492 dat.getline (line, 99);
494 F = &In.facetlist[i];
496 F->numberofpolygons = 1;
497 F->polygonlist =
new tetgenio::polygon[F->numberofpolygons];
498 F->numberofholes = 0;
500 P = &F->polygonlist[0];
502 P->numberofvertices = N_FVert;
503 P->vertexlist =
new int[P->numberofvertices];
505 for(j=0;j<N_FVert;j++)
508 P->vertexlist[j] = k-1;
511 dat >> In.facetmarkerlist[i];
513 if(BDComp_Min > In.facetmarkerlist[i])
514 BDComp_Min=In.facetmarkerlist[i];
521 for(i=0;i<N_Faces; i++)
522 In.facetmarkerlist[i] -= BDComp_Min;
531 void DeleteDomain(
TDomain *&Domain)
533 int i, j, k, N_Cells, N_RootCells;
534 int CurrVertex, N_Faces, N_Vertices, ID;
535 TBaseCell **CellTree, **SurfCellTree, *cell;
545 VertexDel =
new TVertex*[8*N_RootCells];
548 for(i=0;i<N_Cells;i++)
553 for(j=0;j<N_Faces;j++)
557 VertexDel[CurrVertex++] = cell->
GetVertex(j);
562 for(k=0;k<CurrVertex;k++)
569 VertexDel[CurrVertex++] = cell->
GetVertex(j);
574 for(k=0;k<CurrVertex;k++)
575 if(VertexDel[k]==cell->
GetVertex((j+1)%N_Vertices))
581 VertexDel[CurrVertex++] = cell->
GetVertex((j+1)%N_Vertices);
585 for(k=0;k<CurrVertex;k++)
586 if(VertexDel[k]==cell->
GetVertex((j+2)%N_Vertices))
592 VertexDel[CurrVertex++] = cell->
GetVertex((j+2)%N_Vertices);
597 for(k=0;k<CurrVertex;k++)
598 if(VertexDel[k]==cell->
GetVertex((j+4)%N_Vertices))
604 VertexDel[CurrVertex++] = cell->
GetVertex((j+4)%N_Vertices);
610 for(i=0;i<CurrVertex;i++)
613 OutPut(CurrVertex<<
" vertices were deleted"<<endl);
616 for(i=0;i<N_RootCells;i++)
619 OutPut(N_RootCells<<
" cells were deleted"<<endl);
988 void TetrameshCreate(
TDomain *&Domain)
990 int i, j, k, l, dimension, N_Vertices;
991 int N_FVert, N_Faces, *Facemarkerlist, *Facelist, N_RootCells;
992 int v1, v2, v3, v4, CellMarker, RefLevel=0, *Tetrahedrals, *PointNeighb, maxEpV, *Tetrahedrals_loc;
993 int a, b, c, Neib[2], Neighb_tmp, CurrNeib, len1, len2, len3, CurrComp, N_Points ;
994 int N_Cells, MaxLen, jj, N_SourcePts;
995 const int *EdgeVertex, *TmpFV, *TmpLen;
996 double X, Y, Z, DispX, DispY, DispZ;
997 double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
998 double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
999 double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
1000 char *MESH, line[100];
1012 MPI_Comm_rank(Comm, &rank);
1013 MPI_Comm_size(Comm, &size);
1022 std::ifstream dat(MESH);
1026 cerr <<
"cannot open '" << MESH <<
"' for input" << endl;
1035 if ( (!strcmp(line,
"Dimension")) || (!strcmp(line,
"dimension")) || (!strcmp(line,
"DIMENSION")))
1037 dat.getline (line, 99);
1042 dat.getline (line, 99);
1047 cerr <<
"dimension: " << dimension << endl;
1054 if ( (!strcmp(line,
"Vertices")) || (!strcmp(line,
"vertices")) || (!strcmp(line,
"VERTICES")) )
1056 dat.getline (line, 99);
1061 dat.getline (line, 99);
1064 NewVertices =
new TVertex*[N_Vertices];
1067 for(i=0;i<N_Vertices; i++)
1069 dat.getline (line, 99);
1071 NewVertices[i] =
new TVertex(X, Y, Z);
1073 if (X > Xmax) Xmax = X;
1074 if (X < Xmin) Xmin = X;
1075 if (Y > Ymax) Ymax = Y;
1076 if (Y < Ymin) Ymin = Y;
1077 if (Z > Zmax) Zmax = Z;
1078 if (Z < Zmin) Zmin = Z;
1084 BoundX = Xmax - Xmin;
1085 BoundY = Ymax - Ymin;
1086 BoundZ = Zmax - Zmin;
1088 Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
1094 if ( (!strcmp(line,
"Tetrahedra")) || (!strcmp(line,
"tetrahedra")) || (!strcmp(line,
"TETRAHEDRA")) )
1096 dat.getline (line, 99);
1101 dat.getline (line, 99);
1104 cout <<
"number of root cells "<<N_RootCells <<endl;
1105 if (CellTree == NULL)
1106 cout <<
"no celltree"<<endl;
1108 Tetrahedrals =
new int[4*N_RootCells];
1110 for (i=0;i<N_RootCells;i++)
1112 dat.getline (line, 99);
1113 dat >> v1 >> v2 >> v3 >> v4 >> CellMarker;
1114 Tetrahedrals[4*i ] = v1;
1115 Tetrahedrals[4*i + 1] = v2;
1116 Tetrahedrals[4*i + 2] = v3;
1117 Tetrahedrals[4*i + 3] = v4;
1129 CellTree[i]->
SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
1130 CellTree[i]->
SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
1131 CellTree[i]->
SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
1132 CellTree[i]->
SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
1135 ((
TMacroCell *) CellTree[i])->SetSubGridID(0);
1148 PointNeighb =
new int[N_Vertices];
1152 cout<<
"Number of vertices "<<N_Vertices<<endl;
1153 memset(PointNeighb, 0, N_Vertices*SizeOfInt);
1155 for (i=0;i<4*N_RootCells;i++)
1156 PointNeighb[Tetrahedrals[i]]++;
1159 for (i=0;i<N_Vertices;i++)
1160 if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
1161 delete [] PointNeighb;
1166 cout<<
"max edges per vertex "<< maxEpV<<endl;
1167 PointNeighb =
new int[++maxEpV * N_Vertices];
1168 memset(PointNeighb, 0, maxEpV*N_Vertices*SizeOfInt);
1179 for(i=0;i<4*N_RootCells;i++)
1181 j = Tetrahedrals[i]*maxEpV;
1185 PointNeighb[j + PointNeighb[j]] = i / 4;
1194 for(i=0;i<N_Cells;i++)
1203 Tetrahedrals_loc = Tetrahedrals+4*i;
1205 for(jj=0;jj<N_Faces;jj++)
1208 N_Points = TmpLen[jj];
1211 cerr <<
"Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
1216 a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
1217 b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
1218 c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];
1222 len1 = PointNeighb[a*maxEpV];
1223 len2 = PointNeighb[b*maxEpV];
1224 len3 = PointNeighb[c*maxEpV];
1226 for (j=1;j<=len1;j++)
1228 Neighb_tmp = PointNeighb[a*maxEpV + j];
1229 for (k=1;k<=len2;k++)
1231 if(Neighb_tmp == PointNeighb[b*maxEpV + k])
1233 for(l=1;l<=len3;l++)
1234 if(Neighb_tmp == PointNeighb[c*maxEpV + l])
1236 Neib[CurrNeib++] = Neighb_tmp;
1242 if (CurrNeib == 2)
break;
1248 if(bdcomp->
GetTSofXYZ(NewVertices[a]->
GetX(), NewVertices[a]->
GetY(), NewVertices[a]->
GetZ(), T[1], S[1])
1254 cerr<<
"Error: could not set parameter values"<<endl;
1255 OutPut(NewVertices[a]<<endl);
1256 OutPut(NewVertices[b]<<endl);
1257 OutPut(NewVertices[c]<<endl);
1265 cerr <<
"Error !!!!!!!! not enough neighbours!" << endl;
1266 Joint =
new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
1271 if (Tetrahedrals[4*Neib[0]+j] == a)
break;
1274 if (Tetrahedrals[4*Neib[0]+k] == b)
break;
1277 if (Tetrahedrals[4*Neib[0]+l] == c)
break;
1278 l = l*100 + k*10 + j;
1282 case 210:
case 21:
case 102:
1283 case 120:
case 12:
case 201:
1286 case 310:
case 31:
case 103:
1287 case 130:
case 13:
case 301:
1290 case 321:
case 132:
case 213:
1291 case 231:
case 123:
case 312:
1294 case 230:
case 23:
case 302:
1295 case 320:
case 32:
case 203:
1300 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
1303 CellTree[Neib[0]]->
SetJoint(j, Joint);
1310 if (Tetrahedrals[4*Neib[1]+j] == a)
break;
1314 if (Tetrahedrals[4*Neib[1]+k] == b)
break;
1318 if (Tetrahedrals[4*Neib[1]+l] == c)
break;
1320 l = l*100 + k*10 + j;
1323 case 210:
case 21:
case 102:
1324 case 120:
case 12:
case 201:
1327 case 310:
case 31:
case 103:
1328 case 130:
case 13:
case 301:
1331 case 321:
case 132:
case 213:
1332 case 231:
case 123:
case 312:
1335 case 230:
case 23:
case 302:
1336 case 320:
case 32:
case 203:
1341 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
1344 CellTree[Neib[1]]->
SetJoint(j, Joint);
1347 if (Joint->
GetType() == InterfaceJoint3D || Joint->
GetType() == IsoInterfaceJoint3D)
1353 if (Joint->
GetType() == JointEqN)
1358 delete [] PointNeighb;
1359 cout<<
"Tetrameshcreate successful " <<
"\n";
1362 void TetrameshGen(
TDomain *&Domain)
1367 int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
1368 int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
1369 int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
1370 int *Tetrahedrals, *PointNeighb, dimension;
1371 int *Facelist, *Facemarkerlist;
1373 double *Coordinates, N_x, N_y, N_z;
1375 double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
1376 double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
1377 double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
1383 TVertex **VertexDel, **NewVertices, **NewSurfVertices;
1392 char *SMESH, line[100];
1399 MPI_Comm_rank(Comm, &rank);
1400 MPI_Comm_size(Comm, &size);
1415 std::ostringstream opts;
1418 opts.seekp(std::ios::beg);
1439 ReadMeditMesh(SMESH, In);
1447 tetrahedralize((
char*)opts.str().c_str(), &In, &Out);
1462 VertexDel =
new TVertex*[8*N_RootCells];
1465 for(i=0;i<N_Cells;i++)
1470 for(j=0;j<N_Faces;j++)
1474 VertexDel[CurrVertex++] = cell->
GetVertex(j);
1479 for(k=0;k<CurrVertex;k++)
1486 VertexDel[CurrVertex++] = cell->
GetVertex(j);
1491 for(k=0;k<CurrVertex;k++)
1492 if(VertexDel[k]==cell->
GetVertex((j+1)%N_Vertices))
1498 VertexDel[CurrVertex++] = cell->
GetVertex((j+1)%N_Vertices);
1502 for(k=0;k<CurrVertex;k++)
1503 if(VertexDel[k]==cell->
GetVertex((j+2)%N_Vertices))
1509 VertexDel[CurrVertex++] = cell->
GetVertex((j+2)%N_Vertices);
1514 for(k=0;k<CurrVertex;k++)
1515 if(VertexDel[k]==cell->
GetVertex((j+4)%N_Vertices))
1521 VertexDel[CurrVertex++] = cell->
GetVertex((j+4)%N_Vertices);
1527 for(i=0;i<CurrVertex;i++)
1528 delete VertexDel[i];
1529 delete [] VertexDel;
1530 OutPut(CurrVertex<<
" vertices were deleted"<<endl);
1534 for(i=0;i<N_RootCells;i++)
1537 OutPut(N_RootCells<<
" cells were deleted"<<endl);
1544 N_RootCells = Out.numberoftetrahedra;
1545 N_G = Out.numberofpoints;
1548 Coordinates = Out.pointlist;
1549 Tetrahedrals = Out.tetrahedronlist;
1555 MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
1556 MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
1561 Coordinates =
new double [3*N_G];
1562 Tetrahedrals =
new int[4*N_RootCells];
1565 MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);
1566 MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);
1570 NewVertices =
new TVertex*[N_G];
1574 NewVertices[i] =
new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);
1577 if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
1578 if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
1579 if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
1580 if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
1581 if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
1582 if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
1589 BoundX = Xmax - Xmin;
1590 BoundY = Ymax - Ymin;
1591 BoundZ = Zmax - Zmin;
1594 Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
1605 for (i=0;i<N_RootCells;i++)
1610 CellTree[i]->
SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
1611 CellTree[i]->
SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
1612 CellTree[i]->
SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
1613 CellTree[i]->
SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
1616 ((
TMacroCell *) CellTree[i])->SetSubGridID(0);
1631 PointNeighb =
new int[N_G];
1635 cout<<
"numberofpoints "<<N_G<<endl;
1636 memset(PointNeighb, 0, N_G *SizeOfInt);
1638 for (i=0;i<4*N_RootCells;i++)
1639 PointNeighb[Tetrahedrals[i]]++;
1642 if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
1643 delete [] PointNeighb;
1648 cout<<
"maxEpV "<< maxEpV<<endl;
1650 PointNeighb =
new int[++maxEpV * N_G];
1652 memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);
1659 for(i=0;i<4*N_RootCells;i++)
1661 j = Tetrahedrals[i]*maxEpV;
1664 PointNeighb[j + PointNeighb[j]] = i / 4;
1682 N_G = Out.numberoftrifaces;
1683 Facelist = Out.trifacelist;
1684 Facemarkerlist = Out.trifacemarkerlist;
1688 MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
1692 Facelist =
new int [3*N_G];
1693 Facemarkerlist =
new int[N_G];
1697 MPI_Bcast(Facelist, 3*N_G, MPI_INT, 0, Comm);
1698 MPI_Bcast(Facemarkerlist, N_G, MPI_INT, 0, Comm);
1702 cout<<
"numberoftrifaces "<<N_G<<endl;
1707 b = Facelist[3*i+1];
1708 c = Facelist[3*i+2];
1716 len1 = PointNeighb[a*maxEpV];
1717 len2 = PointNeighb[b*maxEpV];
1718 len3 = PointNeighb[c*maxEpV];
1721 for (j=1;j<=len1;j++)
1723 Neighb_tmp = PointNeighb[a*maxEpV + j];
1724 for (k=1;k<=len2;k++)
1726 if (Neighb_tmp == PointNeighb[b*maxEpV + k])
1728 for (l=1;l<=len3;l++)
1729 if (Neighb_tmp == PointNeighb[c*maxEpV + l])
1731 Neib[CurrNeib++] = Neighb_tmp;
1736 if (CurrNeib == 2)
break;
1740 if (Facemarkerlist[i])
1743 CurrComp = Facemarkerlist[i] - 1;
1754 NewVertices[a]->
GetY(), T[1], S[1]) ||
1756 NewVertices[b]->
GetY(), T[2], S[2]) ||
1758 NewVertices[c]->
GetY(), T[3], S[3]) )
1760 cerr<<
"Error: could not set parameter values"<<endl;
1761 OutPut(NewVertices[a]<<endl);
1762 OutPut(NewVertices[b]<<endl);
1763 OutPut(NewVertices[c]<<endl);
1771 CellTree[Neib[0]], CellTree[Neib[1]] );
1774 CellTree[Neib[0]], CellTree[Neib[1]] );
1788 cerr <<
"Error !!!!!!!! not enough neighbours!" << endl;
1790 Joint =
new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
1797 if (Tetrahedrals[4*Neib[0]+j] == a)
break;
1801 if (Tetrahedrals[4*Neib[0]+k] == b)
break;
1805 if (Tetrahedrals[4*Neib[0]+l] == c)
break;
1807 l = l*100 + k*10 + j;
1813 case 210:
case 21:
case 102:
1814 case 120:
case 12:
case 201:
1817 case 310:
case 31:
case 103:
1818 case 130:
case 13:
case 301:
1821 case 321:
case 132:
case 213:
1822 case 231:
case 123:
case 312:
1825 case 230:
case 23:
case 302:
1826 case 320:
case 32:
case 203:
1831 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
1834 CellTree[Neib[0]]->
SetJoint(j, Joint);
1840 if (Tetrahedrals[4*Neib[1]+j] == a)
break;
1844 if (Tetrahedrals[4*Neib[1]+k] == b)
break;
1848 if (Tetrahedrals[4*Neib[1]+l] == c)
break;
1850 l = l*100 + k*10 + j;
1856 case 210:
case 21:
case 102:
1857 case 120:
case 12:
case 201:
1860 case 310:
case 31:
case 103:
1861 case 130:
case 13:
case 301:
1864 case 321:
case 132:
case 213:
1865 case 231:
case 123:
case 312:
1868 case 230:
case 23:
case 302:
1869 case 320:
case 32:
case 203:
1874 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
1877 CellTree[Neib[1]]->
SetJoint(j, Joint);
1880 if (Joint->
GetType() == InterfaceJoint3D ||
1881 Joint->
GetType() == IsoInterfaceJoint3D)
1887 if (Joint->
GetType() == JointEqN)
1906 delete [] Tetrahedrals ;
1907 delete [] Coordinates;
1909 delete [] Facemarkerlist;
1913 delete [] PointNeighb;
1926 void TetraGen(
TDomain *&Domain)
1931 int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
1932 int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
1933 int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
1934 int *Tetrahedrals, *PointNeighb, dimension;
1935 int jj, N_Points, MaxLen, *Tetrahedrals_loc;
1936 const int *EdgeVertex, *TmpFV, *TmpLen;
1938 double *Coordinates, N_x, N_y, N_z;
1940 double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
1941 double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
1942 double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
1948 TVertex **VertexDel, **NewVertices, **NewSurfVertices;
1958 char *SMESH, line[100];
1964 MPI_Comm_rank(Comm, &rank);
1965 MPI_Comm_size(Comm, &size);
1975 std::ostringstream opts;
1978 opts.seekp(std::ios::beg);
1998 ReadMeditMesh(SMESH, In);
2005 tetrahedralize((
char*)opts.str().c_str(), &In, &Out);
2014 DeleteDomain(Domain);
2020 N_RootCells = Out.numberoftetrahedra;
2021 N_G = Out.numberofpoints;
2023 Coordinates = Out.pointlist;
2024 Tetrahedrals = Out.tetrahedronlist;
2028 MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
2029 MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
2032 Coordinates =
new double [3*N_G];
2033 Tetrahedrals =
new int[4*N_RootCells];
2035 MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);
2036 MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);
2040 NewVertices =
new TVertex*[N_G];
2044 NewVertices[i] =
new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);
2046 if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
2047 if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
2048 if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
2049 if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
2050 if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
2051 if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
2058 BoundX = Xmax - Xmin;
2059 BoundY = Ymax - Ymin;
2060 BoundZ = Zmax - Zmin;
2062 Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
2072 for (i=0;i<N_RootCells;i++)
2076 CellTree[i]->
SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
2077 CellTree[i]->
SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
2078 CellTree[i]->
SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
2079 CellTree[i]->
SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
2081 ((
TMacroCell *) CellTree[i])->SetSubGridID(0);
2092 PointNeighb =
new int[N_G];
2096 cout<<
"numberofpoints "<<N_G<<endl;
2097 memset(PointNeighb, 0, N_G *SizeOfInt);
2099 for (i=0;i<4*N_RootCells;i++)
2100 PointNeighb[Tetrahedrals[i]]++;
2103 if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
2104 delete [] PointNeighb;
2109 cout<<
"maxEpV "<< maxEpV<<endl;
2110 PointNeighb =
new int[++maxEpV * N_G];
2111 memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);
2117 for(i=0;i<4*N_RootCells;i++)
2119 j = Tetrahedrals[i]*maxEpV;
2122 PointNeighb[j + PointNeighb[j]] = i / 4;
2128 for (i=0;i<N_Cells;i++)
2134 Tetrahedrals_loc = Tetrahedrals+4*i;
2136 for(jj=0;jj<N_Faces;jj++)
2140 N_Points = TmpLen[jj];
2144 cerr <<
"Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
2149 a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
2150 b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
2151 c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];
2157 len1 = PointNeighb[a*maxEpV];
2158 len2 = PointNeighb[b*maxEpV];
2159 len3 = PointNeighb[c*maxEpV];
2162 for (j=1;j<=len1;j++)
2164 Neighb_tmp = PointNeighb[a*maxEpV + j];
2165 for (k=1;k<=len2;k++)
2167 if (Neighb_tmp == PointNeighb[b*maxEpV + k])
2169 for (l=1;l<=len3;l++)
2170 if (Neighb_tmp == PointNeighb[c*maxEpV + l])
2172 Neib[CurrNeib++] = Neighb_tmp;
2177 if (CurrNeib == 2)
break;
2188 NewVertices[a]->
GetY(), T[1], S[1]) ||
2190 NewVertices[b]->
GetY(), T[2], S[2]) ||
2192 NewVertices[c]->
GetY(), T[3], S[3]) )
2194 cerr<<
"Error: could not set parameter values"<<endl;
2195 OutPut(NewVertices[a]<<endl);
2196 OutPut(NewVertices[b]<<endl);
2197 OutPut(NewVertices[c]<<endl);
2205 CellTree[Neib[0]], CellTree[Neib[1]] );
2208 CellTree[Neib[0]], CellTree[Neib[1]] );
2222 cerr <<
"Error !!!!!!!! not enough neighbours!" << endl;
2224 Joint =
new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
2230 if (Tetrahedrals[4*Neib[0]+j] == a)
break;
2234 if (Tetrahedrals[4*Neib[0]+k] == b)
break;
2238 if (Tetrahedrals[4*Neib[0]+l] == c)
break;
2240 l = l*100 + k*10 + j;
2246 case 210:
case 21:
case 102:
2247 case 120:
case 12:
case 201:
2250 case 310:
case 31:
case 103:
2251 case 130:
case 13:
case 301:
2254 case 321:
case 132:
case 213:
2255 case 231:
case 123:
case 312:
2258 case 230:
case 23:
case 302:
2259 case 320:
case 32:
case 203:
2264 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
2267 CellTree[Neib[0]]->
SetJoint(j, Joint);
2273 if (Tetrahedrals[4*Neib[1]+j] == a)
break;
2277 if (Tetrahedrals[4*Neib[1]+k] == b)
break;
2281 if (Tetrahedrals[4*Neib[1]+l] == c)
break;
2283 l = l*100 + k*10 + j;
2287 case 210:
case 21:
case 102:
2288 case 120:
case 12:
case 201:
2291 case 310:
case 31:
case 103:
2292 case 130:
case 13:
case 301:
2295 case 321:
case 132:
case 213:
2296 case 231:
case 123:
case 312:
2299 case 230:
case 23:
case 302:
2300 case 320:
case 32:
case 203:
2305 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
2308 CellTree[Neib[1]]->
SetJoint(j, Joint);
2311 if (Joint->
GetType() == InterfaceJoint3D ||
2312 Joint->
GetType() == IsoInterfaceJoint3D)
2318 if (Joint->
GetType() == JointEqN)
2326 cout<<
"numberoftrifaces after "<<N_G<<endl;
2338 delete [] Tetrahedrals ;
2339 delete [] Coordinates;
2344 delete [] PointNeighb;
2355 void TetraGenWithInputCells(
TDomain *&Domain)
2360 int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
2361 int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
2362 int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
2363 int *Tetrahedrals, *PointNeighb, dimension;
2364 int jj, N_Points, MaxLen, *Tetrahedrals_loc;
2365 const int *EdgeVertex, *TmpFV, *TmpLen;
2366 double *TetRegionlist;
2368 double *Coordinates, N_x, N_y, N_z;
2370 double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
2371 double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
2372 double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
2374 tetgenio In, Out, InAddPts;
2378 TVertex **VertexDel, **NewVertices, **NewSurfVertices;
2388 char *SMESH, *NODE, line[100];
2395 MPI_Comm_rank(Comm, &rank);
2396 MPI_Comm_size(Comm, &size);
2406 std::ostringstream opts;
2409 opts.seekp(std::ios::beg);
2438 ReadMeditMeshWithCells(SMESH, In);
2450 tetrahedralize((
char*)opts.str().c_str(), &In, &Out);
2460 DeleteDomain(Domain);
2466 N_RootCells = Out.numberoftetrahedra;
2467 N_G = Out.numberofpoints;
2469 Coordinates = Out.pointlist;
2470 Tetrahedrals = Out.tetrahedronlist;
2471 TetRegionlist = Out.tetrahedronattributelist;
2475 MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
2476 MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
2480 Coordinates =
new double [3*N_G];
2481 Tetrahedrals =
new int[4*N_RootCells];
2482 TetRegionlist =
new double [N_RootCells];
2485 MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);
2486 MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);
2487 MPI_Bcast(TetRegionlist, N_RootCells, MPI_DOUBLE, 0, Comm);
2491 NewVertices =
new TVertex*[N_G];
2495 NewVertices[i] =
new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);
2498 if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
2499 if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
2500 if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
2501 if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
2502 if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
2503 if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
2510 BoundX = Xmax - Xmin;
2511 BoundY = Ymax - Ymin;
2512 BoundZ = Zmax - Zmin;
2513 Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
2523 for (i=0;i<N_RootCells;i++)
2526 CellTree[i]->
SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
2527 CellTree[i]->
SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
2528 CellTree[i]->
SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
2529 CellTree[i]->
SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
2531 ((
TMacroCell *) CellTree[i])->SetSubGridID(0);
2545 PointNeighb =
new int[N_G];
2549 cout<<
"numberofpoints "<<N_G<<endl;
2550 memset(PointNeighb, 0, N_G *SizeOfInt);
2552 for (i=0;i<4*N_RootCells;i++)
2553 PointNeighb[Tetrahedrals[i]]++;
2556 if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
2557 delete [] PointNeighb;
2562 cout<<
"maxEpV "<< maxEpV<<endl;
2564 PointNeighb =
new int[++maxEpV * N_G];
2566 memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);
2572 for(i=0;i<4*N_RootCells;i++)
2574 j = Tetrahedrals[i]*maxEpV;
2577 PointNeighb[j + PointNeighb[j]] = i / 4;
2584 for (i=0;i<N_Cells;i++)
2590 Tetrahedrals_loc = Tetrahedrals+4*i;
2592 for(jj=0;jj<N_Faces;jj++)
2596 N_Points = TmpLen[jj];
2600 cerr <<
"Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
2605 a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
2606 b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
2607 c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];
2615 len1 = PointNeighb[a*maxEpV];
2616 len2 = PointNeighb[b*maxEpV];
2617 len3 = PointNeighb[c*maxEpV];
2620 for (j=1;j<=len1;j++)
2622 Neighb_tmp = PointNeighb[a*maxEpV + j];
2623 for (k=1;k<=len2;k++)
2625 if (Neighb_tmp == PointNeighb[b*maxEpV + k])
2627 for (l=1;l<=len3;l++)
2628 if (Neighb_tmp == PointNeighb[c*maxEpV + l])
2630 Neib[CurrNeib++] = Neighb_tmp;
2635 if (CurrNeib == 2)
break;
2646 NewVertices[a]->
GetY(), T[1], S[1]) ||
2648 NewVertices[b]->
GetY(), T[2], S[2]) ||
2650 NewVertices[c]->
GetY(), T[3], S[3]) )
2652 cerr<<
"Error: could not set parameter values"<<endl;
2653 OutPut(NewVertices[a]<<endl);
2654 OutPut(NewVertices[b]<<endl);
2655 OutPut(NewVertices[c]<<endl);
2663 CellTree[Neib[0]], CellTree[Neib[1]] );
2666 CellTree[Neib[0]], CellTree[Neib[1]] );
2680 cerr <<
"Error !!!!!!!! not enough neighbours!" << endl;
2682 Joint =
new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
2688 if (Tetrahedrals[4*Neib[0]+j] == a)
break;
2692 if (Tetrahedrals[4*Neib[0]+k] == b)
break;
2696 if (Tetrahedrals[4*Neib[0]+l] == c)
break;
2698 l = l*100 + k*10 + j;
2703 case 210:
case 21:
case 102:
2704 case 120:
case 12:
case 201:
2707 case 310:
case 31:
case 103:
2708 case 130:
case 13:
case 301:
2711 case 321:
case 132:
case 213:
2712 case 231:
case 123:
case 312:
2715 case 230:
case 23:
case 302:
2716 case 320:
case 32:
case 203:
2721 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
2724 CellTree[Neib[0]]->
SetJoint(j, Joint);
2730 if (Tetrahedrals[4*Neib[1]+j] == a)
break;
2734 if (Tetrahedrals[4*Neib[1]+k] == b)
break;
2738 if (Tetrahedrals[4*Neib[1]+l] == c)
break;
2740 l = l*100 + k*10 + j;
2746 case 210:
case 21:
case 102:
2747 case 120:
case 12:
case 201:
2750 case 310:
case 31:
case 103:
2751 case 130:
case 13:
case 301:
2754 case 321:
case 132:
case 213:
2755 case 231:
case 123:
case 312:
2758 case 230:
case 23:
case 302:
2759 case 320:
case 32:
case 203:
2764 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
2767 CellTree[Neib[1]]->
SetJoint(j, Joint);
2770 if (Joint->
GetType() == InterfaceJoint3D ||
2771 Joint->
GetType() == IsoInterfaceJoint3D)
2777 if (Joint->
GetType() == JointEqN)
2785 cout<<
"numberoftrifaces after "<<N_G<<endl;
2798 delete [] Tetrahedrals ;
2799 delete [] Coordinates;
2800 delete [] TetRegionlist;
2805 delete [] PointNeighb;
void SetMapType()
Definition: Joint.C:107
void GetTreeInfo(TBaseCell **&celltree, int &N_rootcells)
get tree of cells
Definition: Domain.h:176
virtual int GetTSofXYZ(double X, double Y, double Z, double &T, double &S)=0
JointType GetType()
Definition: Joint.h:75
TBaseCell * GetCell(int i) const
return Cell with index i in Cells-array
Definition: Collection.h:50
int GetN_Faces()
return the number of faces of the cell
Definition: BaseCell.h:190
TBoundComp2D * GetBdComp(int i)
Definition: BoundPart.h:49
int SetJoint(int J_i, TJoint *J)
set the pointer to face J_i to J
Definition: BaseCell.h:168
double P0
Definition: Database.h:533
int GetFaceVertex(const int *&TmpFV, const int *&TmpLen, int &MaxLen)
Definition: ShapeDesc.h:124
contains the boundary description, the virtual cell tree and macro grid
Definition: Domain.h:36
void SetClipBoard(int value)
set value in ClipBoard
Definition: BaseCell.h:263
int SetParam(TDomain *domain)
Definition: Iterator.C:17
static TIterator ** IteratorDB
Definition: Database.h:1131
Definition: IsoBoundFace.h:18
TShapeDesc * GetShapeDesc() const
return shape descriptor of refinement descriptor
Definition: BaseCell.h:134
store cells in an array, used by cell iterators
Definition: Collection.h:18
virtual TVertex * GetVertex(int Vert_i)=0
return the pointer to vertex with number i
represent a unit of the macro grid
Definition: MacroCell.h:15
Definition: JointEqN.h:20
int GetN_Cells() const
return number of cells
Definition: Collection.h:46
Definition: ShapeDesc.h:29
TJoint * GetJoint(int J_i)
return the pointer to face with number i
Definition: BaseCell.h:175
MPI_Comm Comm
Definition: Database.h:978
TBoundPart * GetBdPart(int i)
get i-th boundary part
Definition: Domain.h:172
Definition: BoundComp3D.h:17
double GetY() const
Definition: Vertex.h:99
int GetN_Vertices()
return the number of vertices of the cell
Definition: BaseCell.h:179
Definition: BoundPart.h:21
void SetTreeInfo(TBaseCell **celltree, int N_rootcells)
set tree of cells
Definition: Domain.h:183
bool IsFreeBoundary() const
Definition: BoundComp.h:57
TCollection * GetCollection(Iterators it, int level)
produce a collection with all cells returned by iterator it
Definition: Domain.C:1982
Definition: IsoInterfaceJoint3D.h:19
information for finite element data structure
Definition: BaseCell.h:25
static TRefDesc ** RefDescDB
Definition: Database.h:1125
void SetRegionID(int val)
set region number to this cell
Definition: BaseCell.h:347
Definition: BoundFace.h:18
double GetZ() const
Definition: Vertex.h:103
double GetX() const
Definition: Vertex.h:96
virtual int SetVertex(int Vert_i, TVertex *Vert)=0
set the pointer of vertex Vert_i to Vert
represent geometric information of the cell
Definition: GridCell.h:15
void SetAsLayerCell(int val)
set as LayerCell cell
Definition: BaseCell.h:355
Definition: BdSphere.h:18
static TParamDB * ParamDB
Definition: Database.h:1134
Definition: InterfaceJoint3D.h:18