]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added GridOut::write_gmsh.
authorLuca Heltai <luca.heltai@sissa.it>
Fri, 17 Feb 2006 18:09:51 +0000 (18:09 +0000)
committerLuca Heltai <luca.heltai@sissa.it>
Fri, 17 Feb 2006 18:09:51 +0000 (18:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@12399 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_in.h
deal.II/deal.II/include/grid/grid_out.h
deal.II/deal.II/source/grid/grid_generator.cc
deal.II/deal.II/source/grid/grid_out.all_dimensions.cc
deal.II/deal.II/source/grid/grid_out.cc
deal.II/doc/news/changes.html

index b2cff3da8edc47b31a58cdb08fbb3f4e7bc1ced7..333ba455a35e78315366e5a76708eba03ed1d5c4 100644 (file)
@@ -30,7 +30,7 @@ class SubCellData;
 /**
  * This class implements an input mechanism for grid data. It allows
  * to read a grid structure into a triangulation object. At present,
- * only UCD (unstructured cell data) and DB Mesh is supported as input
+ * UCD (unstructured cell data), DB Mesh and Gmsh are supported as input
  * format for grid data. Any numerical data after the block of
  * topological information is ignored.
  *
index 0d515a4429f3ed89df0d022bb18eefb0c4591925..20c33a4a3a5e5ca7789576a7511568e948de8b30 100644 (file)
@@ -77,7 +77,41 @@ namespace GridOutFlags
        const bool write_measure = false,
        const bool write_all_faces = true);  
   };                      
-    
+     
+                                  /**
+                                   * Flags describing the details of
+                                   * output in MSH format.
+                                   *
+                                   * @ingroup output
+                                   */
+  struct Msh 
+  {
+                                      /**
+                                       * When writing a mesh, write boundary
+                                       * faces explicitly if their boundary
+                                       * indicator is not the default
+                                       * boundary indicator, which is zero.
+                                       * This is necessary if you later
+                                       * want to re-read the grid and want
+                                       * to get the same boundary indicators
+                                       * for the different parts of the
+                                       * boundary of the triangulation.
+                                       *
+                                       * It is not necessary if you only want
+                                       * to write the triangulation to
+                                       * view or print it.
+                                       *
+                                       * Default: @p false.
+                                       */
+      bool write_faces;
+      
+                                      /**
+                                       * Constructor.
+                                       */
+      Msh (const bool write_faces    = false);
+  };
+  
                                   /**
                                    * Flags describing the details of
                                    * output in UCD format.
@@ -582,7 +616,7 @@ namespace GridOutFlags
  *
  * @ingroup grid
  * @ingroup output
- * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2003; postscript format based on an implementation by Stefan Nauber, 1999
+ * @author Wolfgang Bangerth, Guido Kanschat, Luca Heltai, 1999, 2003, 2006; postscript format based on an implementation by Stefan Nauber, 1999
  */
 class GridOut 
 {
@@ -608,7 +642,9 @@ class GridOut
                                           /// write() calls write_ucd()
          ucd,
                                           /// write() calls write_xfig()
-         xfig
+         xfig,     
+                                          /// write() calls write_msh()
+         msh
     };
     
                                     /**
@@ -703,6 +739,44 @@ class GridOut
                        std::ostream           &out,
                        const Mapping<1>       *mapping=0);
     
+                                    /**
+                                     * Write the triangulation in the
+                                     * msh format.
+                                     *
+                                     * Msh 
+                                     * is the format used by Gmsh and it is
+                                     * described in the gmsh
+                                     * user's guide. Besides the
+                                     * usual output of the grid
+                                     * only, you can decide through
+                                     * additional flags (see below,
+                                     * and the documentation of the
+                                     * GridOutFlags::Msh() class)
+                                     * whether boundary faces with
+                                     * non-zero boundary indicator
+                                     * shall be written to the file
+                                     * explicitly. This is useful,
+                                     * if you want to re-read the
+                                     * grid later on, since
+                                     * <tt>deal.II</tt> sets the boundary
+                                     * indicator to zero by default;
+                                     * therefore, to obtain the
+                                     * same triangulation as before,
+                                     * you have to specify faces
+                                     * with differing boundary
+                                     * indicators explicitly, which
+                                     * is done by this flag.
+                                     *
+                                     * Names and values of further
+                                     * flags controlling the output
+                                     * can be found in the
+                                     * documentation of the
+                                     * GridOut::Msh() class.
+                                     */
+    template <int dim>
+    void write_msh (const Triangulation<dim> &tria,
+                   std::ostream             &out);
+
                                     /**
                                      * Write the triangulation in the
                                      * ucd format.
@@ -858,6 +932,11 @@ class GridOut
                                      */
     void set_flags (const GridOutFlags::DX &flags);
 
+                                    /**
+                                     * Set flags for GMSH output
+                                     */
+    void set_flags (const GridOutFlags::Msh &flags);
+
                                     /**
                                      * Set flags for UCD output
                                      */
@@ -960,6 +1039,13 @@ class GridOut
                                      */
     GridOutFlags::DX dx_flags;
   
+                                    /**
+                                     * Flags for GMSH output. Can be
+                                     * changed by using the
+                                     * @p set_flags function.
+                                     */
+    GridOutFlags::Msh     msh_flags;
+  
                                     /**
                                      * Flags for UCD output. Can be
                                      * changed by using the
@@ -1002,7 +1088,45 @@ class GridOut
                                      * Flags used for XFig output.
                                      */
     GridOutFlags::XFig xfig_flags;
-    
+     
+                                    /**
+                                     * Write the grid information about
+                                     * faces to @p out. Only those faces
+                                     * are printed which are on the boundary
+                                     * and which have a boundary indicator
+                                     * not equal to zero, since the latter
+                                     * is the default for boundary faces.
+                                     *
+                                     * Since cells and faces are continuously
+                                     * numbered, the @p starting_index for
+                                     * the numbering of the faces is passed
+                                     * also.
+                                     *
+                                     * This function unfortunately can not
+                                     * be included in the regular @p write_msh
+                                     * function, since it needs special
+                                     * treatment for the case <tt>dim==1</tt>, in
+                                     * which case the face iterators are
+                                     * <tt>void*</tt>'s and lack the member functions
+                                     * which are called. We would not actually
+                                     * call these functions, but the compiler
+                                     * would complain anyway when compiling
+                                     * the function for <tt>dim==1</tt>. Bad luck.
+                                     */
+    template <int dim>
+    void write_msh_faces (const Triangulation<dim> &tria,
+                         const unsigned int        starting_index,
+                         std::ostream             &out) const;
+
+                                    /**
+                                     * Declaration of the specialization
+                                     * of above function for 1d. Does
+                                     * nothing.
+                                     */
+    void write_msh_faces (const Triangulation<1> &tria,
+                         const unsigned int      starting_index,
+                         std::ostream           &out) const;
+   
                                     /**
                                      * Write the grid information about
                                      * faces to @p out. Only those faces
index 0b357817d5b85c18b4a1249b4a16207d16999531..f2e46282bef689e76d5e1ca19444aa33038a8b42 100644 (file)
@@ -1394,7 +1394,7 @@ GridGenerator::half_hyper_ball (Triangulation<dim>& tria,
     {
        for (unsigned int i=0;i<GeometryInfo<dim>::faces_per_cell;++i)
        {
-           if (cell->face(i)->boundary_indicator() == 255)
+           if (!cell->at_boundary(i))
                continue;
 
            // If the center is on the plane x=0, this is a planar
@@ -1408,7 +1408,7 @@ GridGenerator::half_hyper_ball (Triangulation<dim>& tria,
        // With this loop we restore back the indicator of the outer lines
        for (unsigned int i=0;i<GeometryInfo<dim>::faces_per_cell;++i)
        {
-           if (cell->face(i)->boundary_indicator() == 255)
+           if (!cell->at_boundary(i))
                continue;
 
            // If the center is not on the plane x=0, this is a curvilinear
index 4cf57c2731a1a2e8d85c80565f1a092def3ac4cf..b2c1d057d435f3fe8efa58cb709e10b49e774afc 100644 (file)
@@ -31,6 +31,11 @@ namespace GridOutFlags
 
   
   
+  Msh::Msh (const bool write_faces) :
+                 write_faces (write_faces)
+  {}
+  
+
   Ucd::Ucd (const bool write_preamble,
            const bool write_faces) :
                  write_preamble (write_preamble),
@@ -138,6 +143,12 @@ void GridOut::set_flags (const GridOutFlags::DX &flags)
 
 
 
+void GridOut::set_flags (const GridOutFlags::Msh &flags) 
+{
+  msh_flags = flags;
+}
+
+
 void GridOut::set_flags (const GridOutFlags::Ucd &flags) 
 {
   ucd_flags = flags;
@@ -197,6 +208,8 @@ GridOut::default_suffix (const OutputFormat output_format)
        return ".eps";
       case xfig:
        return ".fig";
+      case msh:
+       return ".msh";
       default: 
            Assert (false, ExcNotImplemented()); 
            return "";
@@ -226,6 +239,9 @@ GridOut::parse_output_format (const std::string &format_name)
   if (format_name == "xfig")
     return xfig;
   
+  if (format_name == "msh")
+    return msh;
+  
   AssertThrow (false, ExcInvalidState ());
                                   // return something weird
   return OutputFormat(-1);
@@ -235,7 +251,7 @@ GridOut::parse_output_format (const std::string &format_name)
 
 std::string GridOut::get_output_format_names () 
 {
-  return "none|dx|gnuplot|eps|ucd|xfig";
+  return "none|dx|gnuplot|eps|ucd|xfig|msh";
 }
 
 
index e87ec6df0e6500a90710e29cf4f2bf44cdebdd67..7d4952990ec55883e923c79a99603d37a515a037 100644 (file)
@@ -265,7 +265,7 @@ void GridOut::write_dx (const Triangulation<dim> &tria,
 
 
 template <int dim>
-void GridOut::write_ucd (const Triangulation<dim> &tria,
+void GridOut::write_msh (const Triangulation<dim> &tria,
                         std::ostream             &out) 
 {
   AssertThrow (out, ExcIO());
@@ -281,242 +281,380 @@ void GridOut::write_ucd (const Triangulation<dim> &tria,
   typename Triangulation<dim>::active_cell_iterator       cell=tria.begin_active();
   const typename Triangulation<dim>::active_cell_iterator endc=tria.end();
 
-                                  // write preamble
-  if (ucd_flags.write_preamble)
-    {
-                                      // block this to have local
-                                      // variables destroyed after
-                                      // use
-      std::time_t  time1= std::time (0);
-      std::tm     *time = std::localtime(&time1); 
-      out << "# This file was generated by the deal.II library." << '\n'
-         << "# Date =  "
-         << time->tm_year+1900 << "/"
-         << time->tm_mon+1 << "/"
-         << time->tm_mday << '\n'
-         << "# Time =  "
-         << time->tm_hour << ":"
-         << std::setw(2) << time->tm_min << ":"
-         << std::setw(2) << time->tm_sec << '\n'
-         << "#" << '\n'
-         << "# For a description of the UCD format see the AVS Developer's guide."
-         << '\n'
-         << "#" << '\n';
-    };
-
-                                  // start with ucd data
-  out << n_vertices << ' '
-      << tria.n_active_cells() + (ucd_flags.write_faces ?
-                                 n_boundary_faces(tria) :
-                                 0)
-      << " 0 0 0"                  // no data
-      << '\n';
+  // Write Header
+  // The file format is:
+  /* 
 
-                                  // actually write the vertices.
-                                  // note that we shall number them
-                                  // with first index 1 instead of 0
+       
+     $NOD
+     number-of-nodes
+     node-number x-coord y-coord z-coord
+     ...
+     $ENDNOD
+     $ELM
+     number-of-elements
+     elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list
+     ...
+     $ENDELM
+    */
+  out << "$NOD" << std::endl 
+      << n_vertices << std::endl;
+
+  // actually write the vertices.
+  // note that we shall number them
+  // with first index 1 instead of 0
   for (unsigned int i=0; i<vertices.size(); ++i)
-    if (vertex_used[i])
+      if (vertex_used[i])
       {
-       out << i+1                 // vertex index
-           << "  "
-           << vertices[i];
-       for (unsigned int d=dim+1; d<=3; ++d)
-         out << " 0";             // fill with zeroes
-       out << '\n';
+         out << i+1                 // vertex index
+             << "  "
+             << vertices[i];
+         for (unsigned int d=dim+1; d<=3; ++d)
+             out << " 0";             // fill with zeroes
+         out << std::endl;
       };
-       
-                                  // write cells. Enumerate cells
-                                  // consecutively, starting with 1
+
+  // Write cells preamble
+  out << "$ENDNOD" << std::endl
+      << "$ELM" << std::endl
+      << tria.n_active_cells() + (msh_flags.write_faces ?
+               n_boundary_faces(tria) : 0) << std::endl;
+
+  /*
+     elm-type
+     defines the geometrical type of the n-th element:
+     1
+     Line (2 nodes).
+     2
+     Triangle (3 nodes).
+     3
+     Quadrangle (4 nodes).
+     4
+     Tetrahedron (4 nodes).
+     5
+     Hexahedron (8 nodes).
+     6
+     Prism (6 nodes).
+     7
+     Pyramid (5 nodes).
+     8
+     Second order line (3 nodes: 2 associated with the vertices and 1 with the edge).
+     9
+     Second order triangle (6 nodes: 3 associated with the vertices and 3 with the edges).
+     10
+     Second order quadrangle (9 nodes: 4 associated with the vertices, 4 with the edges and 1 with the face).
+     11
+     Second order tetrahedron (10 nodes: 4 associated with the vertices and 6 with the edges).
+     12
+     Second order hexahedron (27 nodes: 8 associated with the vertices, 12 with the edges, 6 with the faces and 1 with the volume).
+     13
+     Second order prism (18 nodes: 6 associated with the vertices, 9 with the edges and 3 with the quadrangular faces).
+     14
+     Second order pyramid (14 nodes: 5 associated with the vertices, 8 with the edges and 1 with the quadrangular face).
+     15
+     Point (1 node).
+     */
+  unsigned int elm_type;
+  switch(dim) {
+      case 1:
+         elm_type = 1;
+         break;
+      case 2:
+         elm_type = 3;
+         break;
+      case 3:
+         elm_type = 5;
+         break;
+      default:
+         Assert(false, ExcNotImplemented());
+  }
+
+  // write cells. Enumerate cells
+  // consecutively, starting with 1
   unsigned int cell_index=1;
   for (cell=tria.begin_active();
-       cell!=endc; ++cell, ++cell_index)
+         cell!=endc; ++cell, ++cell_index)
+  {
+      out << cell_index << ' ' << elm_type << ' '
+         << static_cast<unsigned int>(cell->material_id()) << ' '
+         << cell->subdomain_id() << ' ' 
+         << GeometryInfo<dim>::vertices_per_cell << ' ';
+
+      // Vertex numbering follows UCD conventions.
+
+      for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+             ++vertex)
+         out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex])+1 << ' ';
+      out << std::endl;
+  };
+
+  // write faces with non-zero boundary
+  // indicator
+  if (msh_flags.write_faces)
+      write_msh_faces (tria, cell_index, out);
+
+  out << "$ENDELM" << std::endl;
+
+  // make sure everything now gets to
+  // disk
+  out.flush ();
+
+  AssertThrow (out, ExcIO());
+}
+
+    template <int dim>
+void GridOut::write_ucd (const Triangulation<dim> &tria,
+       std::ostream             &out) 
+{
+    AssertThrow (out, ExcIO());
+
+    // get the positions of the
+    // vertices and whether they are
+    // used.
+    const std::vector<Point<dim> > &vertices    = tria.get_vertices();
+    const std::vector<bool>        &vertex_used = tria.get_used_vertices();
+
+    const unsigned int n_vertices = tria.n_used_vertices();
+
+    typename Triangulation<dim>::active_cell_iterator       cell=tria.begin_active();
+    const typename Triangulation<dim>::active_cell_iterator endc=tria.end();
+
+    // write preamble
+    if (ucd_flags.write_preamble)
     {
-      out << cell_index << ' '
-         << static_cast<unsigned int>(cell->material_id())
-         << " ";
-      switch (dim) 
+       // block this to have local
+       // variables destroyed after
+       // use
+       std::time_t  time1= std::time (0);
+       std::tm     *time = std::localtime(&time1); 
+       out << "# This file was generated by the deal.II library." << '\n'
+           << "# Date =  "
+           << time->tm_year+1900 << "/"
+           << time->tm_mon+1 << "/"
+           << time->tm_mday << '\n'
+           << "# Time =  "
+           << time->tm_hour << ":"
+           << std::setw(2) << time->tm_min << ":"
+           << std::setw(2) << time->tm_sec << '\n'
+           << "#" << '\n'
+           << "# For a description of the UCD format see the AVS Developer's guide."
+           << '\n'
+           << "#" << '\n';
+    };
+
+    // start with ucd data
+    out << n_vertices << ' '
+       << tria.n_active_cells() + (ucd_flags.write_faces ?
+               n_boundary_faces(tria) :
+               0)
+       << " 0 0 0"                  // no data
+       << '\n';
+
+    // actually write the vertices.
+    // note that we shall number them
+    // with first index 1 instead of 0
+    for (unsigned int i=0; i<vertices.size(); ++i)
+       if (vertex_used[i])
        {
-         case 1:  out << "line    "; break;
-         case 2:  out << "quad    "; break;
-         case 3:  out << "hex     "; break;
-         default:
-               Assert (false, ExcNotImplemented());
+           out << i+1                 // vertex index
+               << "  "
+               << vertices[i];
+           for (unsigned int d=dim+1; d<=3; ++d)
+               out << " 0";             // fill with zeroes
+           out << '\n';
        };
 
-                                      // it follows a list of the
-                                      // vertices of each cell. in 1d
-                                      // this is simply a list of the
-                                      // two vertices, in 2d its counter
-                                      // clockwise, as usual in this
-                                      // library. in 3d, the same applies
-                                      // (special thanks to AVS for
-                                      // numbering their vertices in a
-                                      // way compatible to deal.II!)
-                                      //
-                                      // technical reference:
-                                      // AVS Developer's Guide, Release 4,
-                                      // May, 1992, p. E6
-                                      //
-                                      // note: vertex numbers are 1-base      
-      for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
-          ++vertex)
-       out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex])+1 << ' ';
-      out << '\n';
+    // write cells. Enumerate cells
+    // consecutively, starting with 1
+    unsigned int cell_index=1;
+    for (cell=tria.begin_active();
+           cell!=endc; ++cell, ++cell_index)
+    {
+       out << cell_index << ' '
+           << static_cast<unsigned int>(cell->material_id())
+           << " ";
+       switch (dim) 
+       {
+           case 1:  out << "line    "; break;
+           case 2:  out << "quad    "; break;
+           case 3:  out << "hex     "; break;
+           default:
+                    Assert (false, ExcNotImplemented());
+       };
+
+       // it follows a list of the
+       // vertices of each cell. in 1d
+       // this is simply a list of the
+       // two vertices, in 2d its counter
+       // clockwise, as usual in this
+       // library. in 3d, the same applies
+       // (special thanks to AVS for
+       // numbering their vertices in a
+       // way compatible to deal.II!)
+       //
+       // technical reference:
+       // AVS Developer's Guide, Release 4,
+       // May, 1992, p. E6
+       //
+       // note: vertex numbers are 1-base      
+       for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+               ++vertex)
+           out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex])+1 << ' ';
+       out << '\n';
     };
 
