]> https://gitweb.dealii.org/ - dealii.git/commitdiff
One more file moved to compile-once mode.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 9 Nov 2010 17:12:56 +0000 (17:12 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 9 Nov 2010 17:12:56 +0000 (17:12 +0000)
git-svn-id: https://svn.dealii.org/trunk@22651 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/grid/grid_out.all_dimensions.cc
deal.II/source/grid/grid_out.cc [deleted file]
deal.II/source/grid/grid_out.inst.in [new file with mode: 0644]

index 880c5d361fe5c49beedaddcfa7135c90f45b4ea2..d95514c16680a1b19fc091e822b94181d7261442 100644 (file)
 
 #include <grid/grid_out.h>
 #include <base/parameter_handler.h>
+#include <base/point.h>
+#include <base/quadrature.h>
+#include <base/qprojector.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <fe/mapping.h>
+
+#include <cstring>
+#include <iomanip>
+#include <algorithm>
+#include <list>
+#include <set>
+#include <ctime>
+#include <cmath>
+
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -546,4 +562,1895 @@ GridOut::memory_consumption () const
          sizeof(xfig_flags));
 }
 
+
+
+template <>
+void GridOut::write_dx (const Triangulation<1> &,
+                       std::ostream             &) const
+{
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+template <int dim>
+void GridOut::write_dx (const Triangulation<dim> &tria,
+                       std::ostream             &out) const
+{
+//TODO:[GK] allow for boundary faces only
+  Assert(dx_flags.write_all_faces, ExcNotImplemented());
+  AssertThrow (out, ExcIO());
+                                  // Copied and adapted from write_ucd
+  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();
+
+                                  // vertices are implicitly numbered from 0 to
+                                  // n_vertices-1. we have to renumber the
+                                  // vertices, because otherwise we would end
+                                  // up with wrong results, if there are unused
+                                  // vertices
+  std::vector<unsigned int> renumber(vertices.size());
+                                  // fill this vector with new vertex numbers
+                                  // ranging from 0 to n_vertices-1
+  unsigned int new_number=0;
+  for (unsigned int i=0; i<vertices.size(); ++i)
+    if (vertex_used[i])
+      renumber[i]=new_number++;
+  Assert(new_number==n_vertices, ExcInternalError());
+
+  typename Triangulation<dim>::active_cell_iterator       cell;
+  const typename Triangulation<dim>::active_cell_iterator endc=tria.end();
+
+
+                                  // write the vertices
+  out << "object \"vertices\" class array type float rank 1 shape " << dim
+      << " items " << n_vertices << " data follows"
+      << '\n';
+
+  for (unsigned int i=0; i<vertices.size(); ++i)
+    if (vertex_used[i])
+      {
+       out << '\t' << vertices[i] << '\n';
+      };
+
+                                  // write cells or faces
+  const bool write_cells = dx_flags.write_cells;
+  const bool write_faces = (dim>1) ? dx_flags.write_faces : false;
+
+  const unsigned int n_cells = tria.n_active_cells();
+  const unsigned int n_faces = tria.n_active_cells()
+                              * GeometryInfo<dim>::faces_per_cell;
+
+  const unsigned int n_vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
+  const unsigned int n_vertices_per_face = GeometryInfo<dim>::vertices_per_face;
+
+  if (write_cells)
+    {
+      out << "object \"cells\" class array type int rank 1 shape "
+         << n_vertices_per_cell
+         << " items " << n_cells << " data follows" << '\n';
+
+      for (cell = tria.begin_active(); cell != endc; ++cell)
+       {
+         for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+           out << '\t' << renumber[cell->vertex_index(GeometryInfo<dim>::dx_to_deal[v])];
+         out << '\n';
+       }
+      out << "attribute \"element type\" string \"";
+      if (dim==1) out << "lines";
+      if (dim==2) out << "quads";
+      if (dim==3) out << "cubes";
+      out << "\"" << '\n'
+         << "attribute \"ref\" string \"positions\"" << '\n' << '\n';
+
+                                      // Additional cell information
+
+      out << "object \"material\" class array type int rank 0 items "
+         << n_cells << " data follows" << '\n';
+      for (cell = tria.begin_active(); cell != endc; ++cell)
+       out << ' ' << (unsigned int)cell->material_id();
+      out  << '\n'
+          << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
+
+      out << "object \"level\" class array type int rank 0 items "
+         << n_cells << " data follows" << '\n';
+      for (cell = tria.begin_active(); cell != endc; ++cell)
+       out << ' ' << cell->level();
+      out  << '\n'
+          << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
+
+      if (dx_flags.write_measure)
+       {
+         out << "object \"measure\" class array type float rank 0 items "
+             << n_cells << " data follows" << '\n';
+         for (cell = tria.begin_active(); cell != endc; ++cell)
+           out << '\t' << cell->measure();
+         out  << '\n'
+              << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
+       }
+
+      if (dx_flags.write_diameter)
+       {
+         out << "object \"diameter\" class array type float rank 0 items "
+             << n_cells << " data follows" << '\n';
+         for (cell = tria.begin_active(); cell != endc; ++cell)
+           out << '\t' << cell->diameter();
+         out  << '\n'
+              << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
+       }
+    }
+
+  if (write_faces)
+    {
+      out << "object \"faces\" class array type int rank 1 shape "
+         << n_vertices_per_face
+         << " items " << n_faces << " data follows"
+         << '\n';
+
+      for (cell = tria.begin_active(); cell != endc; ++cell)
+       {
+         for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
+           {
+             typename Triangulation<dim>::face_iterator face = cell->face(f);
+
+             for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+               out << '\t' << renumber[face->vertex_index(GeometryInfo<dim-1>::dx_to_deal[v])];
+             out << '\n';
+           }
+       }
+      out << "attribute \"element type\" string \"";
+      if (dim==2) out << "lines";
+      if (dim==3) out << "quads";
+      out << "\"" << '\n'
+         << "attribute \"ref\" string \"positions\"" << '\n' << '\n';
+
+
+                                      // Additional face information
+
+      out << "object \"boundary\" class array type int rank 0 items "
+         << n_faces << " data follows" << '\n';
+      for (cell = tria.begin_active(); cell != endc; ++cell)
+       {
+                                          // Little trick to get -1
+                                          // for the interior
+         for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
+           out << ' ' << (int)(signed char)cell->face(f)->boundary_indicator();
+         out << '\n';
+       }
+      out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
+
+      if (dx_flags.write_measure)
+       {
+         out << "object \"face measure\" class array type float rank 0 items "
+             << n_faces << " data follows" << '\n';
+         for (cell = tria.begin_active(); cell != endc; ++cell)
+           {
+             for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
+               out << ' ' << cell->face(f)->measure();
+             out << '\n';
+           }
+         out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
+       }
+
+      if (dx_flags.write_diameter)
+       {
+         out << "object \"face diameter\" class array type float rank 0 items "
+             << n_faces << " data follows" << '\n';
+         for (cell = tria.begin_active(); cell != endc; ++cell)
+           {
+             for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
+               out << ' ' << cell->face(f)->diameter();
+             out << '\n';
+           }
+         out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
+       }
+
+    }
+
+
+                                  // Write additional face information
+
+  if (write_faces)
+    {
+
+    }
+  else
+    {
+    }
+
+                                  // The wrapper
+  out << "object \"deal data\" class field" << '\n'
+      << "component \"positions\" value \"vertices\"" << '\n'
+      << "component \"connections\" value \"cells\"" << '\n';
+
+  if (write_cells)
+    {
+      out << "object \"cell data\" class field" << '\n'
+         << "component \"positions\" value \"vertices\"" << '\n'
+         << "component \"connections\" value \"cells\"" << '\n';
+      out << "component \"material\" value \"material\"" << '\n';
+      out << "component \"level\" value \"level\"" << '\n';
+      if (dx_flags.write_measure)
+       out << "component \"measure\" value \"measure\"" << '\n';
+      if (dx_flags.write_diameter)
+       out << "component \"diameter\" value \"diameter\"" << '\n';
+    }
+
+  if (write_faces)
+    {
+      out << "object \"face data\" class field" << '\n'
+         << "component \"positions\" value \"vertices\"" << '\n'
+         << "component \"connections\" value \"faces\"" << '\n';
+      out << "component \"boundary\" value \"boundary\"" << '\n';
+      if (dx_flags.write_measure)
+       out << "component \"measure\" value \"face measure\"" << '\n';
+      if (dx_flags.write_diameter)
+       out << "component \"diameter\" value \"face diameter\"" << '\n';
+    }
+
+  out << '\n'
+      << "object \"grid data\" class group" << '\n';
+  if (write_cells)
+    out << "member \"cells\" value \"cell data\"" << '\n';
+  if (write_faces)
+    out << "member \"faces\" value \"face data\"" << '\n';
+  out << "end" << '\n';
+
+                                  // make sure everything now gets to
+                                  // disk
+  out.flush ();
+
+  AssertThrow (out, ExcIO());
+}
+
+
+
+
+template <int dim, int spacedim>
+void GridOut::write_msh (const Triangulation<dim, spacedim> &tria,
+                        std::ostream             &out) const
+{
+  AssertThrow (out, ExcIO());
+
+                                  // get the positions of the
+                                  // vertices and whether they are
+                                  // used.
+  const std::vector<Point<spacedim> > &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,spacedim>::active_cell_iterator       cell=tria.begin_active();
+  const typename Triangulation<dim,spacedim>::active_cell_iterator endc=tria.end();
+
+                                  // Write Header
+                                  // The file format is:
+                                  /*
+
+
+                                  $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])
+      {
+       out << i+1                 // vertex index
+           << "  "
+           << vertices[i];
+       for (unsigned int d=spacedim+1; d<=3; ++d)
+         out << " 0";             // fill with zeroes
+       out << std::endl;
+      };
+
+                                  // Write cells preamble
+  out << "$ENDNOD" << std::endl
+      << "$ELM" << std::endl
+      << tria.n_active_cells() + ( (msh_flags.write_faces ?
+                                   n_boundary_faces(tria) :0 ) +
+                                  ( msh_flags.write_lines ?
+                                    n_boundary_lines(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)
+    {
+      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 and lines with
+                                  // non-zero boundary indicator
+  if (msh_flags.write_faces)
+    write_msh_faces (tria, cell_index, out);
+  if (msh_flags.write_lines)
+    write_msh_lines (tria, cell_index, out);
+
+  out << "$ENDELM" << std::endl;
+
+                                  // make sure everything now gets to
+                                  // disk
+  out.flush ();
+
+  AssertThrow (out, ExcIO());
+}
+
+
+template <int dim, int spacedim>
+void GridOut::write_ucd (const Triangulation<dim,spacedim> &tria,
+                        std::ostream             &out) const
+{
+  AssertThrow (out, ExcIO());
+
+                                  // get the positions of the
+                                  // vertices and whether they are
+                                  // used.
+  const std::vector<Point<spacedim> > &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,spacedim>::active_cell_iterator       cell=tria.begin_active();
+  const typename Triangulation<dim,spacedim>::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) +
+                                  (ucd_flags.write_lines ?
+                                   n_boundary_lines(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])
+      {
+       out << i+1                 // vertex index
+           << "  "
+           << vertices[i];
+       for (unsigned int d=spacedim+1; d<=3; ++d)
+         out << " 0";             // fill with zeroes
+       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 and lines with
+                                  // non-zero boundary indicator
+  if (ucd_flags.write_faces)
+    write_ucd_faces (tria, cell_index, out);
+  if (ucd_flags.write_lines)
+    write_ucd_lines (tria, cell_index, out);
+
+                                  // make sure everything now gets to
+                                  // disk
+  out.flush ();
+
+  AssertThrow (out, ExcIO());
+}
+
+
+
+template <int dim>
+void GridOut::write_xfig (
+  const Triangulation<dim>&,
+  std::ostream&,
+  const Mapping<dim>*) const
+{
+  Assert (false, ExcNotImplemented());
+}
+
+
+//TODO:[GK] Obey parameters
+//TODO:[GK] Flip y-axis?
+template <>
+void GridOut::write_xfig (
+  const Triangulation<2>& tria,
+  std::ostream&             out,
+  const Mapping<2>*       /*mapping*/) const
+{
+  const int dim = 2;
+
+  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)
+       {
+         const Point<dim>& p = cell->vertex(
+           GeometryInfo<dim>::ucd_to_deal[k % nv]);
+         for (unsigned int d=0; d<static_cast<unsigned int>(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<static_cast<unsigned int>(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 ();
+
+  AssertThrow (out, ExcIO());
+}
+
+
+
+
+
+unsigned int GridOut::n_boundary_faces (const Triangulation<1> &) const
+{
+  return 0;
+}
+
+unsigned int GridOut::n_boundary_lines (const Triangulation<1> &) const
+{
+  return 0;
+}
+
+
+unsigned int GridOut::n_boundary_faces (const Triangulation<1,2> &) const
+{
+  return 0;
+}
+
+unsigned int GridOut::n_boundary_lines (const Triangulation<1,2> &) const
+{
+  return 0;
+}
+
+unsigned int GridOut::n_boundary_lines (const Triangulation<2> &) const
+{
+  return 0;
+}
+
+unsigned int GridOut::n_boundary_lines (const Triangulation<2,3> &) const
+{
+  return 0;
+}
+
+
+
+template <int dim, int spacedim>
+unsigned int GridOut::n_boundary_faces (const Triangulation<dim,spacedim> &tria) const
+{
+    typename Triangulation<dim,spacedim>::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++;
+
+  return n_faces;
+}
+
+
+
+template <int dim, int spacedim>
+unsigned int GridOut::n_boundary_lines (const Triangulation<dim, spacedim> &tria) const
+{
+    typename Triangulation<dim, spacedim>::active_line_iterator edge, endedge;
+  unsigned int n_lines = 0;
+
+  for (edge=tria.begin_active_line(), endedge=tria.end_line();
+       edge != endedge; ++edge)
+    if ((edge->at_boundary()) &&
+       (edge->boundary_indicator() != 0))
+      n_lines++;
+
+  return n_lines;
+}
+
+
+
+
+void GridOut::write_msh_faces (const Triangulation<1> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+
+void GridOut::write_msh_faces (const Triangulation<1,2> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+
+void GridOut::write_msh_lines (const Triangulation<1> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+void GridOut::write_msh_lines (const Triangulation<1,2> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+
+void GridOut::write_msh_lines (const Triangulation<2> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+void GridOut::write_msh_lines (const Triangulation<2,3> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+
+
+
+template <int dim, int spacedim>
+void GridOut::write_msh_faces (const Triangulation<dim,spacedim> &tria,
+                              const unsigned int        starting_index,
+                              std::ostream             &out) const
+{
+    typename Triangulation<dim,spacedim>::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())
+           << ' '
+           << static_cast<unsigned int>(face->boundary_indicator())
+           << ' ' << 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;
+      };
+}
+
+
+template <int dim, int spacedim>
+void GridOut::write_msh_lines (const Triangulation<dim, spacedim> &tria,
+                              const unsigned int        starting_index,
+                              std::ostream             &out) const
+{
+    typename Triangulation<dim,spacedim>::active_line_iterator line, endl;
+  unsigned int index=starting_index;
+
+  for (line=tria.begin_active_line(), endl=tria.end_line();
+       line != endl; ++line)
+    if (line->at_boundary() &&
+       (line->boundary_indicator() != 0))
+      {
+       out << index << " 1 ";
+       out << static_cast<unsigned int>(line->boundary_indicator())
+           << ' '
+           << static_cast<unsigned int>(line->boundary_indicator())
+           << " 2 ";
+                                        // note: vertex numbers are 1-base
+       for (unsigned int vertex=0; vertex<2; ++vertex)
+         out << ' '
+             << line->vertex_index(GeometryInfo<dim-2>::ucd_to_deal[vertex])+1;
+       out << '\n';
+
+       ++index;
+      };
+}
+
+
+
+
+void GridOut::write_ucd_faces (const Triangulation<1> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+void GridOut::write_ucd_faces (const Triangulation<1,2> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+void GridOut::write_ucd_lines (const Triangulation<1> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+void GridOut::write_ucd_lines (const Triangulation<1,2> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+
+void GridOut::write_ucd_lines (const Triangulation<2> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+void GridOut::write_ucd_lines (const Triangulation<2,3> &,
+                              const unsigned int,
+                              std::ostream &) const
+{
+  return;
+}
+
+
+
+
+template <int dim, int spacedim>
+void GridOut::write_ucd_faces (const Triangulation<dim, spacedim> &tria,
+                              const unsigned int        starting_index,
+                              std::ostream             &out) const
+{
+    typename Triangulation<dim,spacedim>::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';
+
+       ++index;
+      };
+}
+
+
+
+template <int dim, int spacedim>
+void GridOut::write_ucd_lines (const Triangulation<dim, spacedim> &tria,
+                              const unsigned int        starting_index,
+                              std::ostream             &out) const
+{
+    typename Triangulation<dim, spacedim>::active_line_iterator line, endl;
+  unsigned int index=starting_index;
+
+  for (line=tria.begin_active_line(), endl=tria.end_line();
+       line != endl; ++line)
+    if (line->at_boundary() &&
+       (line->boundary_indicator() != 0))
+      {
+       out << index << "  "
+           << static_cast<unsigned int>(line->boundary_indicator())
+           << "  line    ";
+       // note: vertex numbers are 1-base
+       for (unsigned int vertex=0; vertex<2; ++vertex)
+         out << line->vertex_index(GeometryInfo<dim-2>::ucd_to_deal[vertex])+1 << ' ';
+       out << '\n';
+
+       ++index;
+      };
+}
+
+
+namespace internal
+{
+  namespace
+  {
+    template <int spacedim>
+    void write_gnuplot (const dealii::Triangulation<1,spacedim> &tria,
+                       std::ostream             &out,
+                       const Mapping<1,spacedim>       *,
+                       const GridOutFlags::Gnuplot &gnuplot_flags)
+    {
+      AssertThrow (out, ExcIO());
+
+      const int dim = 1;
+
+      typename dealii::Triangulation<dim,spacedim>::active_cell_iterator
+       cell=tria.begin_active();
+      const typename dealii::Triangulation<dim,spacedim>::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';
+       }
+
+                                      // make sure everything now gets to
+                                      // disk
+      out.flush ();
+
+      AssertThrow (out, ExcIO());
+    }
+
+
+
+    template <int spacedim>
+    void write_gnuplot (const dealii::Triangulation<2,spacedim> &tria,
+                       std::ostream           &out,
+                       const Mapping<2,spacedim>       *mapping,
+                       const GridOutFlags::Gnuplot &gnuplot_flags)
+    {
+      AssertThrow (out, ExcIO());
+
+      const int dim = 2;
+
+      const unsigned int n_additional_points=
+       gnuplot_flags.n_boundary_face_points;
+      const unsigned int n_points=2+n_additional_points;
+
+      typename dealii::Triangulation<dim,spacedim>::active_cell_iterator
+       cell=tria.begin_active();
+      const typename dealii::Triangulation<dim,spacedim>::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<dim-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<dim-1> quadrature(boundary_points, dummy_weights);
+
+         q_projector = new Quadrature<dim> (QProjector<dim>::project_to_all_faces(quadrature));
+       }
+
+      for (; cell!=endc; ++cell)
+       {
+         if (gnuplot_flags.write_cell_numbers)
+           out << "# cell " << cell << '\n';
+
+         if (mapping==0 ||
+             (!cell->at_boundary() && !gnuplot_flags.curved_inner_cells))
+           {
+                                              // 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'
+                 << '\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
+           {
+             for (unsigned int face_no=0;
+                  face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+               {
+                 const typename dealii::Triangulation<dim,spacedim>::face_iterator
+                   face = cell->face(face_no);
+                 if (face->at_boundary() || gnuplot_flags.curved_inner_cells)
+                   {
+                                                      // 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';
+                   }
+                 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 (q_projector != 0)
+       delete q_projector;
+
+                                      // make sure everything now gets to
+                                      // disk
+      out.flush ();
+
+      AssertThrow (out, ExcIO());
+    }
+
+
+
+    template <int spacedim>
+    void write_gnuplot (const dealii::Triangulation<3,spacedim> &tria,
+                       std::ostream           &out,
+                       const Mapping<3,spacedim>       *mapping,
+                       const GridOutFlags::Gnuplot &gnuplot_flags)
+    {
+      AssertThrow (out, ExcIO());
+
+      const int dim = 3;
+
+      const unsigned int n_additional_points=
+       gnuplot_flags.n_boundary_face_points;
+      const unsigned int n_points=2+n_additional_points;
+
+      typename dealii::Triangulation<dim,spacedim>::active_cell_iterator
+       cell=tria.begin_active();
+      const typename dealii::Triangulation<dim,spacedim>::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));
+       }
+
+      for (; cell!=endc; ++cell)
+       {
+         if (gnuplot_flags.write_cell_numbers)
+           out << "# cell " << cell << '\n';
+
+         if (mapping==0 || n_points==2 ||
+             (!cell->has_boundary_lines() && !gnuplot_flags.curved_inner_cells))
+           {
+                                              // 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
+           {
+             for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+               {
+                 const typename dealii::Triangulation<dim,spacedim>::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<spacedim> 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
+                   {
+                     for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
+                       {
+                         const typename dealii::Triangulation<dim,spacedim>::line_iterator
+                           line=face->line(l);
+
+                         const Point<spacedim> &v0=line->vertex(0),
+                                               &v1=line->vertex(1);
+                         if (line->at_boundary() || gnuplot_flags.curved_inner_cells)
+                           {
+                                                              // transform_real_to_unit_cell
+                                                              // could be
+                                                              // replaced
+                                                              // by using
+                                                              // QProjector<dim>::project_to_line
+                                                              // which is
+                                                              // not yet
+                                                              // implemented
+                             const Point<spacedim> 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';
+
+                         out << '\n' << '\n';
+                       }
+                   }
+               }
+           }
+       }
+
+      if (q_projector != 0)
+       delete q_projector;
+
+
+                                      // make sure everything now gets to
+                                      // disk
+      out.flush ();
+
+      AssertThrow (out, ExcIO());
+    }
+  }
+}
+
+
+
+template <int dim, int spacedim>
+void GridOut::write_gnuplot (
+  const Triangulation<dim,spacedim> &tria,
+  std::ostream             &out,
+  const Mapping<dim,spacedim>       *mapping) const
+{
+  internal::write_gnuplot (tria, out, mapping, gnuplot_flags);
+}
+
+
+
+namespace internal
+{
+  namespace
+  {
+    struct LineEntry
+    {
+       Point<2> first;
+       Point<2> second;
+       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)
+         {}
+    };
+
+
+    void write_eps (const dealii::Triangulation<1> &,
+                   std::ostream &,
+                   const Mapping<1> *,
+                   const GridOutFlags::Eps<2> &,
+                   const GridOutFlags::Eps<3> &)
+    {
+      Assert(false, ExcNotImplemented());
+    }
+
+
+
+    template <int dim>
+    void write_eps (const dealii::Triangulation<dim> &tria,
+                   std::ostream             &out,
+                   const Mapping<dim>       *mapping,
+                   const GridOutFlags::Eps<2> &eps_flags_2,
+                   const GridOutFlags::Eps<3> &eps_flags_3)
+    {
+      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<const GridOutFlags::EpsFlagsBase&>(eps_flags_2) :
+                          (dim==3 ?
+                           static_cast<const GridOutFlags::EpsFlagsBase&>(eps_flags_3) :
+                           *static_cast<const 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 dealii::Triangulation<dim>::active_cell_iterator
+             cell=tria.begin_active(),
+             endc=tria.end();
+
+           for (; cell!=endc; ++cell)
+             for (unsigned int line_no=0;
+                  line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
+               {
+                 typename dealii::Triangulation<dim>::line_iterator
+                   line=cell->line(line_no);
+
+                                                  // 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
+                 if (!line->has_children() &&
+                     (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(),
+                                                  cell->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 dealii::Triangulation<dim>::active_cell_iterator
+                 cell=tria.begin_active ();
+               const typename dealii::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 typename dealii::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(),
+                                                            cell->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(),
+                                                        cell->level()));
+                       };
+                   };
+             };
+
+           break;
+         };
+
+         case 3:
+         {
+                                            // curved boundary output
+                                            // presently not supported
+           Assert (mapping == 0, ExcNotImplemented());
+
+           typename dealii::Triangulation<dim>::active_cell_iterator
+             cell=tria.begin_active(),
+             endc=tria.end();
+
+                                            // 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 = 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 (; cell!=endc; ++cell)
+             for (unsigned int line_no=0;
+                  line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
+               {
+                 typename dealii::Triangulation<dim>::line_iterator
+                   line=cell->line(line_no);
+                 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(),
+                                                cell->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;
+
+      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_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_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));
+
+
+                                      // 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))
+           {
+             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';
+       };
+
+                                      // 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 dealii::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';
+           };
+       };
+
+                                      // 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 dealii::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());
+    }
+  }
+}
+
+
+template <int dim>
+void GridOut::write_eps (const Triangulation<dim> &tria,
+                        std::ostream             &out,
+                        const Mapping<dim>       *mapping) const
+{
+  internal::write_eps (tria, out, mapping,
+                      eps_flags_2, eps_flags_3);
+}
+
+
+template <int dim, int spacedim>
+void GridOut::write (const Triangulation<dim, spacedim> &tria,
+                    std::ostream             &out,
+                    const OutputFormat        output_format,
+                    const Mapping<dim,spacedim>       *mapping) const
+{
+  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 msh:
+           write_msh (tria, out);
+           return;
+    }
+
+  Assert (false, ExcInternalError());
+}
+
+
+template <int dim, int spacedim>
+void GridOut::write (
+    const Triangulation<dim, spacedim> &tria,
+  std::ostream             &out,
+  const Mapping<dim,spacedim>       *mapping) const
+{
+  write(tria, out, default_format, mapping);
+}
+
+
+// explicit instantiations
+#include "grid_out.inst"
+
+
 DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/source/grid/grid_out.cc b/deal.II/source/grid/grid_out.cc
deleted file mode 100644 (file)
index e843f63..0000000
+++ /dev/null
@@ -1,1957 +0,0 @@
-//---------------------------------------------------------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2010 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//---------------------------------------------------------------------------
-
-
-#include <base/point.h>
-#include <base/quadrature.h>
-#include <base/qprojector.h>
-#include <grid/grid_out.h>
-#include <grid/tria.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_iterator.h>
-#include <fe/mapping.h>
-
-#include <cstring>
-#include <iomanip>
-#include <algorithm>
-#include <list>
-#include <set>
-#include <ctime>
-#include <cmath>
-
-DEAL_II_NAMESPACE_OPEN
-
-
-#if deal_II_dimension == 1
-
-
-template <int dim>
-void GridOut::write_dx (const Triangulation<dim> &,
-                       std::ostream             &) const
-{
-  Assert (false, ExcNotImplemented());
-}
-
-#else
-
-
-template <int dim>
-void GridOut::write_dx (const Triangulation<dim> &tria,
-                       std::ostream             &out) const
-{
-//TODO:[GK] allow for boundary faces only
-  Assert(dx_flags.write_all_faces, ExcNotImplemented());
-  AssertThrow (out, ExcIO());
-                                  // Copied and adapted from write_ucd
-  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();
-
-                                  // vertices are implicitly numbered from 0 to
-                                  // n_vertices-1. we have to renumber the
-                                  // vertices, because otherwise we would end
-                                  // up with wrong results, if there are unused
-                                  // vertices
-  std::vector<unsigned int> renumber(vertices.size());
-                                  // fill this vector with new vertex numbers
-                                  // ranging from 0 to n_vertices-1
-  unsigned int new_number=0;
-  for (unsigned int i=0; i<vertices.size(); ++i)
-    if (vertex_used[i])
-      renumber[i]=new_number++;
-  Assert(new_number==n_vertices, ExcInternalError());
-
-  typename Triangulation<dim>::active_cell_iterator       cell;
-  const typename Triangulation<dim>::active_cell_iterator endc=tria.end();
-
-
-                                  // write the vertices
-  out << "object \"vertices\" class array type float rank 1 shape " << dim
-      << " items " << n_vertices << " data follows"
-      << '\n';
-
-  for (unsigned int i=0; i<vertices.size(); ++i)
-    if (vertex_used[i])
-      {
-       out << '\t' << vertices[i] << '\n';
-      };
-
-                                  // write cells or faces
-  const bool write_cells = dx_flags.write_cells;
-  const bool write_faces = (dim>1) ? dx_flags.write_faces : false;
-
-  const unsigned int n_cells = tria.n_active_cells();
-  const unsigned int n_faces = tria.n_active_cells()
-                              * GeometryInfo<dim>::faces_per_cell;
-
-  const unsigned int n_vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
-  const unsigned int n_vertices_per_face = GeometryInfo<dim>::vertices_per_face;
-
-  if (write_cells)
-    {
-      out << "object \"cells\" class array type int rank 1 shape "
-         << n_vertices_per_cell
-         << " items " << n_cells << " data follows" << '\n';
-
-      for (cell = tria.begin_active(); cell != endc; ++cell)
-       {
-         for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-           out << '\t' << renumber[cell->vertex_index(GeometryInfo<dim>::dx_to_deal[v])];
-         out << '\n';
-       }
-      out << "attribute \"element type\" string \"";
-      if (dim==1) out << "lines";
-      if (dim==2) out << "quads";
-      if (dim==3) out << "cubes";
-      out << "\"" << '\n'
-         << "attribute \"ref\" string \"positions\"" << '\n' << '\n';
-
-                                      // Additional cell information
-
-      out << "object \"material\" class array type int rank 0 items "
-         << n_cells << " data follows" << '\n';
-      for (cell = tria.begin_active(); cell != endc; ++cell)
-       out << ' ' << (unsigned int)cell->material_id();
-      out  << '\n'
-          << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
-
-      out << "object \"level\" class array type int rank 0 items "
-         << n_cells << " data follows" << '\n';
-      for (cell = tria.begin_active(); cell != endc; ++cell)
-       out << ' ' << cell->level();
-      out  << '\n'
-          << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
-
-      if (dx_flags.write_measure)
-       {
-         out << "object \"measure\" class array type float rank 0 items "
-             << n_cells << " data follows" << '\n';
-         for (cell = tria.begin_active(); cell != endc; ++cell)
-           out << '\t' << cell->measure();
-         out  << '\n'
-              << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
-       }
-
-      if (dx_flags.write_diameter)
-       {
-         out << "object \"diameter\" class array type float rank 0 items "
-             << n_cells << " data follows" << '\n';
-         for (cell = tria.begin_active(); cell != endc; ++cell)
-           out << '\t' << cell->diameter();
-         out  << '\n'
-              << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
-       }
-    }
-
-  if (write_faces)
-    {
-      out << "object \"faces\" class array type int rank 1 shape "
-         << n_vertices_per_face
-         << " items " << n_faces << " data follows"
-         << '\n';
-
-      for (cell = tria.begin_active(); cell != endc; ++cell)
-       {
-         for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
-           {
-             typename Triangulation<dim>::face_iterator face = cell->face(f);
-
-             for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
-               out << '\t' << renumber[face->vertex_index(GeometryInfo<dim-1>::dx_to_deal[v])];
-             out << '\n';
-           }
-       }
-      out << "attribute \"element type\" string \"";
-      if (dim==2) out << "lines";
-      if (dim==3) out << "quads";
-      out << "\"" << '\n'
-         << "attribute \"ref\" string \"positions\"" << '\n' << '\n';
-
-
-                                      // Additional face information
-
-      out << "object \"boundary\" class array type int rank 0 items "
-         << n_faces << " data follows" << '\n';
-      for (cell = tria.begin_active(); cell != endc; ++cell)
-       {
-                                          // Little trick to get -1
-                                          // for the interior
-         for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
-           out << ' ' << (int)(signed char)cell->face(f)->boundary_indicator();
-         out << '\n';
-       }
-      out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
-
-      if (dx_flags.write_measure)
-       {
-         out << "object \"face measure\" class array type float rank 0 items "
-             << n_faces << " data follows" << '\n';
-         for (cell = tria.begin_active(); cell != endc; ++cell)
-           {
-             for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
-               out << ' ' << cell->face(f)->measure();
-             out << '\n';
-           }
-         out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
-       }
-
-      if (dx_flags.write_diameter)
-       {
-         out << "object \"face diameter\" class array type float rank 0 items "
-             << n_faces << " data follows" << '\n';
-         for (cell = tria.begin_active(); cell != endc; ++cell)
-           {
-             for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
-               out << ' ' << cell->face(f)->diameter();
-             out << '\n';
-           }
-         out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
-       }
-
-    }
-
-
-                                  // Write additional face information
-
-  if (write_faces)
-    {
-
-    }
-  else
-    {
-    }
-
-                                  // The wrapper
-  out << "object \"deal data\" class field" << '\n'
-      << "component \"positions\" value \"vertices\"" << '\n'
-      << "component \"connections\" value \"cells\"" << '\n';
-
-  if (write_cells)
-    {
-      out << "object \"cell data\" class field" << '\n'
-         << "component \"positions\" value \"vertices\"" << '\n'
-         << "component \"connections\" value \"cells\"" << '\n';
-      out << "component \"material\" value \"material\"" << '\n';
-      out << "component \"level\" value \"level\"" << '\n';
-      if (dx_flags.write_measure)
-       out << "component \"measure\" value \"measure\"" << '\n';
-      if (dx_flags.write_diameter)
-       out << "component \"diameter\" value \"diameter\"" << '\n';
-    }
-
-  if (write_faces)
-    {
-      out << "object \"face data\" class field" << '\n'
-         << "component \"positions\" value \"vertices\"" << '\n'
-         << "component \"connections\" value \"faces\"" << '\n';
-      out << "component \"boundary\" value \"boundary\"" << '\n';
-      if (dx_flags.write_measure)
-       out << "component \"measure\" value \"face measure\"" << '\n';
-      if (dx_flags.write_diameter)
-       out << "component \"diameter\" value \"face diameter\"" << '\n';
-    }
-
-  out << '\n'
-      << "object \"grid data\" class group" << '\n';
-  if (write_cells)
-    out << "member \"cells\" value \"cell data\"" << '\n';
-  if (write_faces)
-    out << "member \"faces\" value \"face data\"" << '\n';
-  out << "end" << '\n';
-
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
-
-  AssertThrow (out, ExcIO());
-}
-
-
-#endif
-
-
-template <int dim, int spacedim>
-void GridOut::write_msh (const Triangulation<dim, spacedim> &tria,
-                        std::ostream             &out) const
-{
-  AssertThrow (out, ExcIO());
-
-                                  // get the positions of the
-                                  // vertices and whether they are
-                                  // used.
-  const std::vector<Point<spacedim> > &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,spacedim>::active_cell_iterator       cell=tria.begin_active();
-  const typename Triangulation<dim,spacedim>::active_cell_iterator endc=tria.end();
-
-                                  // Write Header
-                                  // The file format is:
-                                  /*
-
-
-                                  $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])
-      {
-       out << i+1                 // vertex index
-           << "  "
-           << vertices[i];
-       for (unsigned int d=spacedim+1; d<=3; ++d)
-         out << " 0";             // fill with zeroes
-       out << std::endl;
-      };
-
-                                  // Write cells preamble
-  out << "$ENDNOD" << std::endl
-      << "$ELM" << std::endl
-      << tria.n_active_cells() + ( (msh_flags.write_faces ?
-                                   n_boundary_faces(tria) :0 ) +
-                                  ( msh_flags.write_lines ?
-                                    n_boundary_lines(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)
-    {
-      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 and lines with
-                                  // non-zero boundary indicator
-  if (msh_flags.write_faces)
-    write_msh_faces (tria, cell_index, out);
-  if (msh_flags.write_lines)
-    write_msh_lines (tria, cell_index, out);
-
-  out << "$ENDELM" << std::endl;
-
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
-
-  AssertThrow (out, ExcIO());
-}
-
-
-template <int dim, int spacedim>
-void GridOut::write_ucd (const Triangulation<dim,spacedim> &tria,
-                        std::ostream             &out) const
-{
-  AssertThrow (out, ExcIO());
-
-                                  // get the positions of the
-                                  // vertices and whether they are
-                                  // used.
-  const std::vector<Point<spacedim> > &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,spacedim>::active_cell_iterator       cell=tria.begin_active();
-  const typename Triangulation<dim,spacedim>::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) +
-                                  (ucd_flags.write_lines ?
-                                   n_boundary_lines(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])
-      {
-       out << i+1                 // vertex index
-           << "  "
-           << vertices[i];
-       for (unsigned int d=spacedim+1; d<=3; ++d)
-         out << " 0";             // fill with zeroes
-       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 and lines with
-                                  // non-zero boundary indicator
-  if (ucd_flags.write_faces)
-    write_ucd_faces (tria, cell_index, out);
-  if (ucd_flags.write_lines)
-    write_ucd_lines (tria, cell_index, out);
-
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
-
-  AssertThrow (out, ExcIO());
-}
-
-
-#if deal_II_dimension != 2
-
-template <int dim>
-void GridOut::write_xfig (
-  const Triangulation<dim>&,
-  std::ostream&,
-  const Mapping<dim>*) const
-{
-  Assert (false, ExcNotImplemented());
-}
-
-#else
-
-//TODO:[GK] Obey parameters
-//TODO:[GK] Flip y-axis?
-template <int dim>
-void GridOut::write_xfig (
-  const Triangulation<dim>& tria,
-  std::ostream&             out,
-  const Mapping<dim>*       /*mapping*/) const
-{
-  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)
-       {
-         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;
-           }
-         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 ();
-
-  AssertThrow (out, ExcIO());
-}
-
-#endif
-
-
-
-#if deal_II_dimension == 1
-
-unsigned int GridOut::n_boundary_faces (const Triangulation<1> &) const
-{
-  return 0;
-}
-
-unsigned int GridOut::n_boundary_lines (const Triangulation<1> &) const
-{
-  return 0;
-}
-
-
-unsigned int GridOut::n_boundary_faces (const Triangulation<1,2> &) const
-{
-  return 0;
-}
-
-unsigned int GridOut::n_boundary_lines (const Triangulation<1,2> &) const
-{
-  return 0;
-}
-#endif
-
-
-#if deal_II_dimension == 2
-
-unsigned int GridOut::n_boundary_lines (const Triangulation<2> &) const
-{
-  return 0;
-}
-
-unsigned int GridOut::n_boundary_lines (const Triangulation<2,3> &) const
-{
-  return 0;
-}
-#endif
-
-
-
-template <int dim, int spacedim>
-unsigned int GridOut::n_boundary_faces (const Triangulation<dim,spacedim> &tria) const
-{
-    typename Triangulation<dim,spacedim>::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++;
-
-  return n_faces;
-}
-
-
-
-template <int dim, int spacedim>
-unsigned int GridOut::n_boundary_lines (const Triangulation<dim, spacedim> &tria) const
-{
-    typename Triangulation<dim, spacedim>::active_line_iterator edge, endedge;
-  unsigned int n_lines = 0;
-
-  for (edge=tria.begin_active_line(), endedge=tria.end_line();
-       edge != endedge; ++edge)
-    if ((edge->at_boundary()) &&
-       (edge->boundary_indicator() != 0))
-      n_lines++;
-
-  return n_lines;
-}
-
-
-
-#if deal_II_dimension == 1
-
-void GridOut::write_msh_faces (const Triangulation<1> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-
-void GridOut::write_msh_faces (const Triangulation<1,2> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-
-void GridOut::write_msh_lines (const Triangulation<1> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-void GridOut::write_msh_lines (const Triangulation<1,2> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-#endif
-
-
-#if deal_II_dimension == 2
-
-void GridOut::write_msh_lines (const Triangulation<2> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-void GridOut::write_msh_lines (const Triangulation<2,3> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-#endif
-
-
-
-template <int dim, int spacedim>
-void GridOut::write_msh_faces (const Triangulation<dim,spacedim> &tria,
-                              const unsigned int        starting_index,
-                              std::ostream             &out) const
-{
-    typename Triangulation<dim,spacedim>::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())
-           << ' '
-           << static_cast<unsigned int>(face->boundary_indicator())
-           << ' ' << 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;
-      };
-}
-
-
-template <int dim, int spacedim>
-void GridOut::write_msh_lines (const Triangulation<dim, spacedim> &tria,
-                              const unsigned int        starting_index,
-                              std::ostream             &out) const
-{
-    typename Triangulation<dim,spacedim>::active_line_iterator line, endl;
-  unsigned int index=starting_index;
-
-  for (line=tria.begin_active_line(), endl=tria.end_line();
-       line != endl; ++line)
-    if (line->at_boundary() &&
-       (line->boundary_indicator() != 0))
-      {
-       out << index << " 1 ";
-       out << static_cast<unsigned int>(line->boundary_indicator())
-           << ' '
-           << static_cast<unsigned int>(line->boundary_indicator())
-           << " 2 ";
-                                        // note: vertex numbers are 1-base
-       for (unsigned int vertex=0; vertex<2; ++vertex)
-         out << ' '
-             << line->vertex_index(GeometryInfo<dim-2>::ucd_to_deal[vertex])+1;
-       out << '\n';
-
-       ++index;
-      };
-}
-
-
-
-#if deal_II_dimension == 1
-
-void GridOut::write_ucd_faces (const Triangulation<1> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-void GridOut::write_ucd_faces (const Triangulation<1,2> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-void GridOut::write_ucd_lines (const Triangulation<1> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-void GridOut::write_ucd_lines (const Triangulation<1,2> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-#endif
-
-
-
-#if deal_II_dimension == 2
-
-void GridOut::write_ucd_lines (const Triangulation<2> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-void GridOut::write_ucd_lines (const Triangulation<2,3> &,
-                              const unsigned int,
-                              std::ostream &) const
-{
-  return;
-}
-
-#endif
-
-
-
-template <int dim, int spacedim>
-void GridOut::write_ucd_faces (const Triangulation<dim, spacedim> &tria,
-                              const unsigned int        starting_index,
-                              std::ostream             &out) const
-{
-    typename Triangulation<dim,spacedim>::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';
-
-       ++index;
-      };
-}
-
-
-
-template <int dim, int spacedim>
-void GridOut::write_ucd_lines (const Triangulation<dim, spacedim> &tria,
-                              const unsigned int        starting_index,
-                              std::ostream             &out) const
-{
-    typename Triangulation<dim, spacedim>::active_line_iterator line, endl;
-  unsigned int index=starting_index;
-
-  for (line=tria.begin_active_line(), endl=tria.end_line();
-       line != endl; ++line)
-    if (line->at_boundary() &&
-       (line->boundary_indicator() != 0))
-      {
-       out << index << "  "
-           << static_cast<unsigned int>(line->boundary_indicator())
-           << "  line    ";
-       // note: vertex numbers are 1-base
-       for (unsigned int vertex=0; vertex<2; ++vertex)
-         out << line->vertex_index(GeometryInfo<dim-2>::ucd_to_deal[vertex])+1 << ' ';
-       out << '\n';
-
-       ++index;
-      };
-}
-
-
-#if deal_II_dimension==1
-
-template <int dim, int spacedim>
-void GridOut::write_gnuplot (
-  const Triangulation<dim,spacedim> &tria,
-  std::ostream             &out,
-  const Mapping<dim,spacedim>       *) const
-{
-  AssertThrow (out, ExcIO());
-
-  typename Triangulation<dim,spacedim>::active_cell_iterator        cell=tria.begin_active();
-  const typename Triangulation<dim,spacedim>::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';
-    }
-
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
-
-  AssertThrow (out, ExcIO());
-}
-
-#endif
-
-
-#if deal_II_dimension==2
-
-template <int dim, int spacedim>
-void GridOut::write_gnuplot (
-  const Triangulation<dim,spacedim> &tria,
-  std::ostream           &out,
-  const Mapping<dim,spacedim>       *mapping) const
-{
-  AssertThrow (out, ExcIO());
-
-  const unsigned int n_additional_points=
-    gnuplot_flags.n_boundary_face_points;
-  const unsigned int n_points=2+n_additional_points;
-
-  typename Triangulation<dim,spacedim>::active_cell_iterator        cell=tria.begin_active();
-  const typename Triangulation<dim,spacedim>::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<dim-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<dim-1> quadrature(boundary_points, dummy_weights);
-
-      q_projector = new Quadrature<dim> (QProjector<dim>::project_to_all_faces(quadrature));
-    }
-
-  for (; cell!=endc; ++cell)
-    {
-      if (gnuplot_flags.write_cell_numbers)
-       out << "# cell " << cell << '\n';
-
-      if (mapping==0 ||
-         (!cell->at_boundary() && !gnuplot_flags.curved_inner_cells))
-       {
-                                          // 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'
-             << '\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
-       {
-         for (unsigned int face_no=0;
-              face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-           {
-             const typename Triangulation<dim,spacedim>::face_iterator
-               face = cell->face(face_no);
-             if (face->at_boundary() || gnuplot_flags.curved_inner_cells)
-               {
-                                                  // 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';
-               }
-             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 (q_projector != 0)
-    delete q_projector;
-
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
-
-  AssertThrow (out, ExcIO());
-}
-
-#endif
-
-
-#if deal_II_dimension==3
-
-template <int dim, int spacedim>
-void GridOut::write_gnuplot (const Triangulation<dim,spacedim> &tria,
-                            std::ostream           &out,
-                            const Mapping<dim,spacedim>       *mapping) const
-{
-  AssertThrow (out, ExcIO());
-
-  const unsigned int n_additional_points=
-    gnuplot_flags.n_boundary_face_points;
-  const unsigned int n_points=2+n_additional_points;
-
-  typename Triangulation<dim,spacedim>::active_cell_iterator        cell=tria.begin_active();
-  const typename Triangulation<dim,spacedim>::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));
-    }
-
-  for (; cell!=endc; ++cell)
-    {
-      if (gnuplot_flags.write_cell_numbers)
-       out << "# cell " << cell << '\n';
-
-      if (mapping==0 || n_points==2 ||
-         (!cell->has_boundary_lines() && !gnuplot_flags.curved_inner_cells))
-       {
-                                          // 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
-       {
-         for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
-           {
-             const typename Triangulation<dim,spacedim>::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<spacedim> 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
-               {
-                 for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
-                   {
-                     const typename Triangulation<dim,spacedim>::line_iterator
-                       line=face->line(l);
-
-                     const Point<spacedim> &v0=line->vertex(0),
-                                           &v1=line->vertex(1);
-                     if (line->at_boundary() || gnuplot_flags.curved_inner_cells)
-                       {
-                                                          // transform_real_to_unit_cell
-                                                          // could be
-                                                          // replaced
-                                                          // by using
-                                                          // QProjector<dim>::project_to_line
-                                                          // which is
-                                                          // not yet
-                                                          // implemented
-                         const Point<spacedim> 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';
-
-                     out << '\n' << '\n';
-                   }
-               }
-           }
-       }
-    }
-
-  if (q_projector != 0)
-    delete q_projector;
-
-
-                                  // make sure everything now gets to
-                                  // disk
-  out.flush ();
-
-  AssertThrow (out, ExcIO());
-}
-
-#endif
-
-struct LineEntry
-{
-    Point<2> first;
-    Point<2> second;
-    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)
-      {}
-};
-
-
-#if deal_II_dimension==1
-
-template <int dim>
-void GridOut::write_eps (
-  const Triangulation<dim> &,
-  std::ostream &,
-  const Mapping<dim> *) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-#else
-
-template <int dim>
-void GridOut::write_eps (
-  const Triangulation<dim> &tria,
-  std::ostream             &out,
-  const Mapping<dim>       *mapping) const
-{
-
-  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<const GridOutFlags::EpsFlagsBase&>(eps_flags_2) :
-                      (dim==3 ?
-                       static_cast<const GridOutFlags::EpsFlagsBase&>(eps_flags_3) :
-                       *static_cast<const 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_cell_iterator
-         cell=tria.begin_active(),
-         endc=tria.end();
-
-       for (; cell!=endc; ++cell)
-         for (unsigned int line_no=0;
-              line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
-           {
-             typename Triangulation<dim>::line_iterator
-               line=cell->line(line_no);
-
-                                              // 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
-             if (!line->has_children() &&
-                 (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(),
-                                              cell->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)
-               {
-                 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(),
-                                                        cell->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(),
-                                                    cell->level()));
-                   };
-               };
-         };
-
-       break;
-      };
-
-      case 3:
-      {
-                                        // curved boundary output
-                                        // presently not supported
-       Assert (mapping == 0, ExcNotImplemented());
-
-       typename Triangulation<dim>::active_cell_iterator
-         cell=tria.begin_active(),
-         endc=tria.end();
-
-                                        // 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 = 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 (; cell!=endc; ++cell)
-         for (unsigned int line_no=0;
-              line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
-           {
-             typename Triangulation<dim>::line_iterator
-               line=cell->line(line_no);
-             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(),
-                                            cell->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;
-
-  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_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_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));
-
-
-                                  // 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))
-       {
-         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';
-    };
-
-                                  // 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 << (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))
-    {
-      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());
-}
-
-#endif
-
-
-template <int dim, int spacedim>
-void GridOut::write (
-    const Triangulation<dim, spacedim> &tria,
-  std::ostream             &out,
-  const OutputFormat        output_format,
-  const Mapping<dim,spacedim>       *mapping) const
-{
-  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 msh:
-           write_msh (tria, out);
-           return;
-    }
-
-  Assert (false, ExcInternalError());
-}
-
-
-template <int dim, int spacedim>
-void GridOut::write (
-    const Triangulation<dim, spacedim> &tria,
-  std::ostream             &out,
-  const Mapping<dim,spacedim>       *mapping) const
-{
-  write(tria, out, default_format, mapping);
-}
-
-
-//TODO: Should all write functions be instantiated?
-
-// explicit instantiations
-template void GridOut::write_dx
-(const Triangulation<deal_II_dimension>&,
- std::ostream&) const;
-template void GridOut::write_msh
-(const Triangulation<deal_II_dimension>&,
- std::ostream&) const;
-template void GridOut::write_xfig
-(const Triangulation<deal_II_dimension>&,
- std::ostream&,
- const Mapping<deal_II_dimension,deal_II_dimension>*) const;
-template void GridOut::write_gnuplot
-(const Triangulation<deal_II_dimension>&,
- std::ostream&,
- const Mapping<deal_II_dimension,deal_II_dimension>*) const;
-template void GridOut::write_ucd<deal_II_dimension>
-(const Triangulation<deal_II_dimension> &,
- std::ostream &) const;
-template void GridOut::write_eps<deal_II_dimension>
-(const Triangulation<deal_II_dimension> &,
- std::ostream &,
- const Mapping<deal_II_dimension,deal_II_dimension> *) const;
-template void GridOut::write<deal_II_dimension>
-(const Triangulation<deal_II_dimension> &,
- std::ostream &, const OutputFormat,
- const Mapping<deal_II_dimension,deal_II_dimension> *) const;
-template void GridOut::write<deal_II_dimension>
-(const Triangulation<deal_II_dimension> &,
- std::ostream &, const Mapping<deal_II_dimension,deal_II_dimension> *) const;
-
-#if deal_II_dimension == 1 || deal_II_dimension ==2
-// explicit instantiations for codimension one
-template void GridOut::write_msh
-(const Triangulation<deal_II_dimension, deal_II_dimension+1>&,
- std::ostream&) const;
-template void GridOut::write_ucd
-(const Triangulation<deal_II_dimension, deal_II_dimension+1> &,
- std::ostream &) const;
-template void GridOut::write_gnuplot
-(const Triangulation<deal_II_dimension,deal_II_dimension+1>&,
- std::ostream&,
- const Mapping<deal_II_dimension,deal_II_dimension+1>*) const;
-#endif
-
-DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/source/grid/grid_out.inst.in b/deal.II/source/grid/grid_out.inst.in
new file mode 100644 (file)
index 0000000..70af179
--- /dev/null
@@ -0,0 +1,60 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS)
+  {
+    template void GridOut::write_dx
+      (const Triangulation<deal_II_dimension>&,
+       std::ostream&) const;
+    template void GridOut::write_msh
+      (const Triangulation<deal_II_dimension>&,
+       std::ostream&) const;
+    template void GridOut::write_xfig
+      (const Triangulation<deal_II_dimension>&,
+       std::ostream&,
+       const Mapping<deal_II_dimension,deal_II_dimension>*) const;
+    template void GridOut::write_gnuplot
+      (const Triangulation<deal_II_dimension>&,
+       std::ostream&,
+       const Mapping<deal_II_dimension,deal_II_dimension>*) const;
+    template void GridOut::write_ucd<deal_II_dimension>
+      (const Triangulation<deal_II_dimension> &,
+       std::ostream &) const;
+    template void GridOut::write_eps<deal_II_dimension>
+      (const Triangulation<deal_II_dimension> &,
+       std::ostream &,
+       const Mapping<deal_II_dimension,deal_II_dimension> *) const;
+    template void GridOut::write<deal_II_dimension>
+      (const Triangulation<deal_II_dimension> &,
+       std::ostream &, const OutputFormat,
+       const Mapping<deal_II_dimension,deal_II_dimension> *) const;
+    template void GridOut::write<deal_II_dimension>
+      (const Triangulation<deal_II_dimension> &,
+       std::ostream &, const Mapping<deal_II_dimension,deal_II_dimension> *) const;
+
+#if deal_II_dimension == 1 || deal_II_dimension ==2
+// explicit instantiations for codimension one
+    template void GridOut::write_msh
+      (const Triangulation<deal_II_dimension, deal_II_dimension+1>&,
+       std::ostream&) const;
+    template void GridOut::write_ucd
+      (const Triangulation<deal_II_dimension, deal_II_dimension+1> &,
+       std::ostream &) const;
+    template void GridOut::write_gnuplot
+      (const Triangulation<deal_II_dimension,deal_II_dimension+1>&,
+       std::ostream&,
+       const Mapping<deal_II_dimension,deal_II_dimension+1>*) const;
+#endif
+  }
+

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.