* Declaration of a name for each of the
* different output formats.
*/
- enum OutputFormat { gnuplot, eps };
+ enum OutputFormat { gnuplot, eps, ucd };
/**
* Write the triangulation in the
* of what happens here.
*/
template <int dim>
- static void write_gnuplot (const Triangulation<dim> &tria,
- ostream &out);
+ void write_gnuplot (const Triangulation<dim> &tria,
+ ostream &out);
+
+ /**
+ * Write the triangulation in the
+ * ucd format. See the general
+ * documentation for a description
+ * of what happens here.
+ */
+ template <int dim>
+ void write_ucd (const Triangulation<dim> &tria,
+ ostream &out);
/**
* Write the triangulation in the
* of what happens here.
*/
template <int dim>
- static void write_eps (const Triangulation<dim> &tria,
- ostream &out);
-
+ void write_eps (const Triangulation<dim> &tria,
+ ostream &out);
+
/**
* Write data and grid to #out# according
* to the given data format. This function
* #write_*# function.
*/
template <int dim>
- static void write (const Triangulation<dim> &tria,
- ostream &out,
- const OutputFormat output_format);
+ void write (const Triangulation<dim> &tria,
+ ostream &out,
+ const OutputFormat output_format);
/**
* Provide a function which tells us which
* formats are defined:
* \begin{itemize}
* \item #gnuplot#: #.gnuplot#
+ * \item #ucd#: #.inp#
* \item #eps#: #.eps#.
* \end{itemize}
*
* Exception
*/
DeclException0 (ExcIO);
+
+ private:
+
+ /**
+ * Write the grid information about
+ * faces to #out#. Only those faces
+ * are printed which are on the boundary
+ * and which have a boundary indicator
+ * not equal to zero, since the latter
+ * is the default for boundary faces.
+ *
+ * Since cells and faces are continuously
+ * numbered, the #starting_index# for
+ * the numbering of the faces is passed
+ * also.
+ *
+ * This function unfortunately can not
+ * be included in the regular #write_ucd#
+ * function, since it needs special
+ * treatment for the case #dim==1#, in
+ * which case the face iterators are
+ * #void*#'s and lack the member functions
+ * which are called. We would not actually
+ * call these functions, but the compiler
+ * would complain anyway when compiling
+ * the function for #dim==1#. Bad luck.
+ */
+ template <int dim>
+ void write_ucd_faces (const Triangulation<dim> &tria,
+ const unsigned int starting_index,
+ ostream &out) const;
+
+ /**
+ * Return the number of faces in the
+ * triangulation which have a boundary
+ * indicator not equal to zero. Only
+ * these faces are explicitely printed
+ * in the #write_*# functions;
+ * all faces with indicator 255 are
+ * interior ones and an indicator with
+ * value zero for faces at the boundary
+ * are considered default.
+ *
+ * This function always returns an empty
+ * list in one dimension.
+ *
+ * The reason for this function is the
+ * same as for #write_ucd_faces#. See
+ * there for more information.
+ */
+ template <int dim>
+ unsigned int n_boundary_faces (const Triangulation<dim> &tria) const;
};
+template <int dim>
+void GridOut::write_ucd (const Triangulation<dim> &tria,
+ ostream &out)
+{
+ AssertThrow (out, ExcIO());
+
+ typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
+ const typename Triangulation<dim>::active_cell_iterator endc=tria.end();
+
+ // first loop over all cells to
+ // find out the vertices which
+ // are in use. copy them for fast
+ // output
+ vector<Point<dim> > vertices (tria.n_vertices());
+ vector<bool> vertex_used (tria.n_vertices(), false);
+ unsigned int n_vertices = 0;
+
+ for (cell=tria.begin_active(); cell!=endc; ++cell)
+ for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex)
+ // if not yet copied
+ if (vertex_used[cell->vertex_index(vertex)] == false)
+ {
+ vertex_used[cell->vertex_index(vertex)] = true;
+ vertices[cell->vertex_index(vertex)] = cell->vertex(vertex);
+ ++n_vertices;
+ };
+
+
+ // write preamble
+ if (true)
+ {
+/*
+ // block this to have local
+ // variables destroyed after
+ // use
+ time_t time1= time (0);
+ tm *time = localtime(&time1);
+ out << "# This file was generated by the deal.II library." << endl
+ << "# Date = "
+ << time->tm_year+1900 << "/"
+ << time->tm_mon+1 << "/"
+ << time->tm_mday << endl
+ << "# Time = "
+ << time->tm_hour << ":"
+ << setw(2) << time->tm_min << ":"
+ << setw(2) << time->tm_sec << endl
+ << "#" << endl
+ << "# For a description of the UCD format see the AVS Developer's guide."
+ << endl
+ << "#" << endl;
+*/
+ };
+
+ // start with ucd data
+ out << n_vertices << ' '
+ << tria.n_active_cells() + n_boundary_faces(tria)
+ << " 0 0 0" // no data
+ << 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=dim+1; d<=3; ++d)
+ out << " 0"; // fill with zeroes
+ out << endl;
+ };
+
+ // 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(vertex)+1 << ' ';
+ out << endl;
+ };
+
+ // write faces with non-zero boundary
+ // indicator
+ write_ucd_faces (tria, cell_index, out);
+
+ AssertThrow (out, ExcIO());
+};
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+unsigned int GridOut::n_boundary_faces (const Triangulation<1> &) const
+{
+ return 0;
+};
+
+#endif
+
+
+
+template <int dim>
+unsigned int GridOut::n_boundary_faces (const Triangulation<dim> &tria) const
+{
+ typename Triangulation<dim>::active_face_iterator face, endf;
+ unsigned long 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;
+};
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+void GridOut::write_ucd_faces (const Triangulation<1> &,
+ const unsigned int,
+ ostream &) const
+{
+ return;
+};
+
+#endif
+
+
+
+template <int dim>
+void GridOut::write_ucd_faces (const Triangulation<dim> &tria,
+ const unsigned int starting_index,
+ ostream &out) const
+{
+ typename Triangulation<dim>::active_face_iterator face, endf;
+ unsigned int index=starting_index;
+
+ for (face=tria.begin_active_face(), endf=tria.end_face();
+ face != endf; ++face)
+ if (face->at_boundary() &&
+ (face->boundary_indicator() != 0))
+ {
+ out << index << " "
+ << static_cast<unsigned int>(face->boundary_indicator())
+ << " ";
+ switch (dim)
+ {
+ case 2: out << "line "; break;
+ case 3: out << "quad "; break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
+ out << face->vertex_index(vertex) << ' ';
+ out << endl;
+
+ ++index;
+ };
+};
+
+
+
template <int dim>
void GridOut::write_gnuplot (const Triangulation<dim> &tria,
ostream &out)
-// explicit instantiations
-template void GridOut::write_gnuplot (const Triangulation<deal_II_dimension> &, ostream &);
-template void GridOut::write_eps (const Triangulation<deal_II_dimension> &, ostream &);
+template <int dim>
+void GridOut::write (const Triangulation<dim> &tria,
+ ostream &out,
+ OutputFormat output_format)
+{
+ switch (output_format)
+ {
+ case ucd:
+ write_ucd (tria, out);
+ return;
+
+ case gnuplot:
+ write_gnuplot (tria, out);
+ return;
+
+ case eps:
+ write_eps (tria, out);
+ return;
+ };
+
+ Assert (false, ExcInternalError());
+};
+
+
+
+
+// explicit instantiations. note that write instantiates all the other
+// functions as needed
template void GridOut::write (const Triangulation<deal_II_dimension> &, ostream &, OutputFormat);
projects / dealii-svn.git / commitdiff