-                                  // write faces with non-zero boundary
-                                  // indicator
-  if (ucd_flags.write_faces)
-    write_ucd_faces (tria, cell_index, out);
-    
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
+    // write faces with non-zero boundary
+    // indicator
+    if (ucd_flags.write_faces)
+       write_ucd_faces (tria, cell_index, out);
 
-  AssertThrow (out, ExcIO());
+    // make sure everything now gets to
+    // disk
+    out.flush ();
+
+    AssertThrow (out, ExcIO());
 }
 
 
 #if deal_II_dimension != 2
 
-template <int dim>
+    template <int dim>
 void GridOut::write_xfig (const Triangulation<dim>&,
-                         std::ostream&,
-                         const Mapping<dim>*)
+       std::ostream&,
+       const Mapping<dim>*)
 {
-  Assert (false, ExcNotImplemented());
+    Assert (false, ExcNotImplemented());
 }
 
 #else
 
 //TODO:[GK] Obey parameters
 //TODO:[GK] Flip y-axis?
-template <int dim>
+    template <int dim>
 void GridOut::write_xfig (const Triangulation<dim>& tria,
-                         std::ostream&             out,
-                         const Mapping<dim>*       /*mapping*/)
+       std::ostream&             out,
+       const Mapping<dim>*       /*mapping*/)
 {
-  const unsigned int nv = GeometryInfo<dim>::vertices_per_cell;
-  const unsigned int nf = GeometryInfo<dim>::faces_per_cell;
-  const unsigned int nvf = GeometryInfo<dim>::vertices_per_face;
-
-                                  // The following text was copied
-                                  // from an existing XFig file.
-  out << "#FIG 3.2\nLandscape\nCenter\nInches" << '\n'
-      << "A4\n100.00\nSingle" << '\n'
-                                    // Background is transparent
-      << "-3" << '\n'
-      << "# generated by deal.II GridOut class" << '\n'
-      << "# reduce first number to scale up image" << '\n'
-      << "1200 2" << '\n';
-
-                                  // We write all cells and cells on
-                                  // coarser levels are behind cells
-                                  // on finer levels. Level 0
-                                  // corresponds to a depth of 900,
-                                  // each level subtracting 1
-  typename Triangulation<dim>::cell_iterator cell = tria.begin();
-  const typename Triangulation<dim>::cell_iterator end = tria.end();
-
-  for (;cell != end; ++cell)
+    const unsigned int nv = GeometryInfo<dim>::vertices_per_cell;
+    const unsigned int nf = GeometryInfo<dim>::faces_per_cell;
+    const unsigned int nvf = GeometryInfo<dim>::vertices_per_face;
+
+    // The following text was copied
+    // from an existing XFig file.
+    out << "#FIG 3.2\nLandscape\nCenter\nInches" << '\n'
+       << "A4\n100.00\nSingle" << '\n'
+       // Background is transparent
+       << "-3" << '\n'
+       << "# generated by deal.II GridOut class" << '\n'
+       << "# reduce first number to scale up image" << '\n'
+       << "1200 2" << '\n';
+
+    // We write all cells and cells on
+    // coarser levels are behind cells
+    // on finer levels. Level 0
+    // corresponds to a depth of 900,
+    // each level subtracting 1
+    typename Triangulation<dim>::cell_iterator cell = tria.begin();
+    const typename Triangulation<dim>::cell_iterator end = tria.end();
+
+    for (;cell != end; ++cell)
     {
        // If depth is not encoded, write finest level only
        if (!xfig_flags.level_depth && !cell->active())
            continue;
-                                      // Code for polygon
-      out << "2 3  "
-         << xfig_flags.line_style << ' '
-         << xfig_flags.line_thickness
-                                        // with black line
-         << " 0 ";
-                                        // Fill color
-      if (xfig_flags.level_color)
-       out << cell->level() + 8;
-      else
-       out << cell->material_id() + 1;
-                                      // Depth, unused, fill
-      out << ' '
-         << (xfig_flags.level_depth
-             ? (900-cell->level())
-             : (900+cell->material_id()))
-         << " 0 "
-         << xfig_flags.fill_style << " 0.0 "
-                                        // some style parameters
-         << " 0 0 -1 0 0 "
-                                        // number of points
-         << nv+1 << '\n';
-
-                                      // For each point, write scaled
-                                      // and shifted coordinates
-                                      // multiplied by 1200
-                                      // (dots/inch)
-      for (unsigned int k=0;k<=nv;++k)
+       // Code for polygon
+       out << "2 3  "
+           << xfig_flags.line_style << ' '
+           << xfig_flags.line_thickness
+           // with black line
+           << " 0 ";
+       // Fill color
+       if (xfig_flags.level_color)
+           out << cell->level() + 8;
+       else
+           out << cell->material_id() + 1;
+       // Depth, unused, fill
+       out << ' '
+           << (xfig_flags.level_depth
+                   ? (900-cell->level())
+                   : (900+cell->material_id()))
+           << " 0 "
+           << xfig_flags.fill_style << " 0.0 "
+           // some style parameters
+           << " 0 0 -1 0 0 "
+           // number of points
+           << nv+1 << '\n';
+
+       // For each point, write scaled
+       // and shifted coordinates
+       // multiplied by 1200
+       // (dots/inch)
+       for (unsigned int k=0;k<=nv;++k)
        {
-         const Point<dim>& p = cell->vertex(
-           GeometryInfo<dim>::ucd_to_deal[k % nv]);
-         for (unsigned int d=0;d<dim;++d)
+           const Point<dim>& p = cell->vertex(
+                   GeometryInfo<dim>::ucd_to_deal[k % nv]);
+           for (unsigned int d=0;d<dim;++d)
            {
-             int val = (int)(1200 * xfig_flags.scaling(d) *
-                             (p(d)-xfig_flags.offset(d)));
-             out << '\t' << val;
+               int val = (int)(1200 * xfig_flags.scaling(d) *
+                       (p(d)-xfig_flags.offset(d)));
+               out << '\t' << val;
            }
-         out << '\n';
+           out << '\n';
        }
-                                      // Now write boundary edges
-      static const unsigned int face_reorder[4]={2,1,3,0};
-      if (xfig_flags.draw_boundary)
-       for (unsigned int f=0;f<nf;++f)
-         {
-           typename Triangulation<dim>::face_iterator
-             face = cell->face(face_reorder[f]);
-           const unsigned char bi = face->boundary_indicator();
-           if (bi != 255)
-             {
-                                                // Code for polyline
-               out << "2 1 "
-                                                  // with line style and thickness
-                   << xfig_flags.boundary_style << ' '
-                   << xfig_flags.boundary_thickness << ' '
-                   << (1 + (unsigned int) bi);
-                                                // Fill color
-               out << " -1 ";
-                                                // Depth 100 less than cells
-               out << (xfig_flags.level_depth
-                       ? (800-cell->level())
-                       : 800+bi)
-                                                  // unused, no fill
-                   << " 0 -1 0.0 "
-                                                  // some style parameters
-                   << " 0 0 -1 0 0 "
-                                                  // number of points
-                   << nvf << '\n';
-               
-                                                // For each point, write scaled
-                                                // and shifted coordinates
-                                                // multiplied by 1200
-                                                // (dots/inch)
-               
-               for (unsigned int k=0;k<nvf;++k)
-                 {
-                   const Point<dim>& p = face->vertex(k % nv);
-                   for (unsigned int d=0;d<dim;++d)
-                     {
-                       int val = (int)(1200 * xfig_flags.scaling(d) *
-                                       (p(d)-xfig_flags.offset(d)));
-                       out << '\t' << val;
-                     }
-                   out << '\n';
-                 } 
-             }
-         }
+       // Now write boundary edges
+       static const unsigned int face_reorder[4]={2,1,3,0};
+       if (xfig_flags.draw_boundary)
+           for (unsigned int f=0;f<nf;++f)
+           {
+               typename Triangulation<dim>::face_iterator
+                   face = cell->face(face_reorder[f]);
+               const unsigned char bi = face->boundary_indicator();
+               if (bi != 255)
+               {
+                   // Code for polyline
+                   out << "2 1 "
+                       // with line style and thickness
+                       << xfig_flags.boundary_style << ' '
+                       << xfig_flags.boundary_thickness << ' '
+                       << (1 + (unsigned int) bi);
+                   // Fill color
+                   out << " -1 ";
+                   // Depth 100 less than cells
+                   out << (xfig_flags.level_depth
+                           ? (800-cell->level())
+                           : 800+bi)
+                       // unused, no fill
+                       << " 0 -1 0.0 "
+                       // some style parameters
+                       << " 0 0 -1 0 0 "
+                       // number of points
+                       << nvf << '\n';
+
+                   // For each point, write scaled
+                   // and shifted coordinates
+                   // multiplied by 1200
+                   // (dots/inch)
+
+                   for (unsigned int k=0;k<nvf;++k)
+                   {
+                       const Point<dim>& p = face->vertex(k % nv);
+                       for (unsigned int d=0;d<dim;++d)
+                       {
+                           int val = (int)(1200 * xfig_flags.scaling(d) *
+                                   (p(d)-xfig_flags.offset(d)));
+                           out << '\t' << val;
+                       }
+                       out << '\n';
+                   
+               }
+           }
     }
 
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
+    // make sure everything now gets to
+    // disk
+    out.flush ();
 
-  AssertThrow (out, ExcIO());
+    AssertThrow (out, ExcIO());
 }
 
 #endif
@@ -527,7 +665,7 @@ void GridOut::write_xfig (const Triangulation<dim>& tria,
 
 unsigned int GridOut::n_boundary_faces (const Triangulation<1> &) const
 {
-  return 0;
+    return 0;
 }
 
 #endif
@@ -537,16 +675,63 @@ unsigned int GridOut::n_boundary_faces (const Triangulation<1> &) const
 template <int dim>
 unsigned int GridOut::n_boundary_faces (const Triangulation<dim> &tria) const
 {
-  typename Triangulation<dim>::active_face_iterator face, endf;
-  unsigned int n_faces = 0;
+    typename Triangulation<dim>::active_face_iterator face, endf;
+    unsigned int n_faces = 0;
 
-  for (face=tria.begin_active_face(), endf=tria.end_face();
-       face != endf; ++face)
-    if ((face->at_boundary()) &&
-       (face->boundary_indicator() != 0))
-      n_faces++;
+    for (face=tria.begin_active_face(), endf=tria.end_face();
+           face != endf; ++face)
+       if ((face->at_boundary()) &&
+               (face->boundary_indicator() != 0))
+           n_faces++;
 
-  return n_faces;
+    return n_faces;
+}
+
+
+
+#if deal_II_dimension == 1
+
+void GridOut::write_msh_faces (const Triangulation<1> &,
+       const unsigned int,
+       std::ostream &) const
+{
+    return;
+}
+
+#endif
+
+
+
+template <int dim>
+void GridOut::write_msh_faces (const Triangulation<dim> &tria,
+       const unsigned int        starting_index,
+       std::ostream             &out) const
+{
+    typename Triangulation<dim>::active_face_iterator face, endf;
+    unsigned int index=starting_index;
+
+    for (face=tria.begin_active_face(), endf=tria.end_face();
+           face != endf; ++face)
+       if (face->at_boundary() &&
+               (face->boundary_indicator() != 0)) 
+       {
+           out << index << " ";
+           switch (dim) 
+           {
+               case 2: out << 1 << ' ';  break;
+               case 3: out << 3 << ' ';  break;
+               default:
+                       Assert (false, ExcNotImplemented());
+           };
+           out << static_cast<unsigned int>(face->boundary_indicator()) 
+               << " 1 " << GeometryInfo<dim>::vertices_per_face;
+           // note: vertex numbers are 1-base
+           for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
+               out << ' ' << face->vertex_index(GeometryInfo<dim-1>::ucd_to_deal[vertex])+1;
+           out << '\n';
+
+           ++index;
+       };        
 }
 
 
@@ -554,10 +739,10 @@ unsigned int GridOut::n_boundary_faces (const Triangulation<dim> &tria) const
 #if deal_II_dimension == 1
 
 void GridOut::write_ucd_faces (const Triangulation<1> &,
-                              const unsigned int,
-                              std::ostream &) const
+       const unsigned int,
+       std::ostream &) const
 {
-  return;
+    return;
 }
 
 #endif
@@ -566,34 +751,34 @@ void GridOut::write_ucd_faces (const Triangulation<1> &,
 
 template <int dim>
 void GridOut::write_ucd_faces (const Triangulation<dim> &tria,
-                              const unsigned int        starting_index,
-                              std::ostream             &out) const
+       const unsigned int        starting_index,
+       std::ostream             &out) const
 {
-  typename Triangulation<dim>::active_face_iterator face, endf;
-  unsigned int index=starting_index;
+    typename Triangulation<dim>::active_face_iterator face, endf;
+    unsigned int index=starting_index;
 
-  for (face=tria.begin_active_face(), endf=tria.end_face();
-       face != endf; ++face)
-    if (face->at_boundary() &&
-       (face->boundary_indicator() != 0)) 
-      {
-       out << index << "  "
-           << static_cast<unsigned int>(face->boundary_indicator())
-           << "  ";
-       switch (dim) 
-         {
-           case 2: out << "line    ";  break;
-           case 3: out << "quad    ";  break;
-           default:
-                 Assert (false, ExcNotImplemented());
-         };
-                                      // note: vertex numbers are 1-base
-       for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
-         out << face->vertex_index(GeometryInfo<dim-1>::ucd_to_deal[vertex])+1 << ' ';
-       out << '\n';
+    for (face=tria.begin_active_face(), endf=tria.end_face();
+           face != endf; ++face)
+       if (face->at_boundary() &&
+               (face->boundary_indicator() != 0)) 
+       {
+           out << index << "  "
+               << static_cast<unsigned int>(face->boundary_indicator())
+               << "  ";
+           switch (dim) 
+           {
+               case 2: out << "line    ";  break;
+               case 3: out << "quad    ";  break;
+               default:
+                       Assert (false, ExcNotImplemented());
+           };
+           // note: vertex numbers are 1-base
+           for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
+               out << face->vertex_index(GeometryInfo<dim-1>::ucd_to_deal[vertex])+1 << ' ';
+           out << '\n';
 
-       ++index;
-      };         
+           ++index;
+       };        
 }
 
 
@@ -601,33 +786,33 @@ void GridOut::write_ucd_faces (const Triangulation<dim> &tria,
 #if deal_II_dimension==1
 
 void GridOut::write_gnuplot (const Triangulation<1> &tria,
-                            std::ostream           &out,
-                            const Mapping<1>       *) 
+       std::ostream           &out,
+       const Mapping<1>       *) 
 {
-  AssertThrow (out, ExcIO());
+    AssertThrow (out, ExcIO());
 
-  Triangulation<1>::active_cell_iterator        cell=tria.begin_active();
-  const Triangulation<1>::active_cell_iterator  endc=tria.end();
-  for (; cell!=endc; ++cell)
+    Triangulation<1>::active_cell_iterator        cell=tria.begin_active();
+    const Triangulation<1>::active_cell_iterator  endc=tria.end();
+    for (; cell!=endc; ++cell)
     {
-      if (gnuplot_flags.write_cell_numbers)
-       out << "# cell " << cell << '\n';
-
-      out << cell->vertex(0)
-         << ' ' << cell->level()
-         << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-         << cell->vertex(1)
-                   << ' ' << cell->level()
-         << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-         << '\n';
-      break;
+       if (gnuplot_flags.write_cell_numbers)
+           out << "# cell " << cell << '\n';
+
+       out << cell->vertex(0)
+           << ' ' << cell->level()
+           << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+           << cell->vertex(1)
+           << ' ' << cell->level()
+           << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+           << '\n';
+       break;
     }
 
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
+    // make sure everything now gets to
+    // disk
+    out.flush ();
 
-  AssertThrow (out, ExcIO());
+    AssertThrow (out, ExcIO());
 }
 
 #endif
@@ -636,124 +821,124 @@ void GridOut::write_gnuplot (const Triangulation<1> &tria,
 #if deal_II_dimension==2
 
 void GridOut::write_gnuplot (const Triangulation<2> &tria,
-                            std::ostream           &out,
-                            const Mapping<2>       *mapping) 
+       std::ostream           &out,
+       const Mapping<2>       *mapping) 
 {
-  AssertThrow (out, ExcIO());
-
-  const unsigned int dim=2;
-  const unsigned int n_additional_points=
-    gnuplot_flags.n_boundary_face_points;
-  const unsigned int n_points=2+n_additional_points;
-  
-  Triangulation<dim>::active_cell_iterator        cell=tria.begin_active();
-  const Triangulation<dim>::active_cell_iterator  endc=tria.end();
-
-                                  // if we are to treat curved
-                                  // boundaries, then generate a
-                                  // quadrature formula which will be
-                                  // used to probe boundary points at
-                                  // curved faces
-  Quadrature<dim> *q_projector=0;
-  std::vector<Point<1> > boundary_points;
-  if (mapping!=0)
+    AssertThrow (out, ExcIO());
+
+    const unsigned int dim=2;
+    const unsigned int n_additional_points=
+       gnuplot_flags.n_boundary_face_points;
+    const unsigned int n_points=2+n_additional_points;
+
+    Triangulation<dim>::active_cell_iterator        cell=tria.begin_active();
+    const Triangulation<dim>::active_cell_iterator  endc=tria.end();
+
+    // if we are to treat curved
+    // boundaries, then generate a
+    // quadrature formula which will be
+    // used to probe boundary points at
+    // curved faces
+    Quadrature<dim> *q_projector=0;
+    std::vector<Point<1> > boundary_points;
+    if (mapping!=0)
     {
-      boundary_points.resize(n_points);
-      boundary_points[0][0]=0;
-      boundary_points[n_points-1][0]=1;      
-      for (unsigned int i=1; i<n_points-1; ++i)
-        boundary_points[i](0)= 1.*i/(n_points-1);
+       boundary_points.resize(n_points);
+       boundary_points[0][0]=0;
+       boundary_points[n_points-1][0]=1;      
+       for (unsigned int i=1; i<n_points-1; ++i)
+           boundary_points[i](0)= 1.*i/(n_points-1);
 
-      std::vector<double> dummy_weights(n_points, 1./n_points);
-      Quadrature<1> quadrature(boundary_points, dummy_weights);
+       std::vector<double> dummy_weights(n_points, 1./n_points);
+       Quadrature<1> quadrature(boundary_points, dummy_weights);
 
-      q_projector = new Quadrature<dim> (QProjector<dim>::project_to_all_faces(quadrature));
+       q_projector = new Quadrature<dim> (QProjector<dim>::project_to_all_faces(quadrature));
     }
-  
-  for (; cell!=endc; ++cell)
+
+    for (; cell!=endc; ++cell)
     {
-      if (gnuplot_flags.write_cell_numbers)
-       out << "# cell " << cell << '\n';
+       if (gnuplot_flags.write_cell_numbers)
+           out << "# cell " << cell << '\n';
 
-      if (mapping==0 || !cell->at_boundary())
+       if (mapping==0 || !cell->at_boundary())
        {
-                                          // write out the four sides
-                                          // of this cell by putting
-                                          // the four points (+ the
-                                          // initial point again) in
-                                          // a row and lifting the
-                                          // drawing pencil at the
-                                          // end
-         for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-           out << cell->vertex(GeometryInfo<dim>::ucd_to_deal[i])
+           // write out the four sides
+           // of this cell by putting
+           // the four points (+ the
+           // initial point again) in
+           // a row and lifting the
+           // drawing pencil at the
+           // end
+           for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+               out << cell->vertex(GeometryInfo<dim>::ucd_to_deal[i])
+                   << ' ' << cell->level()
+                   << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+           out << cell->vertex(0)
                << ' ' << cell->level()
-               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
-         out << cell->vertex(0)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << '\n'  // double new line for gnuplot 3d plots
-             << '\n';
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << '\n'  // double new line for gnuplot 3d plots
+               << '\n';
        }
-      else
-                                        // cell is at boundary and we
-                                        // are to treat curved
-                                        // boundaries. so loop over
-                                        // all faces and draw them as
-                                        // small pieces of lines
+       else
+           // cell is at boundary and we
+           // are to treat curved
+           // boundaries. so loop over
+           // all faces and draw them as
+           // small pieces of lines
        {
-         for (unsigned int face_no=0;
-              face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+           for (unsigned int face_no=0;
+                   face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
            {
-             const Triangulation<dim>::face_iterator
-               face = cell->face(face_no);
-             if (face->at_boundary())
+               const Triangulation<dim>::face_iterator
+                   face = cell->face(face_no);
+               if (face->at_boundary())
                {
-                                                  // compute offset
-                                                  // of quadrature
-                                                  // points within
-                                                  // set of projected
-                                                  // points
-                 const unsigned int offset=face_no*n_points;
-                 for (unsigned int i=0; i<n_points; ++i)
-                   out << (mapping->transform_unit_to_real_cell
-                           (cell, q_projector->point(offset+i)))
-                       << ' ' << cell->level()
-                       << ' ' << static_cast<unsigned int>(cell->material_id())
+                   // compute offset
+                   // of quadrature
+                   // points within
+                   // set of projected
+                   // points
+                   const unsigned int offset=face_no*n_points;
+                   for (unsigned int i=0; i<n_points; ++i)
+                       out << (mapping->transform_unit_to_real_cell
+                               (cell, q_projector->point(offset+i)))
+                           << ' ' << cell->level()
+                           << ' ' << static_cast<unsigned int>(cell->material_id())
+                           << '\n';
+
+                   out << '\n'
                        << '\n';
-                 
-                 out << '\n'
-                     << '\n';
                }
-             else
+               else
                {
-                                                  // if, however, the
-                                                  // face is not at
-                                                  // the boundary,
-                                                  // then draw it as
-                                                  // usual
-                 out << face->vertex(0)
-                     << ' ' << cell->level()
-                     << ' ' << static_cast<unsigned int>(cell->material_id())
-                     << '\n'
-                     << face->vertex(1)
-                     << ' ' << cell->level()
-                     << ' ' << static_cast<unsigned int>(cell->material_id())
-                     << '\n'
-                     << '\n'
-                     << '\n';
+                   // if, however, the
+                   // face is not at
+                   // the boundary,
+                   // then draw it as
+                   // usual
+                   out << face->vertex(0)
+                       << ' ' << cell->level()
+                       << ' ' << static_cast<unsigned int>(cell->material_id())
+                       << '\n'
+                       << face->vertex(1)
+                       << ' ' << cell->level()
+                       << ' ' << static_cast<unsigned int>(cell->material_id())
+                       << '\n'
+                       << '\n'
+                       << '\n';
                }
            }
        }
     }
 
-  if (q_projector != 0)
-    delete q_projector;  
-  
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
+    if (q_projector != 0)
+       delete q_projector;  
 
-  AssertThrow (out, ExcIO());
+    // make sure everything now gets to
+    // disk
+    out.flush ();
+
+    AssertThrow (out, ExcIO());
 }
 
 #endif
@@ -762,207 +947,207 @@ void GridOut::write_gnuplot (const Triangulation<2> &tria,
 #if deal_II_dimension==3
 
 void GridOut::write_gnuplot (const Triangulation<3> &tria,
-                            std::ostream           &out,
-                            const Mapping<3>       *mapping) 
+       std::ostream           &out,
+       const Mapping<3>       *mapping) 
 {
-  AssertThrow (out, ExcIO());
-
-  const unsigned int dim=3;
-  const unsigned int n_additional_points=
-    gnuplot_flags.n_boundary_face_points;
-  const unsigned int n_points=2+n_additional_points;
-  
-  Triangulation<dim>::active_cell_iterator        cell=tria.begin_active();
-  const Triangulation<dim>::active_cell_iterator  endc=tria.end();
-
-                                  // if we are to treat curved
-                                  // boundaries, then generate a
-                                  // quadrature formula which will be
-                                  // used to probe boundary points at
-                                  // curved faces
-  Quadrature<dim> *q_projector=0;
-  std::vector<Point<1> > boundary_points;
-  if (mapping!=0)
+    AssertThrow (out, ExcIO());
+
+    const unsigned int dim=3;
+    const unsigned int n_additional_points=
+       gnuplot_flags.n_boundary_face_points;
+    const unsigned int n_points=2+n_additional_points;
+
+    Triangulation<dim>::active_cell_iterator        cell=tria.begin_active();
+    const Triangulation<dim>::active_cell_iterator  endc=tria.end();
+
+    // if we are to treat curved
+    // boundaries, then generate a
+    // quadrature formula which will be
+    // used to probe boundary points at
+    // curved faces
+    Quadrature<dim> *q_projector=0;
+    std::vector<Point<1> > boundary_points;
+    if (mapping!=0)
     {
-      boundary_points.resize(n_points);
-      boundary_points[0][0]=0;
-      boundary_points[n_points-1][0]=1;      
-      for (unsigned int i=1; i<n_points-1; ++i)
-        boundary_points[i](0)= 1.*i/(n_points-1);
-
-      std::vector<double> dummy_weights(n_points, 1./n_points);
-      Quadrature<1> quadrature1d(boundary_points, dummy_weights);
-
-                                      // tensor product of points,
-                                      // only one copy
-      QIterated<dim-1> quadrature(quadrature1d, 1);
-      q_projector = new Quadrature<dim> (QProjector<dim>::project_to_all_faces(quadrature));
+       boundary_points.resize(n_points);
+       boundary_points[0][0]=0;
+       boundary_points[n_points-1][0]=1;      
+       for (unsigned int i=1; i<n_points-1; ++i)
+           boundary_points[i](0)= 1.*i/(n_points-1);
+
+       std::vector<double> dummy_weights(n_points, 1./n_points);
+       Quadrature<1> quadrature1d(boundary_points, dummy_weights);
+
+       // tensor product of points,
+       // only one copy
+       QIterated<dim-1> quadrature(quadrature1d, 1);
+       q_projector = new Quadrature<dim> (QProjector<dim>::project_to_all_faces(quadrature));
     }
-  
-  for (; cell!=endc; ++cell)
+
+    for (; cell!=endc; ++cell)
     {
-      if (gnuplot_flags.write_cell_numbers)
-       out << "# cell " << cell << '\n';
+       if (gnuplot_flags.write_cell_numbers)
+           out << "# cell " << cell << '\n';
 
-      if (mapping==0 || n_points==2 || !cell->has_boundary_lines())
+       if (mapping==0 || n_points==2 || !cell->has_boundary_lines())
        {
-                                          // front face
-         out << cell->vertex(0)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(1)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(5)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(4)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(0)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << '\n';
-                                          // back face
-         out << cell->vertex(2)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(3)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(7)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(6)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(2)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << '\n';
-         
-                                          // now for the four connecting lines
-         out << cell->vertex(0)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(2)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << '\n';
-         out << cell->vertex(1)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(3)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << '\n';
-         out << cell->vertex(5)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(7)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << '\n';
-         out << cell->vertex(4)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << cell->vertex(6)
-             << ' ' << cell->level()
-             << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-             << '\n';
+           // front face
+           out << cell->vertex(0)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(1)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(5)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(4)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(0)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << '\n';
+           // back face
+           out << cell->vertex(2)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(3)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(7)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(6)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(2)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << '\n';
+
+           // now for the four connecting lines
+           out << cell->vertex(0)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(2)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << '\n';
+           out << cell->vertex(1)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(3)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << '\n';
+           out << cell->vertex(5)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(7)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << '\n';
+           out << cell->vertex(4)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << cell->vertex(6)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+               << '\n';
        }
-      else
+       else
        {
-         for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+           for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
            {
-             const Triangulation<dim>::face_iterator
-               face = cell->face(face_no);
-             
-             if (face->at_boundary())
+               const Triangulation<dim>::face_iterator
+                   face = cell->face(face_no);
+
+               if (face->at_boundary())
                {
-                 const unsigned int offset=face_no*n_points*n_points;
-                 for (unsigned int i=0; i<n_points-1; ++i)
-                   for (unsigned int j=0; j<n_points-1; ++j)
-                     {
-                       const Point<dim> p0=mapping->transform_unit_to_real_cell(
-                         cell, q_projector->point(offset+i*n_points+j));
-                       out << p0
-                           << ' ' << cell->level()
-                           << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
-                       out << (mapping->transform_unit_to_real_cell(
-                                 cell, q_projector->point(offset+(i+1)*n_points+j)))
-                           << ' ' << cell->level()
-                           << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
-                       out << (mapping->transform_unit_to_real_cell(
-                                 cell, q_projector->point(offset+(i+1)*n_points+j+1)))
-                           << ' ' << cell->level()
-                           << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
-                       out << (mapping->transform_unit_to_real_cell(
-                                 cell, q_projector->point(offset+i*n_points+j+1)))
-                           << ' ' << cell->level()
-                           << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
-                                                        // and the
-                                                        // first
-                                                        // point
-                                                        // again
-                       out << p0
-                           << ' ' << cell->level()
-                           << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
-                       out << '\n' << '\n';
-                     }
+                   const unsigned int offset=face_no*n_points*n_points;
+                   for (unsigned int i=0; i<n_points-1; ++i)
+                       for (unsigned int j=0; j<n_points-1; ++j)
+                       {
+                           const Point<dim> p0=mapping->transform_unit_to_real_cell(
+                                   cell, q_projector->point(offset+i*n_points+j));
+                           out << p0
+                               << ' ' << cell->level()
+                               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+                           out << (mapping->transform_unit_to_real_cell(
+                                       cell, q_projector->point(offset+(i+1)*n_points+j)))
+                               << ' ' << cell->level()
+                               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+                           out << (mapping->transform_unit_to_real_cell(
+                                       cell, q_projector->point(offset+(i+1)*n_points+j+1)))
+                               << ' ' << cell->level()
+                               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+                           out << (mapping->transform_unit_to_real_cell(
+                                       cell, q_projector->point(offset+i*n_points+j+1)))
+                               << ' ' << cell->level()
+                               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+                           // and the
+                           // first
+                           // point
+                           // again
+                           out << p0
+                               << ' ' << cell->level()
+                               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+                           out << '\n' << '\n';
+                       }
                }
-             else
+               else
                {
-                 for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
+                   for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
                    {
-                     const Triangulation<dim>::line_iterator
-                       line=face->line(l);
-                     
-                     const Point<dim> &v0=line->vertex(0),
-                                      &v1=line->vertex(1);
-                     if (line->at_boundary())
+                       const Triangulation<dim>::line_iterator
+                           line=face->line(l);
+
+                       const Point<dim> &v0=line->vertex(0),
+                             &v1=line->vertex(1);
+                       if (line->at_boundary())
                        {
-                                                          // transform_real_to_unit_cell
-                                                          // could be
-                                                          // replaced
-                                                          // by using
-                                                          // QProjector<3>::project_to_line
-                                                          // which is
-                                                          // not yet
-                                                          // implemented
-                         const Point<dim> u0=mapping->transform_real_to_unit_cell(cell, v0),
-                                          u1=mapping->transform_real_to_unit_cell(cell, v1);
-                         
-                         for (unsigned int i=0; i<n_points; ++i)
-                           out << (mapping->transform_unit_to_real_cell
-                                   (cell, (1-boundary_points[i][0])*u0+boundary_points[i][0]*u1))
+                           // transform_real_to_unit_cell
+                           // could be
+                           // replaced
+                           // by using
+                           // QProjector<3>::project_to_line
+                           // which is
+                           // not yet
+                           // implemented
+                           const Point<dim> u0=mapping->transform_real_to_unit_cell(cell, v0),
+                                 u1=mapping->transform_real_to_unit_cell(cell, v1);
+
+                           for (unsigned int i=0; i<n_points; ++i)
+                               out << (mapping->transform_unit_to_real_cell
+                                       (cell, (1-boundary_points[i][0])*u0+boundary_points[i][0]*u1))
+                                   << ' ' << cell->level()
+                                   << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+                       }
+                       else
+                           out << v0
+                               << ' ' << cell->level()
+                               << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+                               << v1
                                << ' ' << cell->level()
                                << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
-                       }
-                     else
-                       out << v0
-                           << ' ' << cell->level()
-                           << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
-                           << v1
-                           << ' ' << cell->level()
-                           << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
 
-                     out << '\n' << '\n';
+                       out << '\n' << '\n';
                    }
                }
            }
        }
     }
 
-  if (q_projector != 0)
-    delete q_projector;
-  
-  
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
+    if (q_projector != 0)
+       delete q_projector;
 
-  AssertThrow (out, ExcIO());
+
+    // make sure everything now gets to
+    // disk
+    out.flush ();
+
+    AssertThrow (out, ExcIO());
 }
 
 #endif
@@ -974,493 +1159,497 @@ struct LineEntry
     bool colorize;
     unsigned int level;
     LineEntry (const Point<2>    &f,
-              const Point<2>    &s,
-              const bool         c,
-              const unsigned int l)
-                   :
-                   first(f), second(s),
-                   colorize(c), level(l)
-      {}
+           const Point<2>    &s,
+           const bool         c,
+           const unsigned int l)
+       :
+       first(f), second(s),
+       colorize(c), level(l)
+       {}
 };
 
 
 #if deal_II_dimension==1
 
 void GridOut::write_eps (const Triangulation<1> &,
-                        std::ostream &,
-                        const Mapping<1> *) 
+       std::ostream &,
+       const Mapping<1> *) 
 {
-  Assert(false, ExcNotImplemented());
+    Assert(false, ExcNotImplemented());
 }
 
 
 #else
 
