5 OutPut(
"Example: import3Dmesh.h" << endl);
8 void GridBoundCondition(
double x,
double y,
double z, BoundCond &cond)
13 void SolidBoundCondition(
double x,
double y,
double z, BoundCond &cond)
18 void FluidBoundCondition(
double x,
double y,
double z, BoundCond &cond)
23 void Exact(
double x,
double y,
double z,
double *values)
33 void InitialCondition(
double x,
double y,
double z,
double *values)
40 void BoundCondition(
double x,
double y,
double z, BoundCond &cond)
53 if(x< 240 && y<-72. && z>191)
61 void BoundValue(
double x,
double y,
double z,
double &value)
66 void BilinearCoeffs(
int n_points,
double *X,
double *Y,
double *Z,
double **parameters,
double **coeffs)
74 double rhsfact, alpha, char_L, beam_r, DimlessBeam_r;
75 double xp=0., yp=0., zp=0., SourceCoord, Sourceradius;
80 xp=2.922130000000000e+02;
81 zp=1.583800000000000e+02;
94 Region_ID = coeffs[0][1];
95 DimlessBeam_r = beam_r/char_L;
118 for(i=0;i<n_points;i++)
140 Sourceradius = sqrt( (x-xp)*(x-xp) + (z-zp)*(z-zp));
147 Sourceradius = sqrt( (x-xp)*(x-xp) + (z-zp)*(z-zp));
154 if(Sourceradius<=DimlessBeam_r)
156 coeff[5] = rhsfact*exp(alpha*SourceCoord*char_L);
169 void ReadAdditionalNModes(
char *NODE, tetgenio &InAddPts)
171 int i, N_Vertices, index, N_Atribude, BDMarker;
173 std::ifstream dat(NODE);
177 cerr <<
"cannot open '" << NODE <<
"' for input" << endl;
181 dat.getline (line, 99);
182 dat >> N_Vertices >> N_Atribude >> BDMarker;
183 dat.getline (line, 99);
185 InAddPts.numberofpoints = N_Vertices;
186 InAddPts.pointlist =
new double[3*N_Vertices];
189 for(i=0;i<N_Vertices; i++)
191 dat.getline (line, 99);
192 dat >> index >> InAddPts.pointlist[3*i] >> InAddPts.pointlist[3*i+1] >> InAddPts.pointlist[3*i+2];
193 cout<< i <<
" vert X: " <<InAddPts.pointlist[3*i] <<
" vert Y: " <<InAddPts.pointlist[3*i+1] <<
" vert Z: " <<InAddPts.pointlist[3*i+2] <<endl;
205 void ReadMeditMeshWithCells(
char *SMESH, tetgenio &In)
207 int i, j, k, dimension, N_FVert, N_Faces, N_Vertices, N_CVert, N_Cells;
208 int BDComp_Min=10000000;
215 tetgenio::polygon *P;
217 std::ifstream dat(SMESH);
221 cerr <<
"cannot open '" << SMESH <<
"' for input" << endl;
230 if ( (!strcmp(line,
"Dimension")) || (!strcmp(line,
"dimension")) || (!strcmp(line,
"DIMENSION")))
232 dat.getline (line, 99);
237 dat.getline (line, 99);
242 cerr <<
"dimension: " << dimension << endl;
251 if ( (!strcmp(line,
"Vertices")) || (!strcmp(line,
"vertices")) || (!strcmp(line,
"VERTICES")) )
253 dat.getline (line, 99);
258 dat.getline (line, 99);
261 In.numberofpoints = N_Vertices;
262 In.pointlist =
new double[3*N_Vertices];
265 for(i=0;i<N_Vertices; i++)
267 dat.getline (line, 99);
268 dat >> x >> In.pointlist[3*i+1] >> In.pointlist[3*i+2];
270 In.pointlist[3*i] = x - 3.606596999999999830777142;
273 if(i==1062 || i==916 || i==341 || i==914 || i==162 || i==1063 || i==1061)
274 cout<< i <<
" vert X: " <<In.pointlist[3*i] <<
" vert Y: " <<In.pointlist[3*i+1] <<
" vert Z: " <<In.pointlist[3*i+2] <<endl;
283 if ( (!strcmp(line,
"Triangles")) || (!strcmp(line,
"triangles")) || (!strcmp(line,
"TRIANGLES")) )
286 dat.getline (line, 99);
290 else if ( (!strcmp(line,
"Quadrilaterals")) || (!strcmp(line,
"quadrilaterals")) || (!strcmp(line,
"QUADRILATERALS")) )
293 dat.getline (line, 99);
299 dat.getline (line, 99);
302 In.numberoffacets = N_Faces;
303 In.facetlist =
new tetgenio::facet[In.numberoffacets];
304 In.facetmarkerlist =
new int[In.numberoffacets];
306 for(i=0;i<N_Faces; i++)
308 dat.getline (line, 99);
310 F = &In.facetlist[i];
312 F->numberofpolygons = 1;
313 F->polygonlist =
new tetgenio::polygon[F->numberofpolygons];
314 F->numberofholes = 0;
316 P = &F->polygonlist[0];
318 P->numberofvertices = N_FVert;
319 P->vertexlist =
new int[P->numberofvertices];
321 for(j=0;j<N_FVert;j++)
324 P->vertexlist[j] = k-1;
327 dat >> In.facetmarkerlist[i];
329 if(BDComp_Min > In.facetmarkerlist[i])
330 BDComp_Min=In.facetmarkerlist[i];
337 for(i=0;i<N_Faces; i++)
338 In.facetmarkerlist[i] -= BDComp_Min;
346 if ( (!strcmp(line,
"Tetrahedron")) || (!strcmp(line,
"tetrahedron")) || (!strcmp(line,
"TETRAHEDRON")) )
349 dat.getline (line, 99);
354 dat.getline (line, 99);
357 In.numberoftetrahedra = N_Cells;
358 In.numberofcorners = N_CVert;
359 In.numberoftetrahedronattributes = 1;
360 In.tetrahedronlist =
new int[N_Cells * 4];
361 In.tetrahedronattributelist =
new REAL[N_Cells];
362 In.tetrahedronvolumelist =
new REAL[N_Cells];
364 for(i=0; i<N_Cells; i++)
366 dat.getline (line, 99);
367 plist = &(In.tetrahedronlist[i * 4]);
369 for(j=0;j<N_CVert;j++)
376 dat >> In.tetrahedronattributelist[i];
389 void ReadMeditMesh(
char *SMESH, tetgenio &In)
391 int i, j, k, dimension, N_FVert, N_Faces, N_Vertices;
392 int BDComp_Min=10000000;
396 tetgenio::polygon *P;
398 std::ifstream dat(SMESH);
402 cerr <<
"cannot open '" << SMESH <<
"' for input" << endl;
411 if ( (!strcmp(line,
"Dimension")) || (!strcmp(line,
"dimension")) || (!strcmp(line,
"DIMENSION")))
413 dat.getline (line, 99);
418 dat.getline (line, 99);
423 cerr <<
"dimension: " << dimension << endl;
432 if ( (!strcmp(line,
"Vertices")) || (!strcmp(line,
"vertices")) || (!strcmp(line,
"VERTICES")) )
434 dat.getline (line, 99);
439 dat.getline (line, 99);
442 In.numberofpoints = N_Vertices;
443 In.pointlist =
new double[3*N_Vertices];
446 for(i=0;i<N_Vertices; i++)
448 dat.getline (line, 99);
449 dat >> In.pointlist[3*i] >> In.pointlist[3*i+1] >> In.pointlist[3*i+2];
459 if ( (!strcmp(line,
"Triangles")) || (!strcmp(line,
"triangles")) || (!strcmp(line,
"TRIANGLES")) )
462 dat.getline (line, 99);
466 else if ( (!strcmp(line,
"Quadrilaterals")) || (!strcmp(line,
"quadrilaterals")) || (!strcmp(line,
"QUADRILATERALS")) )
469 dat.getline (line, 99);
475 dat.getline (line, 99);
478 In.numberoffacets = N_Faces;
479 In.facetlist =
new tetgenio::facet[In.numberoffacets];
480 In.facetmarkerlist =
new int[In.numberoffacets];
482 for(i=0;i<N_Faces; i++)
484 dat.getline (line, 99);
486 F = &In.facetlist[i];
488 F->numberofpolygons = 1;
489 F->polygonlist =
new tetgenio::polygon[F->numberofpolygons];
490 F->numberofholes = 0;
492 P = &F->polygonlist[0];
494 P->numberofvertices = N_FVert;
495 P->vertexlist =
new int[P->numberofvertices];
497 for(j=0;j<N_FVert;j++)
500 P->vertexlist[j] = k-1;
503 dat >> In.facetmarkerlist[i];
505 if(BDComp_Min > In.facetmarkerlist[i])
506 BDComp_Min=In.facetmarkerlist[i];
513 for(i=0;i<N_Faces; i++)
514 In.facetmarkerlist[i] -= BDComp_Min;
523 void DeleteDomain(
TDomain *&Domain)
525 int i, j, k, N_Cells, N_RootCells;
526 int CurrVertex, N_Faces, N_Vertices, ID;
527 TBaseCell **CellTree, **SurfCellTree, *cell;
536 VertexDel =
new TVertex*[8*N_RootCells];
539 for(i=0;i<N_Cells;i++)
544 for(j=0;j<N_Faces;j++)
548 VertexDel[CurrVertex++] = cell->
GetVertex(j);
553 for(k=0;k<CurrVertex;k++)
560 VertexDel[CurrVertex++] = cell->
GetVertex(j);
565 for(k=0;k<CurrVertex;k++)
566 if(VertexDel[k]==cell->
GetVertex((j+1)%N_Vertices))
572 VertexDel[CurrVertex++] = cell->
GetVertex((j+1)%N_Vertices);
576 for(k=0;k<CurrVertex;k++)
577 if(VertexDel[k]==cell->
GetVertex((j+2)%N_Vertices))
583 VertexDel[CurrVertex++] = cell->
GetVertex((j+2)%N_Vertices);
588 for(k=0;k<CurrVertex;k++)
589 if(VertexDel[k]==cell->
GetVertex((j+4)%N_Vertices))
595 VertexDel[CurrVertex++] = cell->
GetVertex((j+4)%N_Vertices);
601 for(i=0;i<CurrVertex;i++)
604 OutPut(CurrVertex<<
" vertices were deleted"<<endl);
607 for(i=0;i<N_RootCells;i++)
610 OutPut(N_RootCells<<
" cells were deleted"<<endl);
613 void TetrameshCreate(
TDomain *&Domain)
615 int i, j, k, l, dimension, N_Vertices;
616 int N_Faces, N_RootCells;
617 int v1, v2, v3, v4, CellMarker, RefLevel=0, *Tetrahedrals, *PointNeighb, maxEpV, *Tetrahedrals_loc;
618 int a, b, c, Neib[2], Neighb_tmp, CurrNeib, len1, len2, len3, CurrComp=0, N_Points ;
619 int N_Cells, MaxLen, jj, CompID=0;
620 const int *EdgeVertex, *TmpFV, *TmpLen;
621 double X, Y, Z, DispX, DispY, DispZ;
622 double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
623 double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
624 double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
625 char *MESH, line[100];
634 std::ifstream dat(MESH);
638 cerr <<
"cannot open '" << MESH <<
"' for input" << endl;
646 if ( (!strcmp(line,
"Dimension")) || (!strcmp(line,
"dimension")) || (!strcmp(line,
"DIMENSION")))
648 dat.getline (line, 99);
653 dat.getline (line, 99);
658 cerr <<
"dimension: " << dimension << endl;
665 if ( (!strcmp(line,
"Vertices")) || (!strcmp(line,
"vertices")) || (!strcmp(line,
"VERTICES")) )
667 dat.getline (line, 99);
672 dat.getline (line, 99);
675 cout <<
"N_Vertices "<<N_Vertices<<endl;
676 NewVertices =
new TVertex*[N_Vertices];
678 for(i=0;i<N_Vertices; i++)
680 dat.getline (line, 99);
682 NewVertices[i] =
new TVertex(X, Y, Z);
684 if (X > Xmax) Xmax = X;
685 if (X < Xmin) Xmin = X;
686 if (Y > Ymax) Ymax = Y;
687 if (Y < Ymin) Ymin = Y;
688 if (Z > Zmax) Zmax = Z;
689 if (Z < Zmin) Zmin = Z;
695 BoundX = Xmax - Xmin;
696 BoundY = Ymax - Ymin;
697 BoundZ = Zmax - Zmin;
699 Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
705 if ( (!strcmp(line,
"Tetrahedra")) || (!strcmp(line,
"tetrahedra")) || (!strcmp(line,
"TETRAHEDRA")) )
707 dat.getline (line, 99);
712 dat.getline (line, 99);
715 cout <<
"number of root cells "<<N_RootCells <<endl;
717 Tetrahedrals =
new int[4*N_RootCells];
719 for (i=0;i<N_RootCells;i++)
721 dat.getline (line, 99);
722 dat >> v1 >> v2 >> v3 >> v4 >> CellMarker;
723 Tetrahedrals[4*i ] = v1-1;
724 Tetrahedrals[4*i + 1] = v2-1;
725 Tetrahedrals[4*i + 2] = v3-1;
726 Tetrahedrals[4*i + 3] = v4-1;
731 CellTree[i]->
SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
732 CellTree[i]->
SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
733 CellTree[i]->
SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
734 CellTree[i]->
SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
736 ((
TMacroCell *) CellTree[i])->SetSubGridID(0);
749 PointNeighb =
new int[N_Vertices];
750 memset(PointNeighb, 0, N_Vertices*SizeOfInt);
752 for (i=0;i<4*N_RootCells;i++)
753 PointNeighb[Tetrahedrals[i]]++;
756 for (i=0;i<N_Vertices;i++)
757 if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
758 delete [] PointNeighb;
759 cout<<
"maximum edges per vertex "<< maxEpV<<endl;
760 PointNeighb =
new int[++maxEpV * N_Vertices];
761 memset(PointNeighb, 0, maxEpV*N_Vertices*SizeOfInt);
767 for(i=0;i<4*N_RootCells;i++)
769 j = Tetrahedrals[i]*maxEpV;
771 PointNeighb[j + PointNeighb[j]] = i / 4;
777 for(i=0;i<N_Cells;i++)
786 Tetrahedrals_loc = Tetrahedrals+4*i;
787 for(jj=0;jj<N_Faces;jj++)
790 N_Points = TmpLen[jj];
793 cerr <<
"Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
796 a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
797 b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
798 c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];
802 len1 = PointNeighb[a*maxEpV];
803 len2 = PointNeighb[b*maxEpV];
804 len3 = PointNeighb[c*maxEpV];
806 for (j=1;j<=len1;j++)
808 Neighb_tmp = PointNeighb[a*maxEpV + j];
809 for (k=1;k<=len2;k++)
811 if(Neighb_tmp == PointNeighb[b*maxEpV + k])
814 if(Neighb_tmp == PointNeighb[c*maxEpV + l])
816 Neib[CurrNeib++] = Neighb_tmp;
821 if (CurrNeib == 2)
break;
833 cout <<
"chek 11.5 "<< i<<endl;
834 cerr<<
"Error: could not set parameter values"<<endl;
835 OutPut(NewVertices[a]<<endl);
836 OutPut(NewVertices[b]<<endl);
837 OutPut(NewVertices[c]<<endl);
845 cerr <<
"Error !!!!!!!! not enough neighbours!" << endl;
846 Joint =
new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
851 if (Tetrahedrals[4*Neib[0]+j] == a)
break;
854 if (Tetrahedrals[4*Neib[0]+k] == b)
break;
857 if (Tetrahedrals[4*Neib[0]+l] == c)
break;
858 l = l*100 + k*10 + j;
862 case 210:
case 21:
case 102:
863 case 120:
case 12:
case 201:
866 case 310:
case 31:
case 103:
867 case 130:
case 13:
case 301:
870 case 321:
case 132:
case 213:
871 case 231:
case 123:
case 312:
874 case 230:
case 23:
case 302:
875 case 320:
case 32:
case 203:
880 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
883 CellTree[Neib[0]]->
SetJoint(j, Joint);
890 if (Tetrahedrals[4*Neib[1]+j] == a)
break;
894 if (Tetrahedrals[4*Neib[1]+k] == b)
break;
898 if (Tetrahedrals[4*Neib[1]+l] == c)
break;
900 l = l*100 + k*10 + j;
903 case 210:
case 21:
case 102:
904 case 120:
case 12:
case 201:
907 case 310:
case 31:
case 103:
908 case 130:
case 13:
case 301:
911 case 321:
case 132:
case 213:
912 case 231:
case 123:
case 312:
915 case 230:
case 23:
case 302:
916 case 320:
case 32:
case 203:
921 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
924 CellTree[Neib[1]]->
SetJoint(j, Joint);
927 if (Joint->
GetType() == InterfaceJoint3D || Joint->
GetType() == IsoInterfaceJoint3D)
933 if (Joint->
GetType() == JointEqN)
937 delete [] PointNeighb;
938 cout<<
"tetrameshcreatesuccessful " <<
"\n";
941 void TetrameshGen(
TDomain *&Domain)
946 int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
947 int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
948 int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
949 int *Tetrahedrals, *PointNeighb, dimension;
950 int *Facelist, *Facemarkerlist;
952 double *Coordinates, N_x, N_y, N_z;
954 double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
955 double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
956 double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
962 TVertex **VertexDel, **NewVertices, **NewSurfVertices;
971 char *SMESH, line[100];
978 MPI_Comm_rank(Comm, &rank);
979 MPI_Comm_size(Comm, &size);
992 std::ostringstream opts;
995 opts.seekp(std::ios::beg);
1015 ReadMeditMesh(SMESH, In);
1022 tetrahedralize((
char*)opts.str().c_str(), &In, &Out);
1037 VertexDel =
new TVertex*[8*N_RootCells];
1040 for(i=0;i<N_Cells;i++)
1045 for(j=0;j<N_Faces;j++)
1049 VertexDel[CurrVertex++] = cell->
GetVertex(j);
1054 for(k=0;k<CurrVertex;k++)
1061 VertexDel[CurrVertex++] = cell->
GetVertex(j);
1066 for(k=0;k<CurrVertex;k++)
1067 if(VertexDel[k]==cell->
GetVertex((j+1)%N_Vertices))
1073 VertexDel[CurrVertex++] = cell->
GetVertex((j+1)%N_Vertices);
1077 for(k=0;k<CurrVertex;k++)
1078 if(VertexDel[k]==cell->
GetVertex((j+2)%N_Vertices))
1084 VertexDel[CurrVertex++] = cell->
GetVertex((j+2)%N_Vertices);
1089 for(k=0;k<CurrVertex;k++)
1090 if(VertexDel[k]==cell->
GetVertex((j+4)%N_Vertices))
1096 VertexDel[CurrVertex++] = cell->
GetVertex((j+4)%N_Vertices);
1102 for(i=0;i<CurrVertex;i++)
1103 delete VertexDel[i];
1104 delete [] VertexDel;
1105 OutPut(CurrVertex<<
" vertices were deleted"<<endl);
1109 for(i=0;i<N_RootCells;i++)
1112 OutPut(N_RootCells<<
" cells were deleted"<<endl);
1119 N_RootCells = Out.numberoftetrahedra;
1120 N_G = Out.numberofpoints;
1123 Coordinates = Out.pointlist;
1124 Tetrahedrals = Out.tetrahedronlist;
1130 MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
1131 MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
1136 Coordinates =
new double [3*N_G];
1137 Tetrahedrals =
new int[4*N_RootCells];
1140 MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);
1141 MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);
1145 NewVertices =
new TVertex*[N_G];
1149 NewVertices[i] =
new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);
1152 if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
1153 if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
1154 if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
1155 if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
1156 if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
1157 if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
1163 BoundX = Xmax - Xmin;
1164 BoundY = Ymax - Ymin;
1165 BoundZ = Zmax - Zmin;
1167 Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
1178 for (i=0;i<N_RootCells;i++)
1182 CellTree[i]->
SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
1183 CellTree[i]->
SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
1184 CellTree[i]->
SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
1185 CellTree[i]->
SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
1188 ((
TMacroCell *) CellTree[i])->SetSubGridID(0);
1202 PointNeighb =
new int[N_G];
1206 cout<<
"numberofpoints "<<N_G<<endl;
1207 memset(PointNeighb, 0, N_G *SizeOfInt);
1209 for (i=0;i<4*N_RootCells;i++)
1210 PointNeighb[Tetrahedrals[i]]++;
1213 if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
1214 delete [] PointNeighb;
1219 cout<<
"maxEpV "<< maxEpV<<endl;
1221 PointNeighb =
new int[++maxEpV * N_G];
1223 memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);
1230 for(i=0;i<4*N_RootCells;i++)
1232 j = Tetrahedrals[i]*maxEpV;
1235 PointNeighb[j + PointNeighb[j]] = i / 4;
1251 N_G = Out.numberoftrifaces;
1252 Facelist = Out.trifacelist;
1253 Facemarkerlist = Out.trifacemarkerlist;
1257 MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
1261 Facelist =
new int [3*N_G];
1262 Facemarkerlist =
new int[N_G];
1266 MPI_Bcast(Facelist, 3*N_G, MPI_INT, 0, Comm);
1267 MPI_Bcast(Facemarkerlist, N_G, MPI_INT, 0, Comm);
1271 cout<<
"numberoftrifaces "<<N_G<<endl;
1276 b = Facelist[3*i+1];
1277 c = Facelist[3*i+2];
1285 len1 = PointNeighb[a*maxEpV];
1286 len2 = PointNeighb[b*maxEpV];
1287 len3 = PointNeighb[c*maxEpV];
1290 for (j=1;j<=len1;j++)
1292 Neighb_tmp = PointNeighb[a*maxEpV + j];
1293 for (k=1;k<=len2;k++)
1295 if (Neighb_tmp == PointNeighb[b*maxEpV + k])
1297 for (l=1;l<=len3;l++)
1298 if (Neighb_tmp == PointNeighb[c*maxEpV + l])
1300 Neib[CurrNeib++] = Neighb_tmp;
1305 if (CurrNeib == 2)
break;
1309 if (Facemarkerlist[i])
1312 CurrComp = Facemarkerlist[i] - 1;
1323 NewVertices[a]->
GetY(), T[1], S[1]) ||
1325 NewVertices[b]->
GetY(), T[2], S[2]) ||
1327 NewVertices[c]->
GetY(), T[3], S[3]) )
1329 cerr<<
"Error: could not set parameter values"<<endl;
1330 OutPut(NewVertices[a]<<endl);
1331 OutPut(NewVertices[b]<<endl);
1332 OutPut(NewVertices[c]<<endl);
1340 CellTree[Neib[0]], CellTree[Neib[1]] );
1343 CellTree[Neib[0]], CellTree[Neib[1]] );
1357 cerr <<
"Error !!!!!!!! not enough neighbours!" << endl;
1359 Joint =
new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
1366 if (Tetrahedrals[4*Neib[0]+j] == a)
break;
1370 if (Tetrahedrals[4*Neib[0]+k] == b)
break;
1374 if (Tetrahedrals[4*Neib[0]+l] == c)
break;
1376 l = l*100 + k*10 + j;
1382 case 210:
case 21:
case 102:
1383 case 120:
case 12:
case 201:
1386 case 310:
case 31:
case 103:
1387 case 130:
case 13:
case 301:
1390 case 321:
case 132:
case 213:
1391 case 231:
case 123:
case 312:
1394 case 230:
case 23:
case 302:
1395 case 320:
case 32:
case 203:
1400 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
1403 CellTree[Neib[0]]->
SetJoint(j, Joint);
1409 if (Tetrahedrals[4*Neib[1]+j] == a)
break;
1413 if (Tetrahedrals[4*Neib[1]+k] == b)
break;
1417 if (Tetrahedrals[4*Neib[1]+l] == c)
break;
1419 l = l*100 + k*10 + j;
1425 case 210:
case 21:
case 102:
1426 case 120:
case 12:
case 201:
1429 case 310:
case 31:
case 103:
1430 case 130:
case 13:
case 301:
1433 case 321:
case 132:
case 213:
1434 case 231:
case 123:
case 312:
1437 case 230:
case 23:
case 302:
1438 case 320:
case 32:
case 203:
1443 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
1446 CellTree[Neib[1]]->
SetJoint(j, Joint);
1449 if (Joint->
GetType() == InterfaceJoint3D ||
1450 Joint->
GetType() == IsoInterfaceJoint3D)
1456 if (Joint->
GetType() == JointEqN)
1475 delete [] Tetrahedrals ;
1476 delete [] Coordinates;
1478 delete [] Facemarkerlist;
1482 delete [] PointNeighb;
1495 void TetraGen(
TDomain *&Domain)
1500 int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
1501 int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
1502 int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
1503 int *Tetrahedrals, *PointNeighb, dimension;
1504 int jj, N_Points, MaxLen, *Tetrahedrals_loc;
1505 const int *EdgeVertex, *TmpFV, *TmpLen;
1507 double *Coordinates, N_x, N_y, N_z;
1509 double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
1510 double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
1511 double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
1517 TVertex **VertexDel, **NewVertices, **NewSurfVertices;
1527 char *SMESH, line[100];
1533 MPI_Comm_rank(Comm, &rank);
1534 MPI_Comm_size(Comm, &size);
1544 std::ostringstream opts;
1547 opts.seekp(std::ios::beg);
1567 ReadMeditMesh(SMESH, In);
1574 tetrahedralize((
char*)opts.str().c_str(), &In, &Out);
1583 DeleteDomain(Domain);
1589 N_RootCells = Out.numberoftetrahedra;
1590 N_G = Out.numberofpoints;
1592 Coordinates = Out.pointlist;
1593 Tetrahedrals = Out.tetrahedronlist;
1597 MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
1598 MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
1601 Coordinates =
new double [3*N_G];
1602 Tetrahedrals =
new int[4*N_RootCells];
1604 MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);
1605 MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);
1609 NewVertices =
new TVertex*[N_G];
1613 NewVertices[i] =
new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);
1615 if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
1616 if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
1617 if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
1618 if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
1619 if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
1620 if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
1627 BoundX = Xmax - Xmin;
1628 BoundY = Ymax - Ymin;
1629 BoundZ = Zmax - Zmin;
1631 Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
1641 for (i=0;i<N_RootCells;i++)
1645 CellTree[i]->
SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
1646 CellTree[i]->
SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
1647 CellTree[i]->
SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
1648 CellTree[i]->
SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
1650 ((
TMacroCell *) CellTree[i])->SetSubGridID(0);
1661 PointNeighb =
new int[N_G];
1665 cout<<
"numberofpoints "<<N_G<<endl;
1666 memset(PointNeighb, 0, N_G *SizeOfInt);
1668 for (i=0;i<4*N_RootCells;i++)
1669 PointNeighb[Tetrahedrals[i]]++;
1672 if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
1673 delete [] PointNeighb;
1678 cout<<
"maxEpV "<< maxEpV<<endl;
1679 PointNeighb =
new int[++maxEpV * N_G];
1680 memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);
1686 for(i=0;i<4*N_RootCells;i++)
1688 j = Tetrahedrals[i]*maxEpV;
1691 PointNeighb[j + PointNeighb[j]] = i / 4;
1697 for (i=0;i<N_Cells;i++)
1703 Tetrahedrals_loc = Tetrahedrals+4*i;
1705 for(jj=0;jj<N_Faces;jj++)
1709 N_Points = TmpLen[jj];
1713 cerr <<
"Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
1718 a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
1719 b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
1720 c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];
1726 len1 = PointNeighb[a*maxEpV];
1727 len2 = PointNeighb[b*maxEpV];
1728 len3 = PointNeighb[c*maxEpV];
1731 for (j=1;j<=len1;j++)
1733 Neighb_tmp = PointNeighb[a*maxEpV + j];
1734 for (k=1;k<=len2;k++)
1736 if (Neighb_tmp == PointNeighb[b*maxEpV + k])
1738 for (l=1;l<=len3;l++)
1739 if (Neighb_tmp == PointNeighb[c*maxEpV + l])
1741 Neib[CurrNeib++] = Neighb_tmp;
1746 if (CurrNeib == 2)
break;
1757 NewVertices[a]->
GetY(), T[1], S[1]) ||
1759 NewVertices[b]->
GetY(), T[2], S[2]) ||
1761 NewVertices[c]->
GetY(), T[3], S[3]) )
1763 cerr<<
"Error: could not set parameter values"<<endl;
1764 OutPut(NewVertices[a]<<endl);
1765 OutPut(NewVertices[b]<<endl);
1766 OutPut(NewVertices[c]<<endl);
1774 CellTree[Neib[0]], CellTree[Neib[1]] );
1777 CellTree[Neib[0]], CellTree[Neib[1]] );
1791 cerr <<
"Error !!!!!!!! not enough neighbours!" << endl;
1793 Joint =
new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
1799 if (Tetrahedrals[4*Neib[0]+j] == a)
break;
1803 if (Tetrahedrals[4*Neib[0]+k] == b)
break;
1807 if (Tetrahedrals[4*Neib[0]+l] == c)
break;
1809 l = l*100 + k*10 + j;
1815 case 210:
case 21:
case 102:
1816 case 120:
case 12:
case 201:
1819 case 310:
case 31:
case 103:
1820 case 130:
case 13:
case 301:
1823 case 321:
case 132:
case 213:
1824 case 231:
case 123:
case 312:
1827 case 230:
case 23:
case 302:
1828 case 320:
case 32:
case 203:
1833 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
1836 CellTree[Neib[0]]->
SetJoint(j, Joint);
1842 if (Tetrahedrals[4*Neib[1]+j] == a)
break;
1846 if (Tetrahedrals[4*Neib[1]+k] == b)
break;
1850 if (Tetrahedrals[4*Neib[1]+l] == c)
break;
1852 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)
1895 cout<<
"numberoftrifaces after "<<N_G<<endl;
1907 delete [] Tetrahedrals ;
1908 delete [] Coordinates;
1913 delete [] PointNeighb;
1924 void TetraGenWithInputCells(
TDomain *&Domain)
1929 int i, j, k, l, N_Coord, *N_FVerts, N_Faces, *Facets;
1930 int N, N_RootCells, N_Cells, CurrVertex, N_Vertices, ID, N_G, RefLevel=0;
1931 int CurrNeib, len1, len2, len3, maxEpV = 0, a, b, c, Neib[2], Neighb_tmp, CurrComp;
1932 int *Tetrahedrals, *PointNeighb, dimension;
1933 int jj, N_Points, MaxLen, *Tetrahedrals_loc;
1934 const int *EdgeVertex, *TmpFV, *TmpLen;
1935 double *TetRegionlist;
1937 double *Coordinates, N_x, N_y, N_z;
1939 double Xmin = 1e10, Xmax = -1e10, Ymin = 1e10, Ymax = -1e10;
1940 double Zmin = 1e10, Zmax = -1e10, T[4]={0,0,0,0}, S[4]={0,0,0,0};
1941 double StartX, StartY, StartZ, BoundX, BoundY, BoundZ;
1943 tetgenio In, Out, InAddPts;
1947 TVertex **VertexDel, **NewVertices, **NewSurfVertices;
1957 char *SMESH, *NODE, line[100];
1964 MPI_Comm_rank(Comm, &rank);
1965 MPI_Comm_size(Comm, &size);
1975 std::ostringstream opts;
1978 opts.seekp(std::ios::beg);
2007 ReadMeditMeshWithCells(SMESH, In);
2019 tetrahedralize((
char*)opts.str().c_str(), &In, &Out);
2029 DeleteDomain(Domain);
2035 N_RootCells = Out.numberoftetrahedra;
2036 N_G = Out.numberofpoints;
2038 Coordinates = Out.pointlist;
2039 Tetrahedrals = Out.tetrahedronlist;
2040 TetRegionlist = Out.tetrahedronattributelist;
2044 MPI_Bcast(&N_RootCells, 1, MPI_INT, 0, Comm);
2045 MPI_Bcast(&N_G, 1, MPI_INT, 0, Comm);
2049 Coordinates =
new double [3*N_G];
2050 Tetrahedrals =
new int[4*N_RootCells];
2051 TetRegionlist =
new double [N_RootCells];
2054 MPI_Bcast(Coordinates, 3*N_G, MPI_DOUBLE, 0, Comm);
2055 MPI_Bcast(Tetrahedrals, 4*N_RootCells, MPI_INT, 0, Comm);
2056 MPI_Bcast(TetRegionlist, N_RootCells, MPI_DOUBLE, 0, Comm);
2060 NewVertices =
new TVertex*[N_G];
2064 NewVertices[i] =
new TVertex(Coordinates[3*i], Coordinates[3*i+1], Coordinates[3*i+2]);
2067 if (Coordinates[3*i] > Xmax) Xmax = Coordinates[3*i];
2068 if (Coordinates[3*i] < Xmin) Xmin = Coordinates[3*i];
2069 if (Coordinates[3*i+1] > Ymax) Ymax = Coordinates[3*i+1];
2070 if (Coordinates[3*i+1] < Ymin) Ymin = Coordinates[3*i+1];
2071 if (Coordinates[3*i+2] > Zmax) Zmax = Coordinates[3*i+2];
2072 if (Coordinates[3*i+2] < Zmin) Zmin = Coordinates[3*i+2];
2079 BoundX = Xmax - Xmin;
2080 BoundY = Ymax - Ymin;
2081 BoundZ = Zmax - Zmin;
2082 Domain->SetBoundBox(StartX, StartY, StartZ, BoundX, BoundY, BoundZ);
2092 for (i=0;i<N_RootCells;i++)
2095 CellTree[i]->
SetVertex(0, NewVertices[Tetrahedrals[4*i ]]);
2096 CellTree[i]->
SetVertex(1, NewVertices[Tetrahedrals[4*i + 1]]);
2097 CellTree[i]->
SetVertex(2, NewVertices[Tetrahedrals[4*i + 2]]);
2098 CellTree[i]->
SetVertex(3, NewVertices[Tetrahedrals[4*i + 3]]);
2100 ((
TMacroCell *) CellTree[i])->SetSubGridID(0);
2114 PointNeighb =
new int[N_G];
2118 cout<<
"numberofpoints "<<N_G<<endl;
2119 memset(PointNeighb, 0, N_G *SizeOfInt);
2121 for (i=0;i<4*N_RootCells;i++)
2122 PointNeighb[Tetrahedrals[i]]++;
2125 if (PointNeighb[i] > maxEpV) maxEpV = PointNeighb[i];
2126 delete [] PointNeighb;
2131 cout<<
"maxEpV "<< maxEpV<<endl;
2133 PointNeighb =
new int[++maxEpV * N_G];
2135 memset(PointNeighb, 0, maxEpV*N_G*SizeOfInt);
2141 for(i=0;i<4*N_RootCells;i++)
2143 j = Tetrahedrals[i]*maxEpV;
2146 PointNeighb[j + PointNeighb[j]] = i / 4;
2153 for (i=0;i<N_Cells;i++)
2159 Tetrahedrals_loc = Tetrahedrals+4*i;
2161 for(jj=0;jj<N_Faces;jj++)
2165 N_Points = TmpLen[jj];
2169 cerr <<
"Only Tria faces are allowed!!! N_FVert: " << N_Points << endl;
2174 a = Tetrahedrals_loc[TmpFV[jj*MaxLen]];
2175 b = Tetrahedrals_loc[TmpFV[jj*MaxLen + 1]];
2176 c = Tetrahedrals_loc[TmpFV[jj*MaxLen + 2]];
2184 len1 = PointNeighb[a*maxEpV];
2185 len2 = PointNeighb[b*maxEpV];
2186 len3 = PointNeighb[c*maxEpV];
2189 for (j=1;j<=len1;j++)
2191 Neighb_tmp = PointNeighb[a*maxEpV + j];
2192 for (k=1;k<=len2;k++)
2194 if (Neighb_tmp == PointNeighb[b*maxEpV + k])
2196 for (l=1;l<=len3;l++)
2197 if (Neighb_tmp == PointNeighb[c*maxEpV + l])
2199 Neib[CurrNeib++] = Neighb_tmp;
2204 if (CurrNeib == 2)
break;
2215 NewVertices[a]->
GetY(), T[1], S[1]) ||
2217 NewVertices[b]->
GetY(), T[2], S[2]) ||
2219 NewVertices[c]->
GetY(), T[3], S[3]) )
2221 cerr<<
"Error: could not set parameter values"<<endl;
2222 OutPut(NewVertices[a]<<endl);
2223 OutPut(NewVertices[b]<<endl);
2224 OutPut(NewVertices[c]<<endl);
2232 CellTree[Neib[0]], CellTree[Neib[1]] );
2235 CellTree[Neib[0]], CellTree[Neib[1]] );
2249 cerr <<
"Error !!!!!!!! not enough neighbours!" << endl;
2251 Joint =
new TJointEqN(CellTree[Neib[0]], CellTree[Neib[1]]);
2257 if (Tetrahedrals[4*Neib[0]+j] == a)
break;
2261 if (Tetrahedrals[4*Neib[0]+k] == b)
break;
2265 if (Tetrahedrals[4*Neib[0]+l] == c)
break;
2267 l = l*100 + k*10 + j;
2272 case 210:
case 21:
case 102:
2273 case 120:
case 12:
case 201:
2276 case 310:
case 31:
case 103:
2277 case 130:
case 13:
case 301:
2280 case 321:
case 132:
case 213:
2281 case 231:
case 123:
case 312:
2284 case 230:
case 23:
case 302:
2285 case 320:
case 32:
case 203:
2290 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
2293 CellTree[Neib[0]]->
SetJoint(j, Joint);
2299 if (Tetrahedrals[4*Neib[1]+j] == a)
break;
2303 if (Tetrahedrals[4*Neib[1]+k] == b)
break;
2307 if (Tetrahedrals[4*Neib[1]+l] == c)
break;
2309 l = l*100 + k*10 + j;
2315 case 210:
case 21:
case 102:
2316 case 120:
case 12:
case 201:
2319 case 310:
case 31:
case 103:
2320 case 130:
case 13:
case 301:
2323 case 321:
case 132:
case 213:
2324 case 231:
case 123:
case 312:
2327 case 230:
case 23:
case 302:
2328 case 320:
case 32:
case 203:
2333 Error(
"Unable to set the face !!!!!!!!!!!!" << endl);
2336 CellTree[Neib[1]]->
SetJoint(j, Joint);
2339 if (Joint->
GetType() == InterfaceJoint3D ||
2340 Joint->
GetType() == IsoInterfaceJoint3D)
2346 if (Joint->
GetType() == JointEqN)
2354 cout<<
"numberoftrifaces after "<<N_G<<endl;
2367 delete [] Tetrahedrals ;
2368 delete [] Coordinates;
2369 delete [] TetRegionlist;
2374 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