// // Write JCMwave file // 07.07.2005, Sven Burger, ZIB Berlin // #include <mystdlib.h> #include <myadt.hpp> #include <linalg.hpp> #include <csg.hpp> #include <meshing.hpp> #include <sys/stat.h> namespace netgen { #include "writeuser.hpp" void WriteJCMFormat (const Mesh & mesh, const NetgenGeometry & geom, const string & filename) { if (mesh.GetDimension() != 3) { cout <<"\n Error: Dimension 3 only supported by this output format!"<<endl; return; } int bc_at_infinity = 0; int i, j, jj, ct(0), counter; double dx1, dx2, dx3, dy1, dy2, dy3, dz1, dz2, dz3, vol; // number of points int np = mesh.GetNP(); // Identic points NgArray<int,PointIndex::BASE> identmap1, identmap2, identmap3; mesh.GetIdentifications().GetMap(1, identmap1); mesh.GetIdentifications().GetMap(2, identmap2); mesh.GetIdentifications().GetMap(3, identmap3); // number of volume elements int ne = mesh.GetNE(); int ntets = 0; int nprisms = 0; for (i = 1; i <= ne; i++) { Element el = mesh.VolumeElement(i); if (el.GetNP() == 4) { ntets++; // Check that no two points on a tetrahedron are identified with each other for (j = 1; j <= 4; j++) for (jj = 1; jj <=4; jj++) { if (identmap1.Elem(el.PNum(j)) == el.PNum(jj)) { cout << "\n Error: two points on a tetrahedron identified (1) with each other" << "\n REFINE MESH !" << endl; return; } if (identmap2.Elem(el.PNum(j)) == el.PNum(jj)) { cout << "\n Error: two points on a tetrahedron identified (2) with each other" << "\n REFINE MESH !" << endl; return; } if (identmap3.Elem(el.PNum(j)) == el.PNum(jj)) { cout << "\n Error: two points on a tetrahedron identified (3) with each other" << "\n REFINE MESH !" << endl; return; } } } else if (el.GetNP() == 6) nprisms++; } if ( ne != (ntets+nprisms)) { cout<< "\n Error in determining number of volume elements!\n" << "\n Prisms and tetrahedra only implemented in the JCMwave format!\n"<<endl; return; } if (nprisms > 0) cout << " Please note: Boundaries at infinity have to carry the bc-attribute '-bc=" << bc_at_infinity <<"'."<<endl; // number of surface elements int nse = mesh.GetNSE(); // number of boundary triangles int nbtri = 0; // number of boundary quadrilaterals int nbquad = 0; // array with 1 if point on any tetra, 0 else // this is needed in order to arrange the prism points in the right order NgArray<int,1> pointsOnTetras; pointsOnTetras.SetSize (mesh.GetNP()); pointsOnTetras = 0; for (i = 1; i <= ne; i++) { Element el = mesh.VolumeElement(i); if (el.GetNP() == 4) { for (j = 1; j <= 4; j++) pointsOnTetras.Set(int (el.PNum(j)),1); } } // number of boundary triangles and boundary quadrilaterals for (i = 1; i <= nse; i++) { Element2d el = mesh.SurfaceElement(i); if (el.GetNP() == 3 && ( mesh.GetFaceDescriptor (el.GetIndex()).DomainIn()==0 || mesh.GetFaceDescriptor (el.GetIndex()).DomainOut()==0 ) ) nbtri++; else if (el.GetNP() == 4 && ( mesh.GetFaceDescriptor (el.GetIndex()).DomainIn()==0 || mesh.GetFaceDescriptor (el.GetIndex()).DomainOut()==0 ) ) nbquad++; } ofstream outfile (filename.c_str()); outfile.precision(6); outfile.setf (ios::fixed, ios::floatfield); outfile.setf (ios::showpoint); outfile << "/* <BLOBHead>\n"; outfile << "__BLOBTYPE__=Grid\n"; outfile << "__OWNER__=JCMwave\n"; outfile << "<I>SpaceDim=3\n"; outfile << "<I>ManifoldDim=3\n"; outfile << "<I>NRefinementSteps=0\n"; outfile << "<I>NPoints="<<np<<"\n"; outfile << "<I>NTetrahedra="<<ntets<<"\n"; outfile << "<I>NPrisms="<<nprisms<<"\n"; outfile << "<I>NBoundaryTriangles="<<nbtri<<"\n"; outfile << "<I>NBoundaryQuadrilaterals="<<nbquad<<"\n"; outfile << "*/\n"; outfile << "\n"; outfile << "# output from Netgen\n\n"; int nDomains=mesh.GetNDomains(); for (i=1; i<=nDomains; i++) { if (mesh.GetMaterialPtr(i)) outfile << "#" << mesh.GetMaterial(i) << ": Material ID = " << i << "\n"; } outfile << "# Points\n"; cout << " Please note: The unit of length in the .geo file is assumed to be 'microns'."<<endl; for (i = 1; i <= np; i++) { const Point<3> & p = mesh.Point(i); outfile << i << "\n"; outfile << p(0) << "e-6\n"; outfile << p(1) << "e-6\n"; outfile << p(2) << "e-6\n\n"; } outfile << "\n"; outfile << "# Tetrahedra\n"; counter = 0; for (i = 1; i <= ne; i++) { Element el = mesh.VolumeElement(i); if (el.GetNP() == 4) { counter++; dx1 = mesh.Point(el.PNum(2))(0) - mesh.Point(el.PNum(1))(0); dx2 = mesh.Point(el.PNum(3))(0) - mesh.Point(el.PNum(1))(0); dx3 = mesh.Point(el.PNum(4))(0) - mesh.Point(el.PNum(1))(0); dy1 = mesh.Point(el.PNum(2))(1) - mesh.Point(el.PNum(1))(1); dy2 = mesh.Point(el.PNum(3))(1) - mesh.Point(el.PNum(1))(1); dy3 = mesh.Point(el.PNum(4))(1) - mesh.Point(el.PNum(1))(1); dz1 = mesh.Point(el.PNum(2))(2) - mesh.Point(el.PNum(1))(2); dz2 = mesh.Point(el.PNum(3))(2) - mesh.Point(el.PNum(1))(2); dz3 = mesh.Point(el.PNum(4))(2) - mesh.Point(el.PNum(1))(2); vol = (dy1*dz2-dz1*dy2)*dx3 + (dz1*dx2-dx1*dz2)*dy3 + (dx1*dy2-dy1*dx2)*dz3; if ( vol > 0 ) for (j = 1; j <= 4; j++) outfile << el.PNum(j)<<"\n"; else { for (j = 2; j >= 1; j--) outfile << el.PNum(j)<<"\n"; for (j = 3; j <= 4; j++) outfile << el.PNum(j)<<"\n"; } outfile << el.GetIndex() << "\n\n"; } } if ( counter != ntets) { cout<< "\n Error in determining number of tetras!\n"<<endl; return; } outfile << "\n"; outfile << "# Prisms\n"; counter = 0; for (i = 1; i <= ne; i++) { Element el = mesh.VolumeElement(i); if (el.GetNP() == 6) { counter++; dx1 = mesh.Point(el.PNum(2))(0) - mesh.Point(el.PNum(1))(0); dx2 = mesh.Point(el.PNum(3))(0) - mesh.Point(el.PNum(1))(0); dx3 = mesh.Point(el.PNum(4))(0) - mesh.Point(el.PNum(1))(0); dy1 = mesh.Point(el.PNum(2))(1) - mesh.Point(el.PNum(1))(1); dy2 = mesh.Point(el.PNum(3))(1) - mesh.Point(el.PNum(1))(1); dy3 = mesh.Point(el.PNum(4))(1) - mesh.Point(el.PNum(1))(1); dz1 = mesh.Point(el.PNum(2))(2) - mesh.Point(el.PNum(1))(2); dz2 = mesh.Point(el.PNum(3))(2) - mesh.Point(el.PNum(1))(2); dz3 = mesh.Point(el.PNum(4))(2) - mesh.Point(el.PNum(1))(2); vol = (dy1*dz2-dz1*dy2)*dx3 + (dz1*dx2-dx1*dz2)*dy3 + (dx1*dy2-dy1*dx2)*dz3; if (pointsOnTetras.Get(el.PNum(1)) && pointsOnTetras.Get(el.PNum(2)) && pointsOnTetras.Get(el.PNum(3))) { if (vol > 0) for (j = 1; j <= 6; j++) outfile << el.PNum(j)<<"\n"; else { for (j = 3; j >= 1; j--) outfile << el.PNum(j)<<"\n"; for (j = 6; j >= 4; j--) outfile << el.PNum(j)<<"\n"; } } else if ( pointsOnTetras.Get(el.PNum(4)) && pointsOnTetras.Get(el.PNum(5)) && pointsOnTetras.Get(el.PNum(6)) ) { if ( vol < 0 ) { for (j = 4; j <= 6; j++) outfile << el.PNum(j)<<"\n"; for (j = 1; j <= 3; j++) outfile << el.PNum(j)<<"\n"; } else { for (j = 6; j >= 4; j--) outfile << el.PNum(j)<<"\n"; for (j = 3; j >= 1; j--) outfile << el.PNum(j)<<"\n"; } } else { cout << "\n Error in determining prism point numbering!\n"<<endl; return; } outfile << el.GetIndex() << "\n\n"; } } if ( counter != nprisms) { cout<< "\n Error in determining number of prisms!\n"<<endl; return; } int npid1 = 0; int npid2 = 0; int npid3 = 0; for (i=1; i<=np; i++) { if (identmap1.Elem(i)) npid1++; if (identmap2.Elem(i)) npid2++; if (identmap3.Elem(i)) npid3++; } outfile << "\n"; outfile << "# Boundary triangles\n"; outfile << "# Number of identified points in 1-direction: " << npid1 << "\n"; outfile << "# Number of identified points in 2-direction: " << npid2 << "\n"; outfile << "# Number of identified points in 3-direction: " << npid3 << "\n"; for (i = 1; i <= nse; i++) { Element2d el = mesh.SurfaceElement(i); if (el.GetNP() == 3 && (mesh.GetFaceDescriptor (el.GetIndex()).DomainIn()==0 || mesh.GetFaceDescriptor (el.GetIndex()).DomainOut()==0)) { outfile <<"# T\n"; for (j = 1; j <= 3; j++) outfile << el.PNum(j)<<"\n"; if (mesh.GetFaceDescriptor (el.GetIndex()).BCProperty()==bc_at_infinity) outfile << 1000 << "\n"; else outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << "\n"; if (mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() == bc_at_infinity) outfile << "-2\n\n"; else if (identmap1.Elem(el.PNum(1)) &&identmap1.Elem(el.PNum(2)) &&identmap1.Elem(el.PNum(3))) { outfile << "-1\n"; for (j = 1; j <= 3; j++) outfile << identmap1.Elem(el.PNum(j))<<"\n"; outfile << "\n"; } else if (identmap2.Elem(el.PNum(1)) &&identmap2.Elem(el.PNum(2)) &&identmap2.Elem(el.PNum(3))) { outfile << "-1\n"; for (j = 1; j <= 3; j++) outfile << identmap2.Elem(el.PNum(j))<<"\n"; outfile << "\n"; } else if (identmap3.Elem(el.PNum(1)) &&identmap3.Elem(el.PNum(2)) &&identmap3.Elem(el.PNum(3))) { outfile << "-1\n"; for (j = 1; j <= 3; j++) outfile << identmap3.Elem(el.PNum(j))<<"\n"; outfile << "\n"; } else outfile << "1\n\n"; } } outfile << "\n"; outfile << "# Boundary quadrilaterals\n"; for (i = 1; i <= nse; i++) { Element2d el = mesh.SurfaceElement(i); if (el.GetNP() == 4 && (mesh.GetFaceDescriptor (el.GetIndex()).DomainIn()==0 || mesh.GetFaceDescriptor (el.GetIndex()).DomainOut()==0)) { if (pointsOnTetras.Get(el.PNum(1)) && pointsOnTetras.Get(el.PNum(2))) ct = 0; else if (pointsOnTetras.Get(el.PNum(2)) && pointsOnTetras.Get(el.PNum(3))) ct = 1; else if (pointsOnTetras.Get(el.PNum(3)) && pointsOnTetras.Get(el.PNum(4))) ct = 2; else if (pointsOnTetras.Get(el.PNum(4)) && pointsOnTetras.Get(el.PNum(1))) ct = 3; else cout << "\nWarning: Quadrilateral with inconsistent points found!"<<endl; for (j = 1; j <= 4; j++) { jj = j + ct; if ( jj >= 5 ) jj = jj - 4; outfile << el.PNum(jj)<<"\n"; } outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << "\n"; if (mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() == bc_at_infinity) { outfile << "-2\n\n"; cout << "\nWarning: Quadrilateral at infinity found (this should not occur)!"<<endl; } else if ( identmap1.Elem(el.PNum(1)) && identmap1.Elem(el.PNum(2)) && identmap1.Elem(el.PNum(3)) && identmap1.Elem(el.PNum(4)) ) { outfile << "-1\n"; for (j = 1; j <= 4; j++) { jj = j + ct; if ( jj >= 5 ) jj = jj - 4; outfile << identmap1.Elem(el.PNum(jj))<<"\n"; } outfile << "\n"; } else if ( identmap2.Elem(el.PNum(1)) && identmap2.Elem(el.PNum(2)) && identmap2.Elem(el.PNum(3)) && identmap2.Elem(el.PNum(4)) ) { outfile << "-1\n"; for (j = 1; j <= 4; j++) { jj = j + ct; if ( jj >= 5 ) jj = jj - 4; outfile << identmap2.Elem(el.PNum(jj))<<"\n"; } outfile << "\n"; } else if ( identmap3.Elem(el.PNum(1)) && identmap3.Elem(el.PNum(2)) && identmap3.Elem(el.PNum(3)) && identmap3.Elem(el.PNum(4)) ) { outfile << "-1\n"; for (j = 1; j <= 4; j++) { jj = j + ct; if ( jj >= 5 ) jj = jj - 4; outfile << identmap3.Elem(el.PNum(jj))<<"\n"; } outfile << "\n"; } else outfile << "1\n\n"; } } cout << " JCMwave grid file written." << endl; } }