-template <int dim>
+    template <int dim>
 void GridOut::write_eps (const Triangulation<dim> &tria,
-                        std::ostream             &out,
-                        const Mapping<dim>       *mapping) 
+       std::ostream             &out,
+       const Mapping<dim>       *mapping) 
 {
-  
-  typedef std::list<LineEntry> LineList;
-
-                                  // get a pointer to the flags
-                                  // common to all dimensions,
-                                  // in order to avoid the recurring
-                                  // distinctions between
-                                  // eps_flags_1, eps_flags_2, ...
-  const GridOutFlags::EpsFlagsBase
-    &eps_flags_base = (dim==2 ?
-                      static_cast<GridOutFlags::EpsFlagsBase&>(eps_flags_2) :
-                      (dim==3 ?
-                       static_cast<GridOutFlags::EpsFlagsBase&>(eps_flags_3) :
-                       *static_cast<GridOutFlags::EpsFlagsBase*>(0)));
-  
-  AssertThrow (out, ExcIO());
-  const unsigned int n_points = eps_flags_base.n_boundary_face_points;
-
-                                  // make up a list of lines by which
-                                  // we will construct the triangulation
-                                  //
-                                  // this part unfortunately is a bit
-                                  // dimension dependent, so we have to
-                                  // treat every dimension different.
-                                  // however, by directly producing
-                                  // the lines to be printed, i.e. their
-                                  // 2d images, we can later do the
-                                  // actual output dimension independent
-                                  // again
-  LineList line_list;
-
-  switch (dim)
+
+    typedef std::list<LineEntry> LineList;
+
+    // get a pointer to the flags
+    // common to all dimensions,
+    // in order to avoid the recurring
+    // distinctions between
+    // eps_flags_1, eps_flags_2, ...
+    const GridOutFlags::EpsFlagsBase
+       &eps_flags_base = (dim==2 ?
+               static_cast<GridOutFlags::EpsFlagsBase&>(eps_flags_2) :
+               (dim==3 ?
+                static_cast<GridOutFlags::EpsFlagsBase&>(eps_flags_3) :
+                *static_cast<GridOutFlags::EpsFlagsBase*>(0)));
+
+    AssertThrow (out, ExcIO());
+    const unsigned int n_points = eps_flags_base.n_boundary_face_points;
+
+    // make up a list of lines by which
+    // we will construct the triangulation
+    //
+    // this part unfortunately is a bit
+    // dimension dependent, so we have to
+    // treat every dimension different.
+    // however, by directly producing
+    // the lines to be printed, i.e. their
+    // 2d images, we can later do the
+    // actual output dimension independent
+    // again
+    LineList line_list;
+
+    switch (dim)
     {
-      case 1:
-      {
-       Assert(false, ExcInternalError());
-       break;
-      };
-       
-      case 2:
-      {
-       typename Triangulation<dim>::active_line_iterator
-         line   =tria.begin_active_line (),
-         endline=tria.end_line ();
-
-                                        // first treat all interior
-                                        // lines and make up a list
-                                        // of them. if curved lines
-                                        // shall not be supported
-                                        // (i.e. no mapping is
-                                        // provided), then also treat
-                                        // all other lines
-       for (; line!=endline; ++line)
-         if (mapping==0 || !line->at_boundary())
-                                            // one would expect
-                                            // make_pair(line->vertex(0),
-                                            //           line->vertex(1))
-                                            // here, but that is not
-                                            // dimension independent, since
-                                            // vertex(i) is Point<dim>,
-                                            // but we want a Point<2>.
-                                            // in fact, whenever we're here,
-                                            // the vertex is a Point<dim>,
-                                            // but the compiler does not
-                                            // know this. hopefully, the
-                                            // compiler will optimize away
-                                            // this little kludge
-           line_list.push_back (LineEntry(Point<2>(line->vertex(0)(0),
-                                                   line->vertex(0)(1)),
-                                          Point<2>(line->vertex(1)(0),
-                                                   line->vertex(1)(1)),
-                                          line->user_flag_set(),
-                                          line->level()));
-       
-                                        // next if we are to treat
-                                        // curved boundaries
-                                        // specially, then add lines
-                                        // to the list consisting of
-                                        // pieces of the boundary
-                                        // lines
-       if (mapping!=0)
-         {
-                                            // to do so, first
-                                            // generate a sequence of
-                                            // points on a face and
-                                            // project them onto the
-                                            // faces of a unit cell
-           std::vector<Point<dim-1> > boundary_points (n_points);
-           
-           for (unsigned int i=0; i<n_points; ++i)
-             boundary_points[i](0) = 1.*(i+1)/(n_points+1);
-           
-           Quadrature<dim-1> quadrature (boundary_points);
-           Quadrature<dim>   q_projector (QProjector<dim>::project_to_all_faces(quadrature));
-
-                                            // next loop over all
-                                            // boundary faces and
-                                            // generate the info from
-                                            // them
-           typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active ();
-           const typename Triangulation<dim>::active_cell_iterator end=tria.end ();
-           for (; cell!=end; ++cell)
-             for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+       case 1:
+           {
+               Assert(false, ExcInternalError());
+               break;
+           };
+
+       case 2:
+           {
+               typename Triangulation<dim>::active_line_iterator
+                   line   =tria.begin_active_line (),
+                          endline=tria.end_line ();
+
+               // first treat all interior
+               // lines and make up a list
+               // of them. if curved lines
+               // shall not be supported
+               // (i.e. no mapping is
+               // provided), then also treat
+               // all other lines
+               for (; line!=endline; ++line)
+                   if (mapping==0 || !line->at_boundary())
+                       // one would expect
+                       // make_pair(line->vertex(0),
+                       //           line->vertex(1))
+                       // here, but that is not
+                       // dimension independent, since
+                       // vertex(i) is Point<dim>,
+                       // but we want a Point<2>.
+                       // in fact, whenever we're here,
+                       // the vertex is a Point<dim>,
+                       // but the compiler does not
+                       // know this. hopefully, the
+                       // compiler will optimize away
+                       // this little kludge
+                       line_list.push_back (LineEntry(Point<2>(line->vertex(0)(0),
+                                       line->vertex(0)(1)),
+                                   Point<2>(line->vertex(1)(0),
+                                       line->vertex(1)(1)),
+                                   line->user_flag_set(),
+                                   line->level()));
+
+               // next if we are to treat
+               // curved boundaries
+               // specially, then add lines
+               // to the list consisting of
+               // pieces of the boundary
+               // lines
+               if (mapping!=0)
                {
-                 const typename Triangulation<dim>::face_iterator
-                   face = cell->face(face_no);
-                 
-                 if (face->at_boundary())
-                   {
-                     Point<dim> p0_dim(face->vertex(0));
-                     Point<2>   p0    (p0_dim(0), p0_dim(1));
-
-                                                      // loop over
-                                                      // all pieces
-                                                      // of the line
-                                                      // and generate
-                                                      // line-lets
-                     const unsigned int offset=face_no*n_points;
-                     for (unsigned int i=0; i<n_points; ++i)
+                   // to do so, first
+                   // generate a sequence of
+                   // points on a face and
+                   // project them onto the
+                   // faces of a unit cell
+                   std::vector<Point<dim-1> > boundary_points (n_points);
+
+                   for (unsigned int i=0; i<n_points; ++i)
+                       boundary_points[i](0) = 1.*(i+1)/(n_points+1);
+
+                   Quadrature<dim-1> quadrature (boundary_points);
+                   Quadrature<dim>   q_projector (QProjector<dim>::project_to_all_faces(quadrature));
+
+                   // next loop over all
+                   // boundary faces and
+                   // generate the info from
+                   // them
+                   typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active ();
+                   const typename Triangulation<dim>::active_cell_iterator end=tria.end ();
+                   for (; cell!=end; ++cell)
+                       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
                        {
-                         const Point<dim> p1_dim (mapping->transform_unit_to_real_cell
-                                                  (cell, q_projector.point(offset+i)));
-                         const Point<2>   p1     (p1_dim(0), p1_dim(1));
-                         
-                         line_list.push_back (LineEntry(p0, p1,
-                                                        face->user_flag_set(),
-                                                        face->level() ));
-                         p0=p1;
-                       }
-
-                                                      // generate last piece
-                     const Point<dim> p1_dim (face->vertex(1));
-                     const Point<2>   p1     (p1_dim(0), p1_dim(1));
-                     line_list.push_back (LineEntry(p0, p1,
-                                                    face->user_flag_set(),
-                                                    face->level()));
-                   };
+                           const typename Triangulation<dim>::face_iterator
+                               face = cell->face(face_no);
+
+                           if (face->at_boundary())
+                           {
+                               Point<dim> p0_dim(face->vertex(0));
+                               Point<2>   p0    (p0_dim(0), p0_dim(1));
+
+                               // loop over
+                               // all pieces
+                               // of the line
+                               // and generate
+                               // line-lets
+                               const unsigned int offset=face_no*n_points;
+                               for (unsigned int i=0; i<n_points; ++i)
+                               {
+                                   const Point<dim> p1_dim (mapping->transform_unit_to_real_cell
+                                           (cell, q_projector.point(offset+i)));
+                                   const Point<2>   p1     (p1_dim(0), p1_dim(1));
+
+                                   line_list.push_back (LineEntry(p0, p1,
+                                               face->user_flag_set(),
+                                               face->level() ));
+                                   p0=p1;
+                               }
+
+                               // generate last piece
+                               const Point<dim> p1_dim (face->vertex(1));
+                               const Point<2>   p1     (p1_dim(0), p1_dim(1));
+                               line_list.push_back (LineEntry(p0, p1,
+                                           face->user_flag_set(),
+                                           face->level()));
+                           };
+                       };
                };
-         };
-       
-       break;
-      };
-       
-      case 3:
-      {
-                                        // curved boundary output
-                                        // presently not supported
-//TODO:[RH] curved boundaries in eps for 3d    
-       Assert (mapping == 0, ExcNotImplemented());
-       
-       typename Triangulation<dim>::active_line_iterator
-         line   =tria.begin_active_line (),
-         endline=tria.end_line ();
-       
-                                        // loop over all lines and compute their
-                                        // projection on the plane perpendicular
-                                        // to the direction of sight
-
-                                        // direction of view equals the unit 
-                                        // vector of the position of the
-                                        // spectator to the origin.
-                                        //
-                                        // we chose here the viewpoint as in
-                                        // gnuplot as default.
-                                        //
-//TODO:[WB] Fix a potential problem with viewing angles in 3d Eps GridOut
-                                        // note: the following might be wrong
-                                        // if one of the base vectors below
-                                        // is in direction of the viewer, but
-                                        // I am too tired at present to fix
-                                        // this
-       const double pi = deal_II_numbers::PI;
-       const double z_angle    = eps_flags_3.azimut_angle;
-       const double turn_angle = eps_flags_3.turn_angle;
-       const Point<dim> view_direction(-std::sin(z_angle * 2.*pi / 360.) * std::sin(turn_angle * 2.*pi / 360.),
-                                       +std::sin(z_angle * 2.*pi / 360.) * std::cos(turn_angle * 2.*pi / 360.),
-                                       -std::cos(z_angle * 2.*pi / 360.));
-       
-                                        // decide about the two unit vectors
-                                        // in this plane. we chose the first one
-                                        // to be the projection of the z-axis
-                                        // to this plane
-       const Point<dim> vector1
-         = Point<dim>(0,0,1) - ((Point<dim>(0,0,1) * view_direction) * view_direction);
-       const Point<dim> unit_vector1 = vector1 / std::sqrt(vector1.square());
-       
-                                        // now the third vector is fixed. we
-                                        // chose the projection of a more or
-                                        // less arbitrary vector to the plane
-                                        // perpendicular to the first one
-       const Point<dim> vector2
-         = (Point<dim>(1,0,0)
-            - ((Point<dim>(1,0,0) * view_direction) * view_direction)
-            - ((Point<dim>(1,0,0) * unit_vector1)   * unit_vector1));
-       const Point<dim> unit_vector2 = vector2 / std::sqrt(vector2.square());
-       
-       for (; line!=endline; ++line) 
-         line_list.push_back (LineEntry(Point<2>(line->vertex(0) * unit_vector2,
-                                                 line->vertex(0) * unit_vector1),
-                                        Point<2>(line->vertex(1) * unit_vector2,
-                                                 line->vertex(1) * unit_vector1),
-                                        line->user_flag_set(),
-                                        line->level()));
 
-       break;
-      };
+               break;
+           };
 
-      default:
+       case 3:
+           {
+               // curved boundary output
+               // presently not supported
+               //TODO:[RH] curved boundaries in eps for 3d     
+               Assert (mapping == 0, ExcNotImplemented());
+
+               typename Triangulation<dim>::active_line_iterator
+                   line   =tria.begin_active_line (),
+                          endline=tria.end_line ();
+
+               // loop over all lines and compute their
+               // projection on the plane perpendicular
+               // to the direction of sight
+
+               // direction of view equals the unit 
+               // vector of the position of the
+               // spectator to the origin.
+               //
+               // we chose here the viewpoint as in
+               // gnuplot as default.
+               //
+               //TODO:[WB] Fix a potential problem with viewing angles in 3d Eps GridOut
+               // note: the following might be wrong
+               // if one of the base vectors below
+               // is in direction of the viewer, but
+               // I am too tired at present to fix
+               // this
+               const double pi = deal_II_numbers::PI;
+               const double z_angle    = eps_flags_3.azimut_angle;
+               const double turn_angle = eps_flags_3.turn_angle;
+               const Point<dim> view_direction(-std::sin(z_angle * 2.*pi / 360.) * std::sin(turn_angle * 2.*pi / 360.),
+                       +std::sin(z_angle * 2.*pi / 360.) * std::cos(turn_angle * 2.*pi / 360.),
+                       -std::cos(z_angle * 2.*pi / 360.));
+
+               // decide about the two unit vectors
+               // in this plane. we chose the first one
+               // to be the projection of the z-axis
+               // to this plane
+               const Point<dim> vector1
+                   = Point<dim>(0,0,1) - ((Point<dim>(0,0,1) * view_direction) * view_direction);
+               const Point<dim> unit_vector1 = vector1 / std::sqrt(vector1.square());
+
+               // now the third vector is fixed. we
+               // chose the projection of a more or
+               // less arbitrary vector to the plane
+               // perpendicular to the first one
+               const Point<dim> vector2
+                   = (Point<dim>(1,0,0)
+                           - ((Point<dim>(1,0,0) * view_direction) * view_direction)
+                           - ((Point<dim>(1,0,0) * unit_vector1)   * unit_vector1));
+               const Point<dim> unit_vector2 = vector2 / std::sqrt(vector2.square());
+
+               for (; line!=endline; ++line) 
+                   line_list.push_back (LineEntry(Point<2>(line->vertex(0) * unit_vector2,
+                                   line->vertex(0) * unit_vector1),
+                               Point<2>(line->vertex(1) * unit_vector2,
+                                   line->vertex(1) * unit_vector1),
+                               line->user_flag_set(),
+                               line->level()));
+
+               break;
+           };
+
+       default:
            Assert (false, ExcNotImplemented());
     };
 
 
 
