+++ /dev/null
-//---------------------------- data_io.h ---------------------------
-// $Id$
-// Version: $Name$
-//
-// Copyright (C) 1998, 1999, 2000 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.
-//
-//---------------------------- data_io.h ---------------------------
-#ifndef __deal2__data_io_h
-#define __deal2__data_io_h
-
-
-#include <base/exceptions.h>
-#include <vector>
-#include <string>
-
-template <typename number> class Vector;
-template <int dim> class DoFHandler;
-class EpsOutputData;
-
-/**
- * This class is deprecated. Use the @p{DataOut} class instead.
- *
- * This class implements an output mechanism for grid and simulation data
- * in several formats.
- * At present it supports output in UCD (unstructured cell data) and
- * GNUPLOT format. Partly supported are POVRAY and encapsulated postscript.
- * POVRAY allows for only one data set and does not support cell data.
- * Encapsulated Postscript supports Q1-Elements only.
- *
- * It allows the user to attach a degree of freedom handler object
- * (@ref{DoFHandler}) which also gives access to the geometry data of the
- * underlying triangulation and to add data vectors of which the values
- * are to be written.
- *
- *
- * @sect3{Limitations}
- *
- * At present, no grouping of components to vectors is implemented, i.e.
- * you can only write each component independent of the others. Also, it
- * is not possible to output calculations which were performed on elements
- * with more or less than one degree of freedom per vertex.
- *
- *
- * @sect3{UCD format}
- *
- * The UCD format is described in the AVS developer's guide. Due to
- * limitations in the present format, only node based data can be output,
- * so higher order elements are only written with their node values, no
- * interior or line values are used. No use is made of the possibility
- * to give model data since these are not supported by all
- * UCD aware programs. You may give cell data, but that is not supported by
- * all programs.
- *
- * The ASCII UCD format is used. In future versions, a binary version may
- * follow up.
- *
- * Note that to enumerate the vertices, not the vertex index is used but
- * the index of the degree of freedom located on this vertex. This makes
- * the mapping between the vertices and the entries in the data vectors
- * much easier.
- *
- *
- * @sect3{GNUPLOT draft format}
- *
- * The GNUPLOT format is not able to handle data on unstructured grids
- * directly. Directly would mean that you only give the vertices and
- * the solution values thereon and the program constructs its own grid
- * to represent the data. This is only possible for a structured tensor
- * product grid in two dimensions.
- *
- * In one dimension, the format is obviously @p{x v1 v2 ...}, where @p{x}
- * is the coordinate value of a grid point, while the @p{vi} are the
- * vector elements referring to the present node. Within GNUPLOT,
- * call @p{plot "filename" using 1:x}. @p{x} denotes the number of the data set you
- * want to see plus one. For example @p{using 1:4} would mean to plot the
- * third data vector.
- *
- * For more than one dimension, the @p{DataOut_Old<dim>::write_gnuplot()} somehow
- * duplicates the functionality of the @p{Triangulation<dim>::print_gnuplot()}
- * functions. These, however, offer more functionality in some respect.
- * The grid is represented as a sequence of lines, where each cell is
- * a sequence of five vertices (the first one is appended at the end to close
- * the contour of the cell) with the data appended after each vertex. Each cell
- * is therefore a sequence of five lines @p{x y v1 v2 ...} forming together
- * the bounding line of this cell. After each cell, two newlines are inserted
- * to prevent GNUPLOT from joining the lines bounding two cells.
- *
- * To view the results in two dimensions, use @p{set data style lines}
- * within gnuplot and call @p{plot "filename"} to see the grid. Use
- * @p{set parametric} and @p{splot "filename" using 1:2:x} to get a 3d surface
- * plot of the (@p{x-2})th data set. For example, using @p{x=4} would mean to
- * plot the second data set.
- *
- * This format is somewhat restricted to continuous data and to finite elements
- * of first order only. The reason for the first restriction is that it takes
- * nodal values and can therefore only work if a finite element has degrees of
- * freedom in the vertices of each cell. This is not the case for discontinuous
- * elements. The second restriction is only a problem of quality of output: you
- * actually can print quadratic or higher order elements using this style, but
- * you will only seem the contour of each cell with the bounding lines of each
- * cell being straight lines. You won't see the structure of the solution in the
- * interior of a cell nor on the lines bounding it.
- *
- *
- * @sect3{GNUPLOT 'quality' format}
- *
- * To remedy the abovementioned problems, the GNUPLOT 'quality' format was
- * introduced. As noted above, GNUPLOT can handle tensor grids and displays
- * them quite nicely, including hidden line removal. It can also handle more
- * than one tensor grid, then called a patch. The idea now is to describe
- * each cell as a patch of $N\times N$ points (in 2D). For linear elements,
- * $2\times 2$ is sufficient, but for higher elements you will want to use
- * a higher number. The @p{write_gnuplot} function writes the data in this format
- * and the argument it takes is the number of subintervals it divides each cell
- * into, i.e. $N-1$.
- *
- * This output routine also addresses the problem introduced with discontinuous
- * elements, since it takes its data locally from each cell and displays it
- * as a patch separately; it therefore does not use continuity and GNUPLOT will
- * also plot discontinuities, if they are there.
- *
- *
- * @sect3{POVRAY mesh format}
- *
- * POVRAY is a ray tracing tool in the public domain and poduces high quality
- * images of three dimensional scenes. Since this tool can handle only one
- * scene at a time, only the first data vector is output, but this may change
- * in future versions of this library. For a reference of the data format
- * and the usage of POVRAY see the extensive manual of that program.
- *
- * The POVRAY output format is presently only supported for two spatial
- * dimensions. It displays every quadrilateral as two triangles, which clearly
- * is not an optimal solution, but necessary since POVRAY presently does not
- * support bilinear quadrilaterals and using polygons as vastly suboptimal in
- * term of memory and speed to the triangle mesh supported by POVRAY.
- *
- *
- * @sect3{Encapsulated Postscript format}
- *
- * There is a function for generating encapsulated Postscript
- * (EPS) without the need for another graphics tool. This functionality
- * is provided through @p{write_eps} and @p{EpsOutputData}.
- * @p{write_eps} can use one data vector for height information and one
- * cell vector for shading information. Control is done through an
- * object of class @p{EpsOutputData} as follows or by default.
- *
- * Vectors are added as usual by @p{add_data_vector}. Then one has to
- * decide, wether to produce a 2D or 3D plot. This is done by setting
- * @p{height_info} to
- * @begin{itemize}
- * @item @em{NoHeight} for 2D-Output (or Top-View thats the same by no
- * turning is done) or to
- * @item @em{HeightVector} for 3D-Output. You have to attach a
- * @p{dof_data_vector} to actually get 3D. If you don't then output
- * will be generated in 2D.
- * @item @em{DefaultHeight} is 3D if there is a @p{dof_data_vector} and 2D if
- * none is present.
- * @end{itemize}
- * For 3D-Output one has to set @p{azimuth} and @p{elevation} for the
- * angle of view and @p{height_vector} to the number of the @p{dof_data}
- * vector that provides the height information to be used. The default
- * values are analogous to GNUPlot for azimuth and elevation and
- * vector 0 for the height vector.
- *
- * The cells can be shaded in four different modes, controlled by the
- * attribute @p{cell_shading}:
- * @begin{enumerate}
- * @item @em{NoShading} provides transparent shading.
- * @item @em{ShadingVector} uses a cell vector to do shading. The number
- * of the cell vector to be uses is provided in @p{cell_vector}. To
- * scale the cell vector there is the method @p{color}. It is called
- * with the actual value of the cell, the maximum and the minimum
- * value of a cell in the cell vector. It returns three values for
- * red, green and blue. If there no @p{cell_data} vector than there is
- * transparent shading.
- * @item @em{LightShaded} just shades the plot. This is controlled by
- * the vector @p{light} which stores the direction of the light
- * beams. This is done only if there is height information.
- * @item @em{DefaultShading} is controlled by presence of different
- * vectors. If there no height information then do no
- * shading. Otherwise if there is @p{cell_data} use this for shading.
- * Otherwise do light shading.
- * @end{enumerate}
- *
- * Finnaly one can choose to mark the cell boundaries by setting
- * @p{cell_boundary_shading}. It can take one of four values:
- * @begin{itemize}
- * @item NoBoundary for no cell boundaries,
- * @item DefaultBoundary or
- * @item BlackBoundary for black cell boundaries,
- * @item WhiteBoundary for white cell boundaries,
- * @end{itemize}
- *
- * Another interesting feature is that you can write multiple
- * eps-pictures to one file by just doing several invocations of
- * @p{write_eps}. One than can switch between the different graphics
- * with the @p{>>} Button in GhostView for example.
- *
- *
- * @sect3{GMV format}
- *
- * The @p{write_gmv} function and the @p{write} function through the @p{gmv} parameter
- * write the data in a format understood by the GMV (general mesh viewer)
- * program. This program is able to generate 2d and 3d plots of almost
- * arbitrarily many data sets, along with shading, cuts through data sets and
- * many other nifty features.
- *
- * Data is written in the following format: nodes are considered the vertices
- * of the triangulation. In spatial dimensions less than three, zeroes are
- * inserted for the missing coordinates. The data vectors are written as
- * node or cell data, where for the first the data space is interpolated to
- * (bi-,tri-)linear elements.
- *
- * @author Wolfgang Bangerth, Guido Kanschat, Stefan Nauber, 1998, 1999
- */
-template <int dim>
-class DataOut_Old {
- public:
- /**
- * Provide a data type specifying the
- * presently supported output formats.
- */
- enum OutputFormat { ucd, gnuplot, gnuplot_draft, povray_mesh, eps, gmv };
-
- /**
- * Constructor
- */
- DataOut_Old ();
-
- /**
- * Designate a dof handler to be used
- * to extract geometry data and the
- * mapping between nodes and node values.
- */
- void attach_dof_handler (const DoFHandler<dim> &);
-
- /**
- * Add a data vector together with its
- * name and the physical unit
- * (for example meter, kelvin, etc). By
- * default, "<dimensionless>" is assumed
- * for the units.
- *
- * A pointer to the vector is stored, so
- * you have to make sure the vector
- * exists at that address at least as
- * long as you call the
- * @p{write_*} functions.
- *
- * It is assumed that the vector has the
- * same number of components as there are
- * degrees of freedom in the dof handler,
- * in which case it is assumed to be a
- * vector storing nodal data; or the size
- * may be the number of active cells on
- * the present grid, in which case it is
- * assumed to be a cell data vector.
- *
- * The name and unit of a data vector shall
- * only contain characters which are
- * letters, underscore and a few other
- * ones. Refer to the @p{ExcInvalidCharacter}
- * exception declared in this class to
- * see which characters are valid and which
- * are not.
- */
- void add_data_vector (const Vector<double> &data,
- const string &name,
- const string &units="<dimensionless>");
-
- /**
- * Release the pointers to the data
- * vectors. You have to set all data
- * entries again using the
- * @p{add_data_vector} function. The pointer
- * to the dof handler remains stored,
- * however.
- */
- void clear_data_vectors ();
-
- /**
- * Write the stored data to the given
- * stream in UCD data. You may have
- * written any comment to that stream
- * before calling this function. Comments
- * start with the \# character in the
- * first column of a line and may only
- * appear at the beginning of a file,
- * without non-comment lines inbetween.
- */
- void write_ucd (ostream &out) const;
-
- /**
- * Write data and grid in GNUPLOT format.
- * Write a patch for each grid cell using
- * @p{accuracy} subintervals per dimension.
- */
- void write_gnuplot (ostream &out, const unsigned int accuracy=1) const;
-
- /**
- * Write data and grid in GNUPLOT format.
- * Only the edges between cells are written.
- */
- void write_gnuplot_draft (ostream &out) const;
-
- /**
- * Write data of first vector and grid
- * in POVRAY format. Further data vectors
- * as well as cell data is ignored.
- */
- void write_povray_mesh (ostream &out) const;
-
- /**
- * Write data of first vector and grid in
- * Encapsulated postscript format.
- * Further data vectors
- * as well as cell data is ignored.
- *
- * If no special output data is given, a
- * default constructed object is used.
- */
- void write_eps (ostream &out,
- const EpsOutputData &eod = EpsOutputData()) const;
-
- /**
- * Write data and grid in GMV format.
- */
- void write_gmv (ostream &out) const;
-
- /**
- * Write data and grid to @p{out} according
- * to the given data format. This function
- * simply calls the appropriate
- * @p{write_*} function.
- *
- * In case of gnuplot output, the
- * standard accuracy of 1 is
- * chosen.
- */
- void write (ostream &out, const OutputFormat output_format) const;
-
- /**
- * Provide a function which tells us which
- * suffix with a given output format
- * usually has. At present the following
- * formats are defined:
- * @begin{itemize}
- * @item @p{ucd}: @p{.inp}
- * @item @p{gnuplot} and @p{gnuplot_draft}:
- * @p{.gnuplot}
- * @item @p{povray_mesh}: @p{.pov}
- * @item @p{eps}: @p{.eps}
- * @item @p{gmv}: @p{.gmv}.
- * @end{itemize}
- *
- * Since this function does not need data
- * from this object, it is static and can
- * thus be called without creating an
- * object of this class.
- */
- static string default_suffix (const OutputFormat output_format);
-
- /**
- * Return the @p{OutputFormat} value
- * corresponding to the given string. If
- * the string does not match any known
- * format, an exception is thrown.
- *
- * Since this function does not need data
- * from this object, it is static and can
- * thus be called without creating an
- * object of this class. Its main purpose
- * is to allow a program to use any
- * implemented output format without the
- * need to extend the program's parser
- * each time a new format is implemented.
- *
- * To get a list of presently available
- * format names, e.g. to give it to the
- * @p{ParameterHandler} class, use the
- * function @p{get_output_format_names ()}.
- */
- static OutputFormat parse_output_format (const string &format_name);
-
- /**
- * Return a list of implemented output
- * formats. The different names are
- * separated by vertical bar signs (@p{`|'})
- * as used by the @p{ParameterHandler}
- * classes.
- */
- static string get_output_format_names ();
-
- /**
- * Exception
- */
- DeclException0 (ExcIncorrectDofsPerVertex);
- /**
- * Exception
- */
- DeclException0 (ExcNoDoFHandlerSelected);
- /**
- * Exception
- */
- DeclException3 (ExcInvalidVectorSize,
- int, int, int,
- << "The vector has size " << arg1
- << " but the DoFHandler objects says there are " << arg2
- << " degrees of freedom and there are " << arg3
- << " active cells.");
- /**
- * Exception
- */
- DeclException1 (ExcInvalidCharacter,
- string,
- << "Please use only the characters [a-zA-Z0-9_<>()] for" << endl
- << "description strings since AVS will only accept these." << endl
- << "The string you gave was <" << arg1 << ">.");
- /**
- * Exception
- */
- DeclException0 (ExcIO);
- /**
- * Exception
- */
- DeclException0 (ExcInvalidState);
-
- private:
-
- /**
- * Declare an entry in the list of
- * data elements.
- */
- struct DataEntry {
- /**
- * Default constructor, constructs
- * invalid object, but is needed for
- * the STL classes.
- */
- DataEntry ();
-
- /**
- * Constructor
- */
- DataEntry (const Vector<double> *data, const string name, const string units);
-
- /**
- * Pointer to the data vector.
- */
- const Vector<double> *data;
-
- /**
- * Name of this component.
- */
- string name;
-
- /**
- * Physical unit name of this
- * component.
- */
- string units;
- };
-
- /**
- * Pointer to the dof handler object.
- */
- const DoFHandler<dim> *dofs;
-
- /**
- * List of data elements with vectors of
- * values for each degree of freedom.
- */
- vector<DataEntry> dof_data;
-
- /**
- * List of data elements with vectors of
- * values for each cell.
- */
- vector<DataEntry> cell_data;
-
-
-/**
- * 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 @p{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 @p{write_ucd_faces}. See
- * there for more information.
- */
- unsigned int n_boundary_faces () const;
-
- /**
- * Write the grid information about
- * faces to @p{out}. Only those faces
- * are printed which are on the boundary
- * and which have a boundary indicator
- * not equal to zero, since the latter
- * is the default for boundary faces.
- *
- * Since cells and faces are continuously
- * numbered, the @p{starting_index} for
- * the numbering of the faces is passed
- * also.
- *
- * This function unfortunately can not
- * be included in the regular @p{write_ucd}
- * function, since it needs special
- * treatment for the case @p{dim==1}, in
- * which case the face iterators are
- * @p{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 @p{dim==1}. Bad luck.
- */
- void write_ucd_faces (ostream &out,
- const unsigned int starting_index) const;
-
- /**
- * Class of data for vertices
- * for eps output.
- */
- class EpsVertexData{
- public:
- /**
- * Coordinates
- */
- double x;
- double y;
- double z;
-
- /**
- * Default Constructor;
- * needed for STL.
- */
- EpsVertexData()
- {};
-
- /**
- * Constructor that takes three
- * coordinate values as argument.
- */
- EpsVertexData(double a, double b, double c):x(a),y(b),z(c)
- {};
-
- /**
- * Transformation method
- */
- void turn(double azi, double ele);
- };
-
- /**
- * Class of cell data for output.
- * For eps output some calculations have
- * to be done between transformation and
- * projection and output. Because of this all
- * output data is put into a STL multiset.
- */
- class EpsCellData{
- public:
-
- /**
- * STL-vector of vertices
- */
- EpsVertexData vertices[4];
-
- /**
- * Color values
- */
- float red;
- float green;
- float blue;
-
- /**
- * Transformation method
- */
- void turn(double azi, double ele);
-
- /**
- * Comparison operator for
- * sorting
- */
- bool operator < (const EpsCellData &) const;
- };
-};
-
-
-/**
- * Structure for the control of encapsulated postscript output. See
- * general documentation of class @p{DataOut_Old} for description.
- *
- * @author Stefan Nauber
- */
-class EpsOutputData{
- public:
- /**
- * Types of height info
- */
- enum HeightInfo {
- DefaultHeight, NoHeight, HeightVector
- };
-
- /**
- * Types of cell shading
- */
- enum CellShading {
- DefaultShading, NoShading, ShadingVector, LightShaded
- };
-
- /**
- * Types of cell boundary shading
- */
- enum CellBoundaryShading {
- DefaultBoundary, NoBoundary, BlackBoundary, WhiteBoundary
- };
-
- /**
- * Where height information comes from
- */
- HeightInfo height_info;
-
- /**
- * If and how cells are shaded
- */
- CellShading cell_shading;
-
- /**
- * Color selection when shading cell
- * boundaries.
- */
- CellBoundaryShading cell_boundary_shading;
-
- /**
- * Number of the vector which is to be
- * used for the height information, within
- * the list of DoF data vectors.
- */
- unsigned int height_vector;
-
- /**
- * Number of the vector which is to be
- * used for the cell shading values, within
- * the list of cell data vectors.
- */
- unsigned int cell_vector;
-
- /**
- * Azimuth of the spectators position.
- * This defines the position on the
- * x-y-plane and is an angle of rotation
- * around the z-axis. We define that
- * if the azimuth angle is zero, this
- * means that the observer is sitting on
- * (or above or below) the positive y-axis
- * and looks back to the origin (direction
- * of view is always towards the origin).
- * Positive angles denote rotation in
- * clockwise sense (as viewed from the top),
- * i.e. 90 degrees would be sitting on the
- * positive x-axis, etc.
- *
- * Please note that the angle has to be
- * given in radians rather than in degrees.
- */
- float azimuth;
-
- /**
- * Elevation of the spectators position.
- * This is the angle that the line between
- * the origin and the spectators position
- * forms with the x-y-plane, measured
- * upwards in direction towards the
- * positive z-axis.
- *
- * Please note that the angle has to be
- * given in radians rather than in degrees.
- */
- float elevation;
-
- /**
- * Direction of the light beams.
- */
- double light[3];
-
- /**
- * Default constructor. Sets height and
- * shading flags to their default values,
- * azimut and elevation to 0.2 each, and
- * puts the light source at @p{(-1,-1,0)}.
- */
- EpsOutputData();
-
- /**
- * Function returning a color value in
- * rgb variables corresponding to the
- * given value @p{x}. @p{x} is considered
- * to be a values within the range
- * @p{xmin...xmax}.
- */
- void color (const float x,
- const float xmax,
- const float xmin,
- float &r,
- float &g,
- float &b) const;
-};
-
-
-#endif
-
+++ /dev/null
-//---------------------------- data_io.cc ---------------------------
-// $Id$
-// Version: $Name$
-//
-// Copyright (C) 1998, 1999, 2000 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.
-//
-//---------------------------- data_io.cc ---------------------------
-
-
-#include <numerics/data_io.h>
-#include <dofs/dof_handler.h>
-#include <dofs/dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria.h>
-#include <grid/geometry_info.h>
-#include <base/quadrature_lib.h>
-#include <fe/fe_values.h>
-#include <fe/fe.h>
-#include <lac/vector.h>
-#include <map>
-#include <iostream>
-#include <iomanip>
-#include <algorithm>
-#include <ctime>
-#include <set>
-#include <cmath>
-
-
-template <int dim>
-DataOut_Old<dim>::DataEntry::DataEntry () :
- data(0), name(""), units("") {};
-
-
-template <int dim>
-DataOut_Old<dim>::DataEntry::DataEntry (const Vector<double> *data,
- const string name,
- const string units) :
- data(data), name(name), units(units) {};
-
-
-template <int dim>
-DataOut_Old<dim>::DataOut_Old () :
- dofs(0) {};
-
-
-template <int dim>
-void DataOut_Old<dim>::attach_dof_handler (const DoFHandler<dim> &d) {
- dofs = &d;
-};
-
-
-template <int dim>
-void DataOut_Old<dim>::add_data_vector (const Vector<double> &vec,
- const string &name,
- const string &units) {
- Assert (dofs != 0, ExcNoDoFHandlerSelected ());
- Assert (name.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
- "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
- "0123456789_<>()") == string::npos,
- ExcInvalidCharacter (name));
- Assert (units.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
- "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
- "0123456789_<>()") == string::npos,
- ExcInvalidCharacter (units));
-
- DataEntry new_entry (&vec, name, units);
- if (vec.size() == dofs->n_dofs())
- dof_data.push_back (new_entry);
- else
- if (vec.size() == dofs->get_tria().n_active_cells())
- cell_data.push_back (new_entry);
- else
- Assert (false,
- ExcInvalidVectorSize (vec.size(),
- dofs->n_dofs(),
- dofs->get_tria().n_active_cells()));
-};
-
-
-template <int dim>
-void DataOut_Old<dim>::clear_data_vectors () {
- dof_data.erase (dof_data.begin(), dof_data.end());
- cell_data.erase (cell_data.begin(), cell_data.end());
-};
-
-
-template <int dim>
-void DataOut_Old<dim>::write_ucd (ostream &out) const {
- Assert (dofs != 0, ExcNoDoFHandlerSelected());
- Assert (dofs->get_fe().dofs_per_vertex==1,
- ExcIncorrectDofsPerVertex());
- Assert ((1<=dim) && (dim<=3), ExcNotImplemented());
-
- DoFHandler<dim>::active_cell_iterator cell;
- const DoFHandler<dim>::active_cell_iterator endc = dofs->end();
- unsigned int n_vertex_dofs = 0;
-
- // first loop over all cells to
- // find out how many degrees of
- // freedom there are located on
- // vertices
- if (true)
- {
- // block this to have local
- // variables destroyed after
- // use
- vector<bool> is_vertex_dof (dofs->n_dofs(), false);
- for (cell=dofs->begin_active(); cell!=endc; ++cell)
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
- ++vertex)
- is_vertex_dof[cell->vertex_dof_index(vertex,0)] = true;
-
- for (unsigned i=0; i!=is_vertex_dof.size(); ++i)
- if (is_vertex_dof[i] == true)
- ++n_vertex_dofs;
- };
-
-
- // 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_vertex_dofs << ' '
- << dofs->get_tria().n_active_cells() + n_boundary_faces() << ' '
- << dof_data.size() << ' '
- << cell_data.size() << ' '
- << 0 // no model data
- << endl;
-
- // write used nodes
- if (true)
- {
- // note if a given vertex was
- // already written
- vector<bool> already_written (dofs->n_dofs(), false);
- // write vertices
- for (cell=dofs->begin_active(); cell!=endc; ++cell)
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
- ++vertex)
- if (already_written[cell->vertex_dof_index(vertex,0)]==false)
- {
- out << cell->vertex_dof_index(vertex,0) // vertex index
- << " ";
- for (unsigned int d=0; d<dim; ++d) // vertex coordinates
- out << cell->vertex(vertex)(d) << ' ';
- for (unsigned int d=dim; d<3; ++d)
- out << 0 << ' ';
- out << endl;
-
- already_written[cell->vertex_dof_index(vertex,0)] = true;
- };
- };
-
- // write cells. Enumerate cells
- // consecutively
- if (true)
- {
- unsigned int index;
- for (cell=dofs->begin_active(), index=0; cell!=endc; ++cell, ++index)
- {
- out << 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
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
- ++vertex)
- out << cell->vertex_dof_index(vertex,0) << ' ';
- out << endl;
- };
-
- write_ucd_faces (out, index);
- };
-
- // if data given: write data, else
- // only write grid
- if (dof_data.size() != 0)
- {
- out << dof_data.size() << " "; // number of vectors
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << 1 << ' '; // number of components;
- // only 1 supported presently
- out << endl;
-
- for (unsigned int i=0; i<dof_data.size(); ++i)
- out << dof_data[i].name << ',' << dof_data[i].units << endl;
-
- // AVS requires that the dof values
- // be given in exactly the same order
- // as in the vertex section and only
- // once.
-
- // note if a given vertex was
- // already written
- vector<bool> already_written (dofs->n_dofs(), false);
- // write vertices
- for (cell=dofs->begin_active(); cell!=endc; ++cell)
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
- ++vertex)
- if (already_written[cell->vertex_dof_index(vertex,0)]==false)
- {
- out << cell->vertex_dof_index(vertex,0) // vertex index
- << " ";
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << (*dof_data[i].data)(cell->vertex_dof_index(vertex,0)) << ' ';
- out << endl;
-
- already_written[cell->vertex_dof_index(vertex,0)] = true;
- };
- };
-
- if (cell_data.size() != 0)
- {
- out << cell_data.size() << " "; // number of vectors
- for (unsigned int i=0; i!=cell_data.size(); ++i)
- out << 1 << ' '; // number of components;
- // only 1 supported presently
- out << endl;
-
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << cell_data[i].name << ',' << cell_data[i].units << endl;
-
- unsigned int index;
- for (cell=dofs->begin_active(), index=0; cell!=endc; ++cell, ++index)
- {
- out << index << " ";
- for (unsigned int i=0; i!=cell_data.size(); ++i)
- out << (*cell_data[i].data)(index) << ' ';
- out << endl;
- };
-
- // strange enough, but true: the ucd
- // format requires that the number of
- // cell data entries be the same as
- // there were cells. cells however
- // include those boundary faces with
- // an indicator other than zero, so
- // we may have printed faces as well
- // in the above list of cells. we
- // have to give respective values
- // here as well. since faces have no
- // natural value when cell data is
- // concerned, we assign a zero.
- for (unsigned int i=0; i<n_boundary_faces(); ++i, ++index)
- {
- out << index << " ";
- for (unsigned int j=0; j!=cell_data.size(); ++j)
- out << "0 ";
- out << endl;
- };
-
- };
-
- // no model data
-
- // assert the stream is still ok
- AssertThrow (out, ExcIO());
-};
-
-
-#if deal_II_dimension == 1
-
-template <>
-unsigned int DataOut_Old<1>::n_boundary_faces () const {
- return 0;
-};
-
-#endif
-
-
-template <int dim>
-unsigned int DataOut_Old<dim>::n_boundary_faces () const {
- typename DoFHandler<dim>::active_face_iterator face, endf;
- unsigned long int n_faces = 0;
-
- for (face=dofs->begin_active_face(), endf=dofs->end_face();
- face != endf; ++face)
- if ((face->boundary_indicator() != 0) &&
- (face->boundary_indicator() != 255))
- n_faces++;
-
- return n_faces;
-};
-
-
-#if deal_II_dimension == 1
-
-template <>
-void DataOut_Old<1>::write_ucd_faces (ostream &, const unsigned int) const {
- return;
-};
-
-#endif
-
-
-template <int dim>
-void DataOut_Old<dim>::write_ucd_faces (ostream &out,
- const unsigned int starting_index) const {
- typename DoFHandler<dim>::active_face_iterator face, endf;
- unsigned int index=starting_index;
-
- for (face=dofs->begin_active_face(), endf=dofs->end_face();
- face != endf; ++face)
- if ((face->boundary_indicator() != 0) &&
- (face->boundary_indicator() != 255))
- {
- 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_dof_index(vertex,0) << ' ';
- out << endl;
-
- ++index;
- };
-
- AssertThrow (out, ExcIO());
-};
-
-
-template <int dim>
-void DataOut_Old<dim>::write_gnuplot (ostream &out, unsigned int accuracy) const
-{
- Assert (dofs != 0, ExcNoDoFHandlerSelected());
- Assert ((1<=dim) && (dim<=3), ExcNotImplemented());
-
- if (dim==1)
- Assert (dofs->get_fe().dofs_per_vertex==1,
- ExcIncorrectDofsPerVertex());
-
- // 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 GNUPLOT format see the GNUPLOT manual."
- << endl
- << "#" << endl
- << "# ";
- switch (dim)
- {
- case 1:
- out << "<x> ";
- break;
- case 2:
- out << "<x> <y> ";
- break;
- case 3:
- out << "<x> <y> <z>";
- break;
-
- default:
- Assert (false, ExcNotImplemented());
- };
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << '<'
- << dof_data[i].name << ':' << dof_data[i].units
- << "> ";
- for (unsigned int i=0; i!=cell_data.size(); ++i)
- out << '<'
- << cell_data[i].name << ':' << cell_data[i].units
- << "> ";
- out << endl;
-
- };
-
-
- DoFHandler<dim>::active_cell_iterator cell;
- DoFHandler<dim>::active_cell_iterator endc = dofs->end();
-
- QTrapez<1> q_trapez;
- QIterated<dim> points (q_trapez, accuracy);
-
- FEValues<dim> fe(dofs->get_fe(), points,
- UpdateFlags(update_values | update_q_points));
- vector< vector <Vector<double> > >
- values (dof_data.size(),
- vector< Vector<double> >(points.n_quadrature_points,
- Vector<double>(dofs->get_fe().n_components()
- )));
-
- unsigned int cell_index=0;
- for (cell=dofs->begin_active(); cell!=endc; ++cell, ++cell_index)
- {
- fe.reinit(cell);
-
- for (unsigned i=0; i<dof_data.size(); ++i)
- fe.get_function_values(*dof_data[i].data, values[i]);
-
- unsigned supp_pt;
-
- switch (dim)
- {
- case 1:
- // one dimension: write
- // each vertex, right vertex
- // and data values
- for (supp_pt = 0; supp_pt<points.n_quadrature_points; ++supp_pt)
- {
- Point<dim> pt = fe.quadrature_point(supp_pt);
- out << pt << " ";
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- for (unsigned int j=0; j < dofs->get_fe().n_components(); ++j)
- out << values[i][supp_pt](j)
- << ' ';
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << (*cell_data[i].data)(cell_index)
- << ' ';
-
- out << endl;
- };
-
- break;
-
- case 2:
- // two dimension: output
- // grid and data as a sequence
- // of points in 3d forming
- // a tensor mesh on each cell
- for (unsigned xpt = 0, supp_pt = 0;
- xpt <= accuracy; ++xpt)
- {
- for(unsigned ypt = 0; ypt <= accuracy; ++ypt, ++supp_pt)
- {
- Point<dim> pt = fe.quadrature_point(supp_pt);
- out << pt << " ";
-
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- for (unsigned int j=0; j < dofs->get_fe().n_components(); ++j)
- out << values[i][supp_pt](j)
- << ' ';
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << (*cell_data[i].data)(cell_index)
- << ' ';
- out << endl;
- }
- out << endl;
- }
-
- break;
-
- case 3:
- // three dimension: output
- // grid and data as a tensor
- // mesh in 4d. we very
- // obviously have a problem
- // here since humans can watch
- // 4d lines properly. but one
- // can at least view cuts
- // through this hypercube,
- // which is possible this
- // way
-
- for (unsigned xpt = 0, supp_pt = 0;
- xpt <= accuracy; ++xpt)
- {
- for(unsigned ypt = 0; ypt <= accuracy; ++ypt)
- {
- for(unsigned zpt = 0; zpt <= accuracy; ++zpt, ++supp_pt)
- {
- Point<dim> pt = fe.quadrature_point(supp_pt);
- out << pt << " ";
-
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- for (unsigned int j=0; j < dofs->get_fe().n_components(); ++j)
- out << values[i][supp_pt](j)
- << ' ';
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << (*cell_data[i].data)(cell_index)
- << ' ';
- out << endl;
- }
- out << endl;
- }
- out << endl;
- };
-
-
- break;
-
- default:
- Assert (false, ExcNotImplemented());
- }
- out << endl;
- }
-
- AssertThrow (out, ExcIO());
-}
-
-
-template <int dim>
-void DataOut_Old<dim>::write_gnuplot_draft (ostream &out) const
-{
- Assert (dofs != 0, ExcNoDoFHandlerSelected());
- Assert ((1<=dim) && (dim<=3), ExcNotImplemented());
-
- if (dim==1)
- Assert (dofs->get_fe().dofs_per_vertex==1,
- ExcIncorrectDofsPerVertex());
-
- // 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 GNUPLOT format see the GNUPLOT manual."
- << endl
- << "#" << endl
- << "# ";
- switch (dim)
- {
- case 1:
- out << "<x> ";
- break;
- case 2:
- out << "<x> <y> ";
- break;
- case 3:
- out << "<x> <y> <z>";
- break;
- default:
- Assert (false, ExcNotImplemented());
- };
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << '<'
- << dof_data[i].name << ':' << dof_data[i].units
- << "> ";
- out << endl;
-
- };
-
-
- DoFHandler<dim>::active_cell_iterator cell;
- DoFHandler<dim>::active_cell_iterator endc = dofs->end();
-
- unsigned int cell_index=0;
- for (cell=dofs->begin_active(); cell!=endc; ++cell, ++cell_index)
- {
- switch (dim)
- {
- case 1:
- // one dimension: write
- // left vertex, right vertex
- // and data values
- for (unsigned int vertex=0; vertex<2; ++vertex)
- {
- out << cell->vertex(vertex)
- << " ";
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << (*dof_data[i].data)(cell->vertex_dof_index(vertex,0))
- << ' ';
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << (*cell_data[i].data)(cell_index)
- << ' ';
- out << endl;
- };
-
- break;
-
- case 2:
- // two dimension: output
- // grid and data as a sequence
- // of lines in 3d
- for (unsigned int vertex=0; vertex<4; ++vertex)
- {
- out << cell->vertex(vertex) << " ";
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << (*dof_data[i].data)(cell->vertex_dof_index(vertex,0))
- << ' ';
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << (*cell_data[i].data)(cell_index)
- << ' ';
- out << endl;
- };
- // first vertex again
- out << cell->vertex(0) << " ";
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << (*dof_data[i].data)(cell->vertex_dof_index(0,0))
- << ' ';
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << (*cell_data[i].data)(cell_index)
- << ' ';
- out << endl
- << endl
- << endl; // end of cell; two newlines, since this
- // stops continuous drawing
- // of lines
-
- break;
-
- case 3:
- // three dimension: output
- // grid and data as a sequence
- // of lines in 4d
-
- // first we plot the front face
- for (unsigned int vertex=0; vertex<4; ++vertex)
- {
- out << cell->vertex(vertex) << " ";
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << (*dof_data[i].data)(cell->vertex_dof_index(vertex,0))
- << ' ';
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << (*cell_data[i].data)(cell_index)
- << ' ';
- out << endl;
- };
- // first vertex again
- out << cell->vertex(0) << " ";
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << (*dof_data[i].data)(cell->vertex_dof_index(0,0))
- << ' ';
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << (*cell_data[i].data)(cell_index)
- << ' ';
- out << endl
- << endl
- << endl; // end of front face; two newlines, since this
- // stops continuous drawing
- // of lines
-
- // first we back the front face
- for (unsigned int vertex=4; vertex<8; ++vertex)
- {
- out << cell->vertex(vertex) << " ";
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << (*dof_data[i].data)(cell->vertex_dof_index(vertex,0))
- << ' ';
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << (*cell_data[i].data)(cell_index)
- << ' ';
- out << endl;
- };
- // first vertex again
- out << cell->vertex(4) << " ";
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << (*dof_data[i].data)(cell->vertex_dof_index(4,0))
- << ' ';
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << (*cell_data[i].data)(cell_index)
- << ' ';
- out << endl
- << endl
- << endl; // end of back face; two newlines, since this
- // stops continuous drawing
- // of lines
-
- // finally, we plot a
- // continuous line along
- // the vertices 0, 4, 7, 3, 2,
- // 6, 5, 1 to show the other
- // four lines
- for (unsigned int vertex=0; vertex<8; ++vertex)
- {
- static const unsigned int vertex_list[8]
- = { 0, 4, 7, 3, 2, 6, 5, 1 };
-
- out << cell->vertex(vertex_list[vertex]) << " ";
- for (unsigned int i=0; i!=dof_data.size(); ++i)
- out << (*dof_data[i].data)(cell->vertex_dof_index(vertex_list[vertex],0))
- << ' ';
- for (unsigned int i=0; i<cell_data.size(); ++i)
- out << (*cell_data[i].data)(cell_index)
- << ' ';
- out << endl;
- };
- // again: stop drawing
- out << endl << endl;
-
-
- break;
-
- default:
- Assert (false, ExcNotImplemented());
- };
- };
-
- AssertThrow (out, ExcIO());
-};
-
-
-#if deal_II_dimension == 2
-
-template <>
-void DataOut_Old<2>::write_povray_mesh (ostream &out) const {
- Assert (dofs != 0, ExcNoDoFHandlerSelected());
-
- // write preamble
- if (true)
- {
- // block this to have local
- // variables destroyed after
- // use
- time_t time1= time (0);
- tm *time = localtime(&time1);
- out << "/*" << endl
- << "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 POVRAY format see the POVRAY manual."
- << endl
- << endl
- << "*/"
- << endl << endl;
- };
-
- // write list of include files
- out << "#include \"colors.inc\"" << endl
- << endl;
-
- // declare standard texture
- out << "#declare default_texture = texture {" << endl
- << " pigment { color White }" << endl
- << " finish { ambient 0.2 diffuse 0.6 specular 0.5 }" << endl
- << "}" << endl << endl;
-
- // write camera and light sources
- out << "camera {" << endl
- << " location <8, 10, -20>" << endl
- << " angle 7" << endl
- << " look_at <0., 0., 0.>" << endl
- << " sky <0., 0., 1.>" << endl
- << "}" << endl
- << endl
- << "light_source {" << endl
- << " <20, 20, -20> White" << endl
- << "}" << endl
- << endl;
-
-
- // write frame object
- out << "mesh {" << endl;
-
- DoFHandler<2>::active_cell_iterator cell;
- DoFHandler<2>::active_cell_iterator endc = dofs->end();
-
- for (cell=dofs->begin_active(); cell!=endc; ++cell)
- {
- // write cell as two triangles
- // y and z coordinates are switched
- out << " triangle { ";
- out << '<' << cell->vertex(0)(0) << ','
- << (*dof_data[0].data)(cell->vertex_dof_index(0,0)) << ','
- << cell->vertex(0)(1) << '>'
- << ", "
- << '<' << cell->vertex(1)(0) << ','
- << (*dof_data[0].data)(cell->vertex_dof_index(1,0)) << ','
- << cell->vertex(1)(1) << '>'
- << ", "
- << '<' << cell->vertex(2)(0) << ','
- << (*dof_data[0].data)(cell->vertex_dof_index(2,0)) << ','
- << cell->vertex(2)(1) << '>'
- << endl
- << " texture { default_texture }" << endl
- << " }"
- << endl;
- out << " triangle { ";
- out << '<' << cell->vertex(0)(0) << ','
- << (*dof_data[0].data)(cell->vertex_dof_index(0,0)) << ','
- << cell->vertex(0)(1) << '>'
- << ", "
- << '<' << cell->vertex(2)(0) << ','
- << (*dof_data[0].data)(cell->vertex_dof_index(2,0)) << ','
- << cell->vertex(2)(1) << '>'
- << ", "
- << '<' << cell->vertex(3)(0) << ','
- << (*dof_data[0].data)(cell->vertex_dof_index(3,0)) << ','
- << cell->vertex(3)(1) << '>'
- << endl
- << " texture { default_texture }" << endl
- << " }"
- << endl;
- };
- out << "}";
-
- AssertThrow (out, ExcIO());
-};
-
-
-#endif
-
-
-template <int dim>
-void DataOut_Old<dim>::write_povray_mesh (ostream &) const {
- // this is for all other dimensions that
- // are not explicitely specialized
- Assert (false, ExcNotImplemented());
-};
-
-
-#if deal_II_dimension == 2
-
-template <>
-void DataOut_Old<2>::write_eps (ostream &out, const EpsOutputData &eod) const {
- Assert (dofs != 0, ExcNoDoFHandlerSelected());
-
- {
- // write preamble
- // block this to have local
- // variables destroyed after
- // use
- time_t time1= time (0);
- tm *time = localtime(&time1);
- out << "%!PS-Adobe-2.0 EPSF-1.2" << endl
- << "%%Title: deal.II Output" << endl
- << "%%Creator: the deal.II library" << endl
- << "%%Creation Date: "
- << time->tm_year+1900 << "/"
- << time->tm_mon+1 << "/"
- << time->tm_mday << " - "
- << time->tm_hour << ":"
- << setw(2) << time->tm_min << ":"
- << setw(2) << time->tm_sec << endl
- << "%%BoundingBox: -5 -5 305 305" << endl;
- };
-
-
- // Make output values local by
- // copying them to a multiset.
- // Perform the necessary turn.
- const DoFHandler<2>::active_cell_iterator endc = dofs->end();
- multiset<DataOut_Old<2>::EpsCellData> cells;
- multiset<DataOut_Old<2>::EpsCellData> cells2;
-
- bool height_data_p = (
- ((dof_data.size())>0)
- &&
- (
- ( eod.height_info == EpsOutputData::DefaultHeight)
- ||
- ( eod.height_info == EpsOutputData::HeightVector)
- )
- );
-
- // Cells are colored, if there is
- // cell data and the mode is
- // ShadingVector or DefaultShading
- bool cell_data_p = (
- ((cell_data.size())>0)
- &&
- (
- (eod.cell_shading == EpsOutputData::ShadingVector)
- ||
- (eod.cell_shading == EpsOutputData::DefaultShading)
- )
- );
-
- // Cells are shaded, i.e. light
- // shading, if they are not
- // colored, there is height
- // information and the mode is
- // Default or Light
- bool cell_shade_p = (
- (!cell_data_p)
- &&
- (height_data_p)
- &&
- (
- (eod.cell_shading == EpsOutputData::DefaultShading)
- ||
- (eod.cell_shading == EpsOutputData::LightShaded)
- )
- );
-
- unsigned cell_index;
- DoFHandler<2>::active_cell_iterator cell;
- for(cell_index=0, cell=dofs->begin_active();
- cell!=endc;
- ++cell, ++cell_index)
- {
- EpsCellData cd;
- for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
- {
- (cd.vertices[i]).x=cell->vertex(i)(0);
- (cd.vertices[i]).y=cell->vertex(i)(1);
- (cd.vertices[i]).z=(height_data_p ?
- (*dof_data[eod.height_vector].data)(cell->vertex_dof_index(i,0))
- : 0);
- };
- if (height_data_p)
- cd.turn(eod.azimuth,eod.elevation);
-
- if (cell_data_p)
- cd.red=(*cell_data[eod.cell_vector].data)(cell_index);
- cells.insert(cd);
- };
-
- // Now we proceed with the
- // multiset cells. First we look
- // for extrema.
-
- double xmin=cells.begin()->vertices[0].x;
- double xmax=xmin;
- double ymin=cells.begin()->vertices[0].y;
- double ymax=ymin;
- float cell_vector_min=cells.begin()->red;
- float cell_vector_max=cell_vector_min;
-
- for(multiset<DataOut_Old<2>::EpsCellData>::iterator c=cells.begin();
- c!=cells.end(); ++c, ++cell_index)
- {
- for (unsigned int i=0; i<4; ++i)
- {
- double xvv,yvv;
- xvv=c->vertices[i].x;
- xmin=(xmin < xvv ? xmin : xvv);
- xmax=(xmax > xvv ? xmax : xvv);
-
- yvv=c->vertices[i].y;
- ymin=(ymin < yvv ? ymin : yvv);
- ymax=(ymax > yvv ? ymax : yvv);
- }
- if (cell_data_p)
- {
- double cvv;
- cvv = c->red;
- cell_vector_max = (cell_vector_max > cvv ? cell_vector_max : cvv);
- cell_vector_min = (cell_vector_min < cvv ? cell_vector_min : cvv);
- }
- };
- cells2.clear();
-
-
- // If we want shaded output we can
- // do the shading now.
- if (cell_shade_p)
- {
- double spann1[3], spann2[3], normal[3];
- double light_norm, normal_norm;
- float color;
-
- for (multiset<DataOut_Old<2>::EpsCellData>::iterator c=cells.begin();
- c!=cells.end(); ++c)
- {
- EpsCellData cd(*c);
-
- spann1[0]=spann2[0]=cd.vertices[0].x;
- spann1[1]=spann2[1]=cd.vertices[0].y;
- spann1[2]=spann2[2]=cd.vertices[0].z;
-
- spann1[0]-=cd.vertices[1].x;
- spann1[1]-=cd.vertices[1].y;
- spann1[2]-=cd.vertices[1].z;
-
- spann2[0]-=cd.vertices[2].x;
- spann2[1]-=cd.vertices[2].y;
- spann2[2]-=cd.vertices[2].z;
-
- normal[0] = spann1[1]*spann2[2]-spann1[2]*spann2[1];
- normal[1] = spann1[2]*spann2[0]-spann1[0]*spann2[2];
- normal[2] = spann1[0]*spann2[1]-spann1[1]*spann2[0];
-
- normal_norm = sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
- light_norm = sqrt(eod.light[0]*eod.light[0]+eod.light[1]*eod.light[1]+eod.light[2]*eod.light[2]);
-
- color = eod.light[0]*normal[0]+eod.light[1]*normal[1]+eod.light[2]*normal[2];
- color /= light_norm * normal_norm;
-
- cd.red=color;
- cd.green=color;
- cd.blue=color;
-
- cells2.insert(cd);
- };
-
- // since we don't need @p{cells} any
- // more, delete it now
- cells.clear ();
- }
- else
- // copy @p{cells} to @p{cells2}. since
- // we don't need @p{cells} any
- // more, we use a trick for copying
- // that is significantly faster
- cells2.swap (cells);
-
-
- // Next we have to shift and scale
- // a bit so that everything still
- // arrives in our bounding box of
- // 310x310.
- // If cell_data_p we also scale
- // this so that it is in the range
- // between 0 and 1.
-
- const double scale = 300 / (xmax-xmin > ymax-ymin ? xmax-xmin : ymax-ymin);
-
- for (multiset<DataOut_Old<2>::EpsCellData>::iterator c=cells2.begin();
- c!=cells2.end(); ++c)
- {
- EpsCellData cd (*c);
- for (unsigned int i=0; i<4; ++i)
- {
- cd.vertices[i].x=(cd.vertices[i].x-xmin)*scale;
- cd.vertices[i].y=(cd.vertices[i].y-ymin)*scale;
- };
-
- if (cell_data_p)
- eod.color(cd.red,cell_vector_max,cell_vector_min,cd.red,cd.green,cd.blue);
-
- cells.insert(cd);
- };
-
-
- // Now we are ready to output...
- for (multiset<DataOut_Old<2>::EpsCellData>::iterator c=cells.begin();
- c!=cells.end(); ++c)
- {
- if (cell_data_p || cell_shade_p)
- {
- out << c->red << " " << c->green << " " << c->blue << " setrgbcolor "
- << c->vertices[0].x << " " << c->vertices[0].y << " moveto "
- << c->vertices[1].x << " " << c->vertices[1].y << " lineto "
- << c->vertices[2].x << " " << c->vertices[2].y << " lineto "
- << c->vertices[3].x << " " << c->vertices[3].y << " lineto "
- << " closepath fill" << endl;
- };
-
- if (eod.cell_boundary_shading != EpsOutputData::NoBoundary)
- {
- switch (eod.cell_boundary_shading)
- {
- case EpsOutputData::BlackBoundary:
- case EpsOutputData::DefaultBoundary:
- out << "0";
- break;
- case EpsOutputData::WhiteBoundary:
- out << "1";
- break;
- case EpsOutputData::NoBoundary:
- break;
- };
- out << " setgray "
- << c->vertices[0].x << " " << c->vertices[0].y << " moveto "
- << c->vertices[1].x << " " << c->vertices[1].y << " lineto "
- << c->vertices[2].x << " " << c->vertices[2].y << " lineto "
- << c->vertices[3].x << " " << c->vertices[3].y << " lineto closepath stroke" << endl;
- };
- };
- out << "showpage" << endl;
-};
-
-#endif
-
-
-template <int dim>
-void DataOut_Old<dim>::write_eps (ostream &,
- const EpsOutputData &) const{
- // this is for all other dimensions that
- // are not explicitely specialized
- Assert (false, ExcNotImplemented());
-};
-
-
-template <int dim>
-void DataOut_Old<dim>::write_gmv (ostream &out) const
-{
- // this function is mostly copied from
- // the ucd format function
-
- Assert (dofs != 0, ExcNoDoFHandlerSelected());
- Assert (dofs->get_fe().dofs_per_vertex==1,
- ExcIncorrectDofsPerVertex());
- Assert ((1<=dim) && (dim<=3), ExcNotImplemented());
-
- const unsigned int dofs_per_vertex = dofs->get_fe().dofs_per_vertex;
-
- ///////////////////////
- // preamble
- out << "gmvinput ascii"
- << endl
- << endl;
-
-
- ///////////////////////////////
- // first make up a list of used
- // vertices along with their
- // coordinates. since dofs have
- // no separate numbering in
- // gmv format, we need to generate
- // a separate numbering as well
- // for those dofs on vertices.
- // note that we can't simply use
- // the numbers of the vertices,
- // since this numbering may not
- // always be continuous if the
- // triangulation has undergone
- // coarsening somewhen (well, we
- // could place all unused vertices
- // to the origin, but that's not
- // really elegant...)
- //
- // therefore: have a list with
- // one slot for each dof, where
- // the value indicates the
- // respective number of the vertex
- // if the dof is on a vertex,
- // and -1 otherwise.
- //
- // also create a list of vertices
- // storing their coordinates in
- // the order defined by the
- // above map, which we write to
- // the file immediately after
- // creation.
- //
- // note: since GMV seems to be
- // made by Fortranists, the count
- // their indices from 1 onwards,
- // so whenever we actually output
- // a vertex number, we should add
- // a one.
- DoFHandler<dim>::active_cell_iterator cell;
- const DoFHandler<dim>::active_cell_iterator endc = dofs->end();
-
- vector<unsigned int> dof_to_vertex_map (dofs->n_dofs(), DoFHandler<dim>::invalid_dof_index);
- unsigned int used_vertices = 0;
- if (true)
- {
- vector<Point<dim> > vertices (dofs->get_tria().n_used_vertices());
- unsigned int next_free_vertex = 0;
-
- for (cell=dofs->begin_active(); cell!=endc; ++cell)
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
- // check whether we
- // already treated this vertex
- if (dof_to_vertex_map[cell->vertex_dof_index(vertex,0)] ==
- DoFHandler<dim>::invalid_dof_index)
- {
- vertices[next_free_vertex] = cell->vertex(vertex);
-
- // associate all dofs on
- // this vertex with the
- // respective vertex
- for (unsigned int d=0; d<dofs_per_vertex; ++d)
- dof_to_vertex_map[cell->vertex_dof_index(vertex,d)] = next_free_vertex;
-
- ++next_free_vertex;
- };
- Assert (next_free_vertex == vertices.size(),
- ExcInternalError());
- used_vertices = vertices.size();
-
- // now write out the vertices in
- // this order
- out << "nodes " << vertices.size() << endl;
- for (unsigned int component=0; component<dim; ++component)
- {
- for (unsigned int v=0; v<vertices.size(); ++v)
- out << vertices[v](component) << ' ';
- out << endl;
- };
-
- // and write the missing components
- // y (for 1d) and z (for 1d and 2d)
- // by simply setting them to zero
- for (unsigned int d=dim+1; d<=3; ++d)
- {
- fill_n (ostream_iterator<double>(out, " "), vertices.size(), 0.0);
- out << endl;
- }
-
- out << endl;
- };
-
-
- /////////////////////////////////////
- // now for the cells. this is simpler
- // than the above task
- if (true)
- {
- out << "cells " << dofs->get_tria().n_active_cells() << endl;
-
- const char *cell_description[3] = { "line 2\n ",
- "quad 4\n ",
- "hex 8\n "};
-
- for (cell=dofs->begin_active(); cell!=endc; ++cell)
- {
- out << cell_description[dim-1];
- // output vertex indices,
- // counted from 1 onwards
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- out << dof_to_vertex_map[cell->vertex_dof_index(v,0)]+1 << ' ';
- out << endl << endl;
- };
-
- out << endl;
- };
-
-
- ///////////////////////////////////////
- // data output.
- out << "variable" << endl;
-
- // first go for node data
- //
- // since here as with the vertex
- // coordinates the order is a bit
- // unpleasant (first all data of
- // variable 1, then variable 2, etc)
- // we have to copy them a bit around
- //
- // note that we copy vectors when
- // looping over the cells since we
- // have to write them one variable
- // at a time and don't want to use
- // more than one loop
- if (true)
- {
- vector<vector<double> > data_vectors (dof_data.size(),
- vector<double> (used_vertices));
-
- // loop over all cells and copy
- // the data into the other vector
- // note if a vertex has already
- // been visited
- vector<bool> vertex_copied (dofs->n_dofs(), false);
- for (cell=dofs->begin_active(); cell!=endc; ++cell)
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- if (vertex_copied[cell->vertex_dof_index(v,0)] == false)
- {
- // global dof number
- const unsigned int dof_index = cell->vertex_dof_index(v,0);
-
- // this is this vertex's
- // number to GMV
- const int vertex_index = dof_to_vertex_map[dof_index];
- Assert (vertex_index >= 0, ExcInternalError());
-
- for (unsigned vec=0; vec<dof_data.size(); ++vec)
- data_vectors[vec][vertex_index] = (*dof_data[vec].data)(dof_index);
-
- vertex_copied[dof_index] = true;
- };
-
- // indicator 0 = node data
- for (unsigned int vec=0; vec<dof_data.size(); ++vec)
- {
- out << dof_data[vec].name << " 1" << endl;
- copy(data_vectors[vec].begin(),
- data_vectors[vec].end(),
- ostream_iterator<double>(out, " "));
- out << endl
- << endl;
- };
- };
-
-
- // Cell data is the simplest since the order
- // already is correct
- if (true)
- {
- for (unsigned int vec=0; vec<cell_data.size(); ++vec)
- {
- out << cell_data[vec].name << " 0" << endl;
- copy((*cell_data[vec].data).begin(),
- (*cell_data[vec].data).end(),
- ostream_iterator<double>(out, " "));
- out << endl
- << endl;
- };
- };
-
- // end of variable section
- out << "endvars" << endl;
-
- // end of output
- out << "endgmv"
- << endl;
-
- // assert the stream is still ok
- AssertThrow (out, ExcIO());
-};
-
-
-template <int dim>
-void DataOut_Old<dim>::write (ostream &out,
- const OutputFormat output_format) const {
- switch (output_format)
- {
- case ucd:
- write_ucd (out);
- break;
-
- case gnuplot:
- write_gnuplot (out);
- break;
-
- case gnuplot_draft:
- write_gnuplot_draft (out);
- break;
-
- case povray_mesh:
- write_povray_mesh (out);
- break;
-
- case eps:
- write_eps(out);
- break;
-
- case gmv:
- write_gmv (out);
- break;
-
- default:
- Assert (false, ExcNotImplemented());
- };
-};
-
-
-template <int dim>
-string DataOut_Old<dim>::default_suffix (const OutputFormat output_format)
-{
- switch (output_format)
- {
- case ucd:
- return ".inp";
-
- case gnuplot:
- case gnuplot_draft:
- return ".gnuplot";
-
- case povray_mesh:
- return ".pov";
-
- case eps:
- return ".eps";
-
- case gmv:
- return ".gmv";
-
- default:
- Assert (false, ExcNotImplemented());
- return "";
- };
-};
-
-
-template <int dim>
-DataOut_Old<dim>::OutputFormat
-DataOut_Old<dim>::parse_output_format (const string &format_name) {
- if (format_name == "ucd")
- return ucd;
-
- if (format_name == "gnuplot")
- return gnuplot;
-
- if (format_name == "gnuplot draft")
- return gnuplot_draft;
-
- if (format_name == "povray mesh")
- return povray_mesh;
-
- if (format_name == "eps")
- return eps;
-
- if (format_name == "gmv")
- return gmv;
-
- AssertThrow (false, ExcInvalidState ());
-
- // return something invalid
- return OutputFormat(-1);
-};
-
-
-template <int dim>
-string DataOut_Old<dim>::get_output_format_names () {
- return "ucd|gnuplot|gnuplot draft|povray mesh|eps|gmv";
-};
-
-
-template<int dim>
-bool DataOut_Old<dim>::EpsCellData::operator < (const EpsCellData &other) const
-{
- double maxz = vertices[0].z,
- othermaxz = other.vertices[0].z;
- unsigned i;
-
- for (i=1; i<4; ++i)
- {
- maxz = (maxz > vertices[i].z ? maxz : vertices[i].z);
- othermaxz = (othermaxz > other.vertices[i].z ? othermaxz : other.vertices[i].z);
- };
-
- return maxz > othermaxz;
-};
-
-
-template <int dim>
-void DataOut_Old<dim>::EpsVertexData::turn(double azi, double ele)
-{
- double nx,ny,nz;
-
- double cx=cos(ele), cz=cos(azi), sx=sin(ele), sz=sin(azi);
-
- nx = - cz*x+ sz*y;
- ny = -cx*sz*x-cx*cz*y-sx*z;
- nz = -sx*sz*x-sx*cz*y+cx*z;
-
- x=nx; z=ny; y=nz;
-};
-
-// ( 1 0 0 )
-// Dx = ( 0 cx -sx )
-// ( 0 sx cx )
-
-// ( cy 0 sy )
-// Dy = ( 0 1 0 )
-// (-sy 0 cy )
-
-// ( cz -sz 0 )
-// Dz = ( sz cz 0 )
-// ( 0 0 1 )
-
-// ( cz -sz 0 )( 1 0 0 )(x) ( cz*x-sz*(cx*y-sx*z)+0*(sx*y+cx*z) )
-// Dxz = ( sz cz 0 )( 0 cx -sx )(y) = ( sz*x+cz*(cx*y-sx*z)+0*(sx*y+cx*z) )
-// ( 0 0 1 )( 0 sx cx )(z) ( 0*x+ *(cx*y-sx*z)+1*(sx*y+cx*z) )
-
-
-template <int dim>
-void DataOut_Old<dim>::EpsCellData::turn(double azi, double ele)
-{
- for (unsigned i=0; i<4; ++i)
- vertices[i].turn(azi,ele);
-};
-
-
-EpsOutputData::EpsOutputData()
- : height_info(DefaultHeight),
- cell_shading(DefaultShading),
- cell_boundary_shading (DefaultBoundary),
- height_vector(0),
- cell_vector(0),
- azimuth(2*3.1415926* (180 - 30)/360),
- elevation(2*3.1415926* (90-60)/360)
-{
- light[0]=-1;
- light[1]=-1;
- light[2]=1;
-};
-
-
-void EpsOutputData::color(const float x,
- const float xmax,
- const float xmin,
- float &r,
- float &g,
- float &b) const
-{
-// A difficult color scale:
-// xmin = black (1)
-// 3/4*xmin+1/4*xmax = blue (2)
-// 1/2*xmin+1/2*xmax = green (3)
-// 1/4*xmin+3/4*xmax = red (4)
-// xmax = white (5)
-// Makes the following color functions:
-//
-// red green blue
-// __
-// / /\ / /\ /
-// ____/ __/ \/ / \__/
-
-// { 0 (1) - (3)
-// r = { ( 4*x-2*xmin+2*xmax)/(xmax-xmin) (3) - (4)
-// { 1 (4) - (5)
-//
-// { 0 (1) - (2)
-// g = { ( 4*x-3*xmin- xmax)/(xmax-xmin) (2) - (3)
-// { (-4*x+ xmin+3*xmax)/(xmax-xmin) (3) - (4)
-// { ( 4*x- xmin-3*xmax)/(xmax-xmin) (4) - (5)
-//
-// { ( 4*x-4*xmin )/(xmax-xmin) (1) - (2)
-// b = { (-4*x+2*xmin+2*xmax)/(xmax-xmin) (2) - (3)
-// { 0 (3) - (4)
-// { ( 4*x- xmin-3*xmax)/(xmax-xmin) (4) - (5)
-
- float sum = xmax+ xmin;
- float sum13 = xmin+3*xmax;
- float sum22 = 2*xmin+2*xmax;
- float sum31 = 3*xmin+ xmax;
- float dif = xmax-xmin;
- float rezdif = 1.0/dif;
-
- int where;
-
- if (x<(sum31)/4)
- where = 0;
- else if (x<(sum22)/4)
- where = 1;
- else if (x<(sum13)/4)
- where = 2;
- else
- where = 3;
-
- if (dif!=0)
- {
- switch (where)
- {
- case 0:
- r=0; g=0; b=(x-xmin)*4.*rezdif;
- break;
- case 1:
- r=0; g=(4*x-3*xmin-xmax)*rezdif; b=(sum22-4.*x)*rezdif;
- break;
- case 2:
- r=(4*x-2*sum)*rezdif; g=(xmin+3*xmax-4*x)*rezdif; b=0;
- break;
- case 3:
- r=1; g=(4*x-xmin-3*xmax)*rezdif; b=(4.*x-sum13)*rezdif;
- default:
- break;
- };
- }
- else // White
- {
- r=1;
- g=1;
- b=1;
- };
-};
-
-
-//explicit instantiations
-template class DataOut_Old<deal_II_dimension>;
projects / dealii-svn.git / commitdiff