Adding Gmsh v2.xx Mesh Export Capability - Currently upto 2nd Order Triangles, Quadrangles and Tetrahedra

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
 
 

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 <mystdlib.h>
+
+#include <myadt.hpp>
+#include <linalg.hpp>
+#include <csg.hpp>
+#include <meshing.hpp>
+
+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
+
+

libsrc/interface/writeuser.cpp

@@ -22,7 +22,7 @@ void RegisterUserFormats (Array<const char*> & 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<const char*> & names)
       "STL Format",
       "VRML Format",
       "Gmsh Format",
+      "Gmsh2 Format"
       "JCMwave Format",
       "TET Format",
       //      { "Chemnitz Format" },
@@ -46,10 +47,10 @@ void RegisterUserFormats (Array<const char*> & 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 <sawickib@ee.pw.edu.pl>
 
       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.width(8);
-	  outfile << el.PNum(1)-1; 
+	  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<int> 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);

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<const char*> & names);
 
 extern bool WriteUserFormat (const string & format,
 			     const Mesh & mesh,
-			     const CSGeometry & geom, 
+			     const CSGeometry & geom,
 			     const string & filename);