-                                  // find out minimum and maximum x and
-                                  // y coordinates to compute offsets
-                                  // and scaling factors
-  double x_min = tria.begin_active_line()->vertex(0)(0);
-  double x_max = x_min;
-  double y_min = tria.begin_active_line()->vertex(0)(1);
-  double y_max = y_min;
-  unsigned int  max_level = line_list.begin()->level;
+    // find out minimum and maximum x and
+    // y coordinates to compute offsets
+    // and scaling factors
+    double x_min = tria.begin_active_line()->vertex(0)(0);
+    double x_max = x_min;
+    double y_min = tria.begin_active_line()->vertex(0)(1);
+    double y_max = y_min;
+    unsigned int  max_level = line_list.begin()->level;
 
-  for (LineList::const_iterator line=line_list.begin();
-       line!=line_list.end(); ++line)
+    for (LineList::const_iterator line=line_list.begin();
+           line!=line_list.end(); ++line)
     {
-      x_min = std::min (x_min, line->first(0));
-      x_min = std::min (x_min, line->second(0));
+       x_min = std::min (x_min, line->first(0));
+       x_min = std::min (x_min, line->second(0));
 
-      x_max = std::max (x_max, line->first(0));
-      x_max = std::max (x_max, line->second(0));
+       x_max = std::max (x_max, line->first(0));
+       x_max = std::max (x_max, line->second(0));
 
-      y_min = std::min (y_min, line->first(1));
-      y_min = std::min (y_min, line->second(1));
+       y_min = std::min (y_min, line->first(1));
+       y_min = std::min (y_min, line->second(1));
 
-      y_max = std::max (y_max, line->first(1));
-      y_max = std::max (y_max, line->second(1));
-      
-      max_level = std::max (max_level,  line->level);      
+       y_max = std::max (y_max, line->first(1));
+       y_max = std::max (y_max, line->second(1));
+
+       max_level = std::max (max_level,  line->level);      
     };
 
