From af9d554ecb2fc157fc4cfbef78eb40f32039899b Mon Sep 17 00:00:00 2001 From: Philippose Rajan Date: Thu, 29 Jan 2009 20:28:30 +0000 Subject: [PATCH] Adding Gmsh v2.xx Mesh Export Capability - Currently upto 2nd Order Triangles, Quadrangles and Tetrahedra --- libsrc/interface/Makefile.am | 2 +- libsrc/interface/writegmsh2.cpp | 261 ++++++++++++++++++++++++++++++++ libsrc/interface/writeuser.cpp | 138 +++++++++-------- libsrc/interface/writeuser.hpp | 16 +- 4 files changed, 346 insertions(+), 71 deletions(-) create mode 100644 libsrc/interface/writegmsh2.cpp diff --git a/libsrc/interface/Makefile.am b/libsrc/interface/Makefile.am index 9ee9b199..f1675e70 100644 --- a/libsrc/interface/Makefile.am +++ b/libsrc/interface/Makefile.am @@ -6,6 +6,6 @@ noinst_LTLIBRARIES = libinterface.la libinterface_la_SOURCES = read_fnf_mesh.cpp readtetmesh.cpp readuser.cpp writeabaqus.cpp writediffpack.cpp \ writedolfin.cpp writeelmer.cpp writefeap.cpp writefluent.cpp writegmsh.cpp writejcm.cpp \ writepermas.cpp writetecplot.cpp writetet.cpp writetochnog.cpp writeuser.cpp \ - wuchemnitz.cpp + wuchemnitz.cpp writegmsh2.cpp diff --git a/libsrc/interface/writegmsh2.cpp b/libsrc/interface/writegmsh2.cpp new file mode 100644 index 00000000..31452ddd --- /dev/null +++ b/libsrc/interface/writegmsh2.cpp @@ -0,0 +1,261 @@ +/*! \file writegmsh2.cpp + * \brief Export Netgen Mesh in the GMSH v2.xx File format + * \author Philippose Rajan + * \date 02 November 2008 + * + * This function extends the export capabilities of + * Netgen to include the GMSH v2.xx File Format. + * + * Current features of this function include: + * + * 1. Exports Triangles, Quadrangles and Tetrahedra \n + * 2. Supports upto second order elements of each type + * + */ + + +#include + +#include +#include +#include +#include + +namespace netgen +{ +#include "writeuser.hpp" + +// Mapping of entities from Netgen definitions to GMSH definitions +enum GMSH_ELEMENTS {GMSH_TRIG = 2, GMSH_TRIG6 = 9, + GMSH_QUAD = 3, GMSH_QUAD8 = 16, + GMSH_TET = 4, GMSH_TET10 = 11}; +const int triGmsh[7] = {0,1,2,3,6,4,5}; +const int quadGmsh[9] = {0,1,2,3,4,5,8,6,7}; +const int tetGmsh[11] = {0,1,2,3,4,5,8,6,7,10,9}; + + +/*! GMSH v2.xx mesh format export function + * + * This function extends the export capabilities of + * Netgen to include the GMSH v2.xx File Format. + * + * Current features of this function include: + * + * 1. Exports Triangles, Quadrangles and Tetrahedra \n + * 2. Supports upto second order elements of each type + * + */ +void WriteGmsh2Format (const Mesh & mesh, + const CSGeometry & geom, + const string & filename) +{ + ofstream outfile (filename.c_str()); + outfile.precision(6); + outfile.setf (ios::fixed, ios::floatfield); + outfile.setf (ios::showpoint); + + int np = mesh.GetNP(); /// number of points in mesh + int ne = mesh.GetNE(); /// number of 3D elements in mesh + int nse = mesh.GetNSE(); /// number of surface elements (BC) + int i, j, k, l; + + + /* + * 3D section : Volume elements (currently only tetrahedra) + */ + + if ((ne > 0) + && (mesh.VolumeElement(1).GetNP() <= 10) + && (mesh.SurfaceElement(1).GetNP() <= 6)) + { + cout << "Write GMSH v2.xx Format \n"; + cout << "The GMSH v2.xx export is currently available for elements upto 2nd Order\n" << endl; + + int inverttets = mparam.inverttets; + int invertsurf = mparam.inverttrigs; + + /// Prepare GMSH 2.0 file (See GMSH 2.0 Documentation) + outfile << "$MeshFormat\n"; + outfile << (float)2.0 << " " + << (int)0 << " " + << (int)sizeof(double) << "\n"; + outfile << "$EndMeshFormat\n"; + + /// Write nodes + outfile << "$Nodes\n"; + outfile << np << "\n"; + + for (i = 1; i <= np; i++) + { + const Point3d & p = mesh.Point(i); + outfile << i << " "; /// node number + outfile << p.X() << " "; + outfile << p.Y() << " "; + outfile << p.Z() << "\n"; + } + + outfile << "$EndNodes\n"; + + /// write elements (both, surface elements and volume elements) + outfile << "$Elements\n"; + outfile << ne + nse << "\n"; //// number of elements + number of surfaces BC + + for (i = 1; i <= nse; i++) + { + int elType = 0; + + Element2d el = mesh.SurfaceElement(i); + if(invertsurf) el.Invert(); + + if(el.GetNP() == 3) elType = GMSH_TRIG; //// GMSH Type for a 3 node triangle + if(el.GetNP() == 6) elType = GMSH_TRIG6; //// GMSH Type for a 6 node triangle + if(elType == 0) + { + cout << " Invalid surface element type for Gmsh 2.0 3D-Mesh Export Format !\n"; + return; + } + + outfile << i; + outfile << " "; + outfile << elType; + outfile << " "; + outfile << "2"; //// Number of tags (2 => Physical and elementary entities) + outfile << " "; + outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " "; + /// that means that physical entity = elementary entity (arbitrary approach) + outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " "; + for (j = 1; j <= el.GetNP(); j++) + { + outfile << " "; + outfile << el.PNum(triGmsh[j]); + } + outfile << "\n"; + } + + + for (i = 1; i <= ne; i++) + { + int elType = 0; + + Element el = mesh.VolumeElement(i); + if (inverttets) el.Invert(); + + if(el.GetNP() == 4) elType = GMSH_TET; //// GMSH Element type for 4 node tetrahedron + if(el.GetNP() == 10) elType = GMSH_TET10; //// GMSH Element type for 10 node tetrahedron + if(elType == 0) + { + cout << " Invalid volume element type for Gmsh 2.0 3D-Mesh Export Format !\n"; + return; + } + + outfile << nse + i; //// element number (Remember to add on surface elements) + outfile << " "; + outfile << elType; + outfile << " "; + outfile << "2"; //// Number of tags (2 => Physical and elementary entities) + outfile << " "; + outfile << 100000 + el.GetIndex(); + /// that means that physical entity = elementary entity (arbitrary approach) + outfile << " "; + outfile << 100000 + el.GetIndex(); /// volume number + outfile << " "; + for (j = 1; j <= el.GetNP(); j++) + { + outfile << " "; + outfile << el.PNum(tetGmsh[j]); + } + outfile << "\n"; + } + outfile << "$EndElements\n"; + } + /* + * End of 3D section + */ + + + /* + * 2D section : available for triangles and quadrangles + * upto 2nd Order + */ + else if(ne == 0) /// means that there's no 3D element + { + cout << "\n Write Gmsh v2.xx Surface Mesh (triangle and/or quadrangles upto 2nd Order)" << endl; + + /// Prepare GMSH 2.0 file (See GMSH 2.0 Documentation) + outfile << "$MeshFormat\n"; + outfile << (float)2.0 << " " + << (int)0 << " " + << (int)sizeof(double) << "\n"; + outfile << "$EndMeshFormat\n"; + + /// Write nodes + outfile << "$Nodes\n"; + outfile << np << "\n"; + + for (i = 1; i <= np; i++) + { + const Point3d & p = mesh.Point(i); + outfile << i << " "; /// node number + outfile << p.X() << " "; + outfile << p.Y() << " "; + outfile << p.Z() << "\n"; + } + outfile << "$EndNodes\n"; + + /// write triangles & quadrangles + outfile << "$Elements\n"; + outfile << nse << "\n"; + + for (k = 1; k <= nse; k++) + { + int elType = 0; + + const Element2d & el = mesh.SurfaceElement(k); + + if(el.GetNP() == 3) elType = GMSH_TRIG; //// GMSH Type for a 3 node triangle + if(el.GetNP() == 6) elType = GMSH_TRIG6; //// GMSH Type for a 6 node triangle + if(el.GetNP() == 4) elType = GMSH_QUAD; //// GMSH Type for a 4 node quadrangle + if(el.GetNP() == 8) elType = GMSH_QUAD8; //// GMSH Type for an 8 node quadrangle + if(elType == 0) + { + cout << " Invalid surface element type for Gmsh 2.0 2D-Mesh Export Format !\n"; + return; + } + + outfile << k; + outfile << " "; + outfile << elType; + outfile << " "; + outfile << "2"; + outfile << " "; + outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " "; + /// that means that physical entity = elementary entity (arbitrary approach) + outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " "; + for (l = 1; l <= el.GetNP(); l++) + { + outfile << " "; + if((elType == GMSH_TRIG) || (elType == GMSH_TRIG6)) + { + outfile << el.PNum(triGmsh[l]); + } + else if((elType == GMSH_QUAD) || (elType == GMSH_QUAD8)) + { + outfile << el.PNum(quadGmsh[l]); + } + } + outfile << "\n"; + } + outfile << "$EndElements\n"; + } + /* + * End of 2D section + */ + + else + { + cout << " Invalid element type for Gmsh v2.xx Export Format !\n"; + } +} // End: WriteGmsh2Format +} // End: namespace netgen + + diff --git a/libsrc/interface/writeuser.cpp b/libsrc/interface/writeuser.cpp index 77f8ab94..78414fde 100644 --- a/libsrc/interface/writeuser.cpp +++ b/libsrc/interface/writeuser.cpp @@ -22,7 +22,7 @@ void RegisterUserFormats (Array & names) "Neutral Format", "Surface Mesh Format" , "DIFFPACK Format", - "TecPlot Format", + "TecPlot Format", "Tochnog Format", "Abaqus Format", "Fluent Format", @@ -32,6 +32,7 @@ void RegisterUserFormats (Array & names) "STL Format", "VRML Format", "Gmsh Format", + "Gmsh2 Format" "JCMwave Format", "TET Format", // { "Chemnitz Format" }, @@ -46,10 +47,10 @@ void RegisterUserFormats (Array & names) bool WriteUserFormat (const string & format, const Mesh & mesh, - const CSGeometry & geom, + const CSGeometry & geom, const string & filename) { - PrintMessage (1, "Export mesh to file ", filename, + PrintMessage (1, "Export mesh to file ", filename, ", format is ", format); if (format == "Neutral Format") @@ -100,7 +101,12 @@ bool WriteUserFormat (const string & format, else if (format == "Gmsh Format") WriteGmshFormat (mesh, geom, filename); - + + // Philippose - 29/01/2009 + // Added Gmsh v2.xx export capability + else if (format == "Gmsh2 Format") + WriteGmsh2Format (mesh, geom, filename); + else if (format == "JCMwave Format") WriteJCMFormat (mesh, geom, filename); @@ -109,7 +115,7 @@ bool WriteUserFormat (const string & format, WriteTETFormat( mesh, filename);//, "High Frequency" ); #endif - else + else { return 1; } @@ -135,7 +141,7 @@ void WriteNeutralFormat (const Mesh & mesh, int nse = mesh.GetNSE(); int nseg = mesh.GetNSeg(); int i, j; - + int inverttets = mparam.inverttets; int invertsurf = mparam.inverttrigs; @@ -144,13 +150,13 @@ void WriteNeutralFormat (const Mesh & mesh, outfile.precision(6); outfile.setf (ios::fixed, ios::floatfield); outfile.setf (ios::showpoint); - + outfile << np << "\n"; - + for (i = 1; i <= np; i++) { const Point3d & p = mesh.Point(i); - + outfile.width(10); outfile << p.X() << " "; outfile.width(9); @@ -235,13 +241,13 @@ void WriteSurfaceFormat (const Mesh & mesh, { // surface mesh int i, j; - + cout << "Write Surface Mesh" << endl; - + ofstream outfile (filename.c_str()); - + outfile << "surfacemesh" << endl; - + outfile << mesh.GetNP() << endl; for (i = 1; i <= mesh.GetNP(); i++) { @@ -276,37 +282,37 @@ void WriteSTLFormat (const Mesh & mesh, const string & filename) { cout << "\nWrite STL Surface Mesh" << endl; - + ofstream outfile (filename.c_str()); - + int i; - + outfile.precision(10); - + outfile << "solid" << endl; - + for (i = 1; i <= mesh.GetNSE(); i++) { outfile << "facet normal "; const Point3d& p1 = mesh.Point(mesh.SurfaceElement(i).PNum(1)); const Point3d& p2 = mesh.Point(mesh.SurfaceElement(i).PNum(2)); const Point3d& p3 = mesh.Point(mesh.SurfaceElement(i).PNum(3)); - + Vec3d normal = Cross(p2-p1,p3-p1); if (normal.Length() != 0) { - normal /= (normal.Length()); + normal /= (normal.Length()); } - + outfile << normal.X() << " " << normal.Y() << " " << normal.Z() << "\n"; outfile << "outer loop\n"; - + outfile << "vertex " << p1.X() << " " << p1.Y() << " " << p1.Z() << "\n"; outfile << "vertex " << p2.X() << " " << p2.Y() << " " << p2.Z() << "\n"; outfile << "vertex " << p3.X() << " " << p3.Y() << " " << p3.Z() << "\n"; - + outfile << "endloop\n"; - outfile << "endfacet\n"; + outfile << "endfacet\n"; } outfile << "endsolid" << endl; } @@ -329,7 +335,7 @@ void WriteVRMLFormat (const Mesh & mesh, if (faces) { - // Output in VRML, IndexedFaceSet is used + // Output in VRML, IndexedFaceSet is used // Bartosz Sawicki int np = mesh.GetNP(); @@ -351,8 +357,8 @@ void WriteVRMLFormat (const Mesh & mesh, "Shape{ \n" "appearance Appearance { material Material { }} \n" "geometry IndexedFaceSet { \n" - "coord Coordinate { point [ \n"; - + "coord Coordinate { point [ \n"; + for (i = 1; i <= np; i++) { @@ -364,12 +370,12 @@ void WriteVRMLFormat (const Mesh & mesh, } outfile << " ] } \n" - "coordIndex [ \n"; - + "coordIndex [ \n"; + for (i = 1; i <= nse; i++) { const Element2d & el = mesh.SurfaceElement(i); - + for (j = 1; j <= 3; j++) { outfile.width(8); @@ -377,19 +383,19 @@ void WriteVRMLFormat (const Mesh & mesh, } outfile << " -1 \n"; } - + outfile << " ] \n"; //define number and RGB definitions of colors outfile << "color Color { color [1 0 0, 0 1 0, 0 0 1, 1 1 0]} \n" "colorIndex [\n"; - + for (i = 1; i <= nse; i++) { - outfile << mesh.GetFaceDescriptor(mesh.SurfaceElement(i).GetIndex ()).BCProperty(); + outfile << mesh.GetFaceDescriptor(mesh.SurfaceElement(i).GetIndex ()).BCProperty(); outfile << endl; } - + outfile << " ] \n" "colorPerVertex FALSE \n" "creaseAngle 0 \n" @@ -398,7 +404,7 @@ void WriteVRMLFormat (const Mesh & mesh, "convex TRUE \n" "} } # end of Shape\n" "] }\n"; - + } /* end of VRMLFACES */ @@ -427,8 +433,8 @@ void WriteVRMLFormat (const Mesh & mesh, "Shape{ \n" "appearance Appearance { material Material { }} \n" "geometry IndexedLineSet { \n" - "coord Coordinate { point [ \n"; - + "coord Coordinate { point [ \n"; + for (i = 1; i <= np; i++) { @@ -440,40 +446,40 @@ void WriteVRMLFormat (const Mesh & mesh, } outfile << " ] } \n" - "coordIndex [ \n"; - + "coordIndex [ \n"; + for (i = 1; i <= nse; i++) { const Element2d & el = mesh.SurfaceElement(i); - + for (j = 1; j <= 3; j++) { outfile.width(8); outfile << el.PNum(j)-1; } - outfile.width(8); - outfile << el.PNum(1)-1; + outfile.width(8); + outfile << el.PNum(1)-1; outfile << " -1 \n"; } - + outfile << " ] \n"; -/* Uncomment if you want color mesh +/* Uncomment if you want color mesh outfile << "color Color { color [1 1 1, 0 1 0, 0 0 1, 1 1 0]} \n" "colorIndex [\n"; - + for (i = 1; i <= nse; i++) { - outfile << mesh.GetFaceDescriptor(mesh.SurfaceElement(i).GetIndex ()).BCProperty(); + outfile << mesh.GetFaceDescriptor(mesh.SurfaceElement(i).GetIndex ()).BCProperty(); outfile << endl; } - + outfile << " ] \n" -*/ +*/ outfile << "colorPerVertex FALSE \n" "} } #end of Shape\n" "] } \n"; - + } } @@ -490,7 +496,7 @@ void WriteFEPPFormat (const Mesh & mesh, const CSGeometry & geom, const string & filename) { - + ofstream outfile (filename.c_str()); if (mesh.GetDimension() == 3) @@ -498,7 +504,7 @@ void WriteFEPPFormat (const Mesh & mesh, { // output for FEPP - + int np = mesh.GetNP(); int ne = mesh.GetNE(); int nse = mesh.GetNSE(); @@ -508,7 +514,7 @@ void WriteFEPPFormat (const Mesh & mesh, outfile.precision(5); outfile.setf (ios::fixed, ios::floatfield); outfile.setf (ios::showpoint); - + outfile << "volumemesh4" << endl; outfile << nse << endl; for (i = 1; i <= nse; i++) @@ -561,7 +567,7 @@ void WriteFEPPFormat (const Mesh & mesh, outfile << p.Z() << "\n"; } - /* + /* if (typ == WRITE_FEPPML) { int nbn = mesh.mlbetweennodes.Size(); @@ -569,7 +575,7 @@ void WriteFEPPFormat (const Mesh & mesh, for (i = 1; i <= nbn; i++) outfile << mesh.mlbetweennodes.Get(i).I1() << " " << mesh.mlbetweennodes.Get(i).I2() << "\n"; - + // int ncon = mesh.connectedtonode.Size(); // outfile << ncon << "\n"; @@ -586,15 +592,15 @@ void WriteFEPPFormat (const Mesh & mesh, for (i = 1; i <= ns; i++) geom.GetSurface(mesh.GetFaceDescriptor(i).SurfNr())->Print(outfile); } - else + else outfile << "0" << endl; } - + else - + { // 2D fepp format - + ; /* extern SplineGeometry2d * geometry2d; @@ -628,7 +634,7 @@ void WriteEdgeElementFormat (const Mesh & mesh, int nsurfelem = mesh.GetNSE(); int nedges = top->GetNEdges(); int i, j; - + int inverttets = mparam.inverttets; int invertsurf = mparam.inverttrigs; Array edges; @@ -640,12 +646,12 @@ void WriteEdgeElementFormat (const Mesh & mesh, outfile.setf (ios::showpoint); - // vertices with coordinates + // vertices with coordinates outfile << npoints << "\n"; for (i = 1; i <= npoints; i++) { const Point3d & p = mesh.Point(i); - + outfile.width(10); outfile << p.X() << " "; outfile.width(9); @@ -761,7 +767,7 @@ void WriteFile (int typ, double h) { - + int inverttets = mparam.inverttets; int invertsurf = mparam.inverttrigs; @@ -792,7 +798,7 @@ void WriteFile (int typ, // edge (point) on boundary ? BitArray bedge, bpoint(mesh.GetNP()); - + static int eledges[6][2] = { { 1, 2 } , { 1, 3 } , { 1, 4 }, { 2, 3 } , { 2, 4 } , { 3, 4 } }; @@ -815,7 +821,7 @@ void WriteFile (int typ, } } - + // set bedges, bpoints bedge.SetSize (edgelist.Size()); bedge.Clear(); @@ -844,7 +850,7 @@ void WriteFile (int typ, for (i = 1; i <= mesh.GetNE(); i++) { const Element & el = mesh.VolumeElement(i); - + outfile.width(8); outfile << i; for (j = 1; j <= 4; j++) @@ -878,7 +884,7 @@ void WriteFile (int typ, for (i = 1; i <= mesh.GetNP(); i++) { const Point3d & p = mesh.Point(i); - + for (j = 1; j <= 3; j++) { outfile.width(8); diff --git a/libsrc/interface/writeuser.hpp b/libsrc/interface/writeuser.hpp index e11d5370..c841aad0 100644 --- a/libsrc/interface/writeuser.hpp +++ b/libsrc/interface/writeuser.hpp @@ -8,7 +8,7 @@ /**************************************************************************/ -extern +extern void WriteFile (int typ, const Mesh & mesh, const CSGeometry & geom, @@ -18,7 +18,7 @@ void WriteFile (int typ, -extern +extern void ReadFile (Mesh & mesh, const string & filename); @@ -55,6 +55,14 @@ void WriteGmshFormat (const Mesh & mesh, const CSGeometry & geom, const string & filename); + +// Philippose - 29/01/2009 +// Added GMSH v2.xx Mesh Export support +void WriteGmsh2Format (const Mesh & mesh, + const CSGeometry & geom, + const string & filename); + + extern void WriteUserChemnitz (const Mesh & mesh, const string & filename); @@ -65,7 +73,7 @@ void WriteJCMFormat (const Mesh & mesh, const string & filename); -extern +extern void WriteDiffPackFormat (const Mesh & mesh, const CSGeometry & geom, const string & filename); @@ -131,7 +139,7 @@ extern void RegisterUserFormats (Array & names); extern bool WriteUserFormat (const string & format, const Mesh & mesh, - const CSGeometry & geom, + const CSGeometry & geom, const string & filename);