-                                  // scale in x-direction such that
-                                  // in the output 0 <= x <= 300.
-                                  // don't scale in y-direction to
-                                  // preserve the shape of the
-                                  // triangulation
-  const double scale = (eps_flags_base.size /
-                       (eps_flags_base.size_type==GridOutFlags::EpsFlagsBase::width ?
-                        x_max - x_min :
-                        y_min - y_max));
+    // scale in x-direction such that
+    // in the output 0 <= x <= 300.
+    // don't scale in y-direction to
+    // preserve the shape of the
+    // triangulation
+    const double scale = (eps_flags_base.size /
+           (eps_flags_base.size_type==GridOutFlags::EpsFlagsBase::width ?
+            x_max - x_min :
+            y_min - y_max));
 
 
-                                  // now write preamble
-  if (true) 
+    // now write preamble
+    if (true) 
     {
-                                      // block this to have local
-                                      // variables destroyed after
-                                      // use
-      std::time_t  time1= std::time (0);
-      std::tm     *time = std::localtime(&time1); 
-      out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
-         << "%%Title: deal.II Output" << '\n'
-         << "%%Creator: the deal.II library" << '\n'
-         << "%%Creation Date: " 
-         << time->tm_year+1900 << "/"
-         << time->tm_mon+1 << "/"
-         << time->tm_mday << " - "
-         << time->tm_hour << ":"
-         << std::setw(2) << time->tm_min << ":"
-         << std::setw(2) << time->tm_sec << '\n'
-         << "%%BoundingBox: "
-                                        // lower left corner
-         << "0 0 "
-                                        // upper right corner
-         << static_cast<unsigned int>(std::floor(( (x_max-x_min) * scale )+1))
-         << ' '
-          << static_cast<unsigned int>(std::floor(( (y_max-y_min) * scale )+1))
-         << '\n';
-
-                                      // define some abbreviations to keep
-                                      // the output small:
-                                      // m=move turtle to
-                                      // x=execute line stroke
-                                      // b=black pen
-                                      // r=red pen
-      out << "/m {moveto} bind def" << '\n'
-         << "/x {lineto stroke} bind def" << '\n'
-         << "/b {0 0 0 setrgbcolor} def" << '\n'
-         << "/r {1 0 0 setrgbcolor} def" << '\n';
-
-                                      // calculate colors for level
-                                      // coloring; level 0 is black,
-                                      // other levels are blue
-                                      // ... red
-      if (eps_flags_base.color_lines_level)
-       out  << "/l  { neg "
-            << (max_level)
-            << " add "
-            << (0.66666/std::max(1U,(max_level-1)))
-            << " mul 1 0.8 sethsbcolor} def" << '\n';
-
-                                      // in 2d, we can also plot cell
-                                      // and vertex numbers, but this
-                                      // requires a somewhat more
-                                      // lengthy preamble. please
-                                      // don't ask me what most of
-                                      // this means, it is reverse
-                                      // engineered from what GNUPLOT
-                                      // uses in its output
-      if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
-                        eps_flags_2.write_vertex_numbers))
+       // block this to have local
+       // variables destroyed after
+       // use
+       std::time_t  time1= std::time (0);
+       std::tm     *time = std::localtime(&time1); 
+       out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
+           << "%%Title: deal.II Output" << '\n'
+           << "%%Creator: the deal.II library" << '\n'
+           << "%%Creation Date: " 
+           << time->tm_year+1900 << "/"
+           << time->tm_mon+1 << "/"
+           << time->tm_mday << " - "
+           << time->tm_hour << ":"
+           << std::setw(2) << time->tm_min << ":"
+           << std::setw(2) << time->tm_sec << '\n'
+           << "%%BoundingBox: "
+           // lower left corner
+           << "0 0 "
+           // upper right corner
+           << static_cast<unsigned int>(std::floor(( (x_max-x_min) * scale )+1))
+           << ' '
+           << static_cast<unsigned int>(std::floor(( (y_max-y_min) * scale )+1))
+           << '\n';
+
+       // define some abbreviations to keep
+       // the output small:
+       // m=move turtle to
+       // x=execute line stroke
+       // b=black pen
+       // r=red pen
+       out << "/m {moveto} bind def" << '\n'
+           << "/x {lineto stroke} bind def" << '\n'
+           << "/b {0 0 0 setrgbcolor} def" << '\n'
+           << "/r {1 0 0 setrgbcolor} def" << '\n';
+
+       // calculate colors for level
+       // coloring; level 0 is black,
+       // other levels are blue
+       // ... red
+       if (eps_flags_base.color_lines_level)
+           out  << "/l  { neg "
+               << (max_level)
+               << " add "
+               << (0.66666/std::max(1U,(max_level-1)))
+               << " mul 1 0.8 sethsbcolor} def" << '\n';
+
+       // in 2d, we can also plot cell
+       // and vertex numbers, but this
+       // requires a somewhat more
+       // lengthy preamble. please
+       // don't ask me what most of
+       // this means, it is reverse
+       // engineered from what GNUPLOT
+       // uses in its output
+       if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
+                   eps_flags_2.write_vertex_numbers))
        {
-         out << ("/R {rmoveto} bind def\n"
-                 "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
-                 "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
-                 "currentdict end definefont\n"
-                 "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
-                 "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
-                 "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
-                 "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
-                 "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
-                 "5 get stringwidth pop add}\n"
-                 "{pop} ifelse} forall} bind def\n"
-                 "/MCshow { currentpoint stroke m\n"
-                 "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
-             << '\n';
+           out << ("/R {rmoveto} bind def\n"
+                   "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
+                   "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
+                   "currentdict end definefont\n"
+                   "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
+                   "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
+                   "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
+                   "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
+                   "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
+                   "5 get stringwidth pop add}\n"
+                   "{pop} ifelse} forall} bind def\n"
+                   "/MCshow { currentpoint stroke m\n"
+                   "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
+               << '\n';
        };
-      
-      out << "%%EndProlog" << '\n'
-         << '\n';
 
-                                      // set fine lines
-      out << eps_flags_base.line_width << " setlinewidth" << '\n';
+       out << "%%EndProlog" << '\n'
+           << '\n';
+
+       // set fine lines
+       out << eps_flags_base.line_width << " setlinewidth" << '\n';
     };
 
-                                  // now write the lines
-  const Point<2> offset(x_min, y_min);
-  
-  for (LineList::const_iterator line=line_list.begin();
-       line!=line_list.end(); ++line)
-    if (eps_flags_base.color_lines_level && (line->level > 0))
-                                      // lines colored according to
-                                      // refinement level,
-                                      // contributed by Jörg
-                                      // R. Weimar
-      out << line->level
-         << " l "
-         << (line->first  - offset) * scale << " m "
-         << (line->second - offset) * scale << " x" << '\n';
-    else
-      out << ((line->colorize && eps_flags_base.color_lines_on_user_flag) ? "r " : "b ")
-         << (line->first  - offset) * scale << " m "
-         << (line->second - offset) * scale << " x" << '\n';
-    
-                                  // finally write the cell numbers
-                                  // in 2d, if that is desired
-  if ((dim == 2) && (eps_flags_2.write_cell_numbers == true))
+    // now write the lines
+    const Point<2> offset(x_min, y_min);
+
+    for (LineList::const_iterator line=line_list.begin();
+           line!=line_list.end(); ++line)
+       if (eps_flags_base.color_lines_level && (line->level > 0))
+           // lines colored according to
+           // refinement level,
+           // contributed by Jörg
+           // R. Weimar
+           out << line->level
+               << " l "
+               << (line->first  - offset) * scale << " m "
+               << (line->second - offset) * scale << " x" << '\n';
+       else
+           out << ((line->colorize && eps_flags_base.color_lines_on_user_flag) ? "r " : "b ")
+               << (line->first  - offset) * scale << " m "
+               << (line->second - offset) * scale << " x" << '\n';
+
+    // finally write the cell numbers
+    // in 2d, if that is desired
+    if ((dim == 2) && (eps_flags_2.write_cell_numbers == true))
     {
-      out << "(Helvetica) findfont 140 scalefont setfont"
-         << '\n';
-      
-      typename Triangulation<dim>::active_cell_iterator
-       cell = tria.begin_active (),
-       endc = tria.end ();
-      for (; cell!=endc; ++cell)
+       out << "(Helvetica) findfont 140 scalefont setfont"
+           << '\n';
+
+       typename Triangulation<dim>::active_cell_iterator
+           cell = tria.begin_active (),
+                endc = tria.end ();
+       for (; cell!=endc; ++cell)
        {
-         out << (cell->center()(0)-offset(0))*scale << ' '
-             << (cell->center()(1)-offset(1))*scale
-             << " m" << '\n'
-             << "[ [(Helvetica) 12.0 0.0 true true (";
-         if (eps_flags_2.write_cell_number_level)
-           out << cell;
-         else
-           out << cell->index();
-
-         out << ")] "
-             << "] -6 MCshow"
-             << '\n';
+           out << (cell->center()(0)-offset(0))*scale << ' '
+               << (cell->center()(1)-offset(1))*scale
+               << " m" << '\n'
+               << "[ [(Helvetica) 12.0 0.0 true true (";
+           if (eps_flags_2.write_cell_number_level)
+               out << cell;
+           else
+               out << cell->index();
+
+           out << ")] "
+               << "] -6 MCshow"
+               << '\n';
        };
     };
 
-                                  // and the vertex numbers
-  if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
+    // and the vertex numbers
+    if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
     {
-      out << "(Helvetica) findfont 140 scalefont setfont"
-         << '\n';
-
-                                      // have a list of those
-                                      // vertices which we have
-                                      // already tracked, to avoid
-                                      // doing this multiply
-      std::set<unsigned int> treated_vertices;
-      typename Triangulation<dim>::active_cell_iterator
-       cell = tria.begin_active (),
-       endc = tria.end ();
-      for (; cell!=endc; ++cell)
-       for (unsigned int vertex=0;
-            vertex<GeometryInfo<dim>::vertices_per_cell;
-            ++vertex)
-         if (treated_vertices.find(cell->vertex_index(vertex))
-             ==
-             treated_vertices.end())
-           {
-             treated_vertices.insert (cell->vertex_index(vertex));
-             
-             out << (cell->vertex(vertex)(0)-offset(0))*scale << ' '
-                 << (cell->vertex(vertex)(1)-offset(1))*scale
-                 << " m" << '\n'
-                 << "[ [(Helvetica) 10.0 0.0 true true ("
-                 << cell->vertex_index(vertex)
-                 << ")] "
-                 << "] -6 MCshow"
-                 << '\n';
-           };
+       out << "(Helvetica) findfont 140 scalefont setfont"
+           << '\n';
+
+       // have a list of those
+       // vertices which we have
+       // already tracked, to avoid
+       // doing this multiply
+       std::set<unsigned int> treated_vertices;
+       typename Triangulation<dim>::active_cell_iterator
+           cell = tria.begin_active (),
+                endc = tria.end ();
+       for (; cell!=endc; ++cell)
+           for (unsigned int vertex=0;
+                   vertex<GeometryInfo<dim>::vertices_per_cell;
+                   ++vertex)
+               if (treated_vertices.find(cell->vertex_index(vertex))
+                       ==
+                       treated_vertices.end())
+               {
+                   treated_vertices.insert (cell->vertex_index(vertex));
+
+                   out << (cell->vertex(vertex)(0)-offset(0))*scale << ' '
+                       << (cell->vertex(vertex)(1)-offset(1))*scale
+                       << " m" << '\n'
+                       << "[ [(Helvetica) 10.0 0.0 true true ("
+                       << cell->vertex_index(vertex)
+                       << ")] "
+                       << "] -6 MCshow"
+                       << '\n';
+               };
     };
-  
-  out << "showpage" << '\n';
-  
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
 
-  AssertThrow (out, ExcIO());
+    out << "showpage" << '\n';
+
+    // make sure everything now gets to
+    // disk
+    out.flush ();
+
+    AssertThrow (out, ExcIO());
 }
 
 #endif
 
 
-template <int dim>
+    template <int dim>
 void GridOut::write (const Triangulation<dim> &tria,
-                    std::ostream             &out,
-                    const OutputFormat        output_format,
-                    const Mapping<dim>       *mapping)
+       std::ostream             &out,
+       const OutputFormat        output_format,
+       const Mapping<dim>       *mapping)
 {
-  switch (output_format)
+    switch (output_format)
     {
-      case none:
-       return;
-       
-      case dx:
-       write_dx (tria, out);
-       return;
-       
-      case ucd:
-       write_ucd (tria, out);
-       return;
-       
-      case gnuplot:
-       write_gnuplot (tria, out, mapping);
-       return;
-       
-      case eps:
-       write_eps (tria, out, mapping);
-       return;
-       
-      case xfig:
-       write_xfig (tria, out, mapping);
-       return;
+       case none:
+           return;
+
+       case dx:
+           write_dx (tria, out);
+           return;
+
+       case ucd:
+           write_ucd (tria, out);
+           return;
+
+       case gnuplot:
+           write_gnuplot (tria, out, mapping);
+           return;
+
+       case eps:
+           write_eps (tria, out, mapping);
+           return;
+
+       case xfig:
+           write_xfig (tria, out, mapping);
+           return;
+
+       case msh:
+           write_msh (tria, out);
+           return;
     }
-  
-  Assert (false, ExcInternalError());
+
+    Assert (false, ExcInternalError());
 }
 
 
index cc9d12813edb7322b1995ea216f203edcdf11f19..0f2d4d961f2f98b11b3978481860b81fe3a2411f 100644 (file)
@@ -426,6 +426,14 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       New: There is now a <code class="class">GridOut</code>::<code
+       class="member">write_gmsh</code> function to write in the gmsh 
+       format. It has the same features of the write_ucd one. 
+       <br> 
+       (Luca Heltai 2006/02/17)
+       </p>
+
   <li> <p>
        Changed: The <code class="class">hp::FECollection</code>::<code
        class="member">get_fe</code>, <code

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.