//---------------------------- data_out_base.cc ---------------------------
-#include \"data.inc\" " << endl;
+#include <base/data_out_base.h>
+#include <base/parameter_handler.h>
+#include <iomanip>
+#include <ctime>
+#include <cmath>
+#include <set>
+
+
+// egcs does not understand this at present.
+//
+// template <int dim>
+// DataOut::Patch<dim>::Patch () :
+// n_subdivisions (1)
+// // all the rest has a constructor of its own
+// {};
+
+
+template <int dim>
+void DataOutBase::write_ucd (const vector<Patch<dim> > &patches,
+ const vector<string> &data_names,
+ const UcdFlags &flags,
+ ostream &out)
+{
+ AssertThrow (out, ExcIO());
+
+ Assert (patches.size() > 0, ExcNoPatches());
+
+ const unsigned int n_data_sets = data_names.size();
+
+ // first count the number of cells
+ // and cells for later use
+ unsigned int n_cells = 0,
+ n_nodes = 0;
+ for (vector<Patch<dim> >::const_iterator patch=patches.begin();
+ patch!=patches.end(); ++patch)
+ switch (dim)
+ {
+ case 1:
+ n_cells += patch->n_subdivisions;
+ n_nodes += patch->n_subdivisions+1;
+ break;
+ case 2:
+ n_cells += patch->n_subdivisions *
+ patch->n_subdivisions;
+ n_nodes += (patch->n_subdivisions+1) *
+ (patch->n_subdivisions+1);
+ break;
+ case 3:
+ n_cells += patch->n_subdivisions *
+ patch->n_subdivisions *
+ patch->n_subdivisions;
+ n_nodes += (patch->n_subdivisions+1) *
+ (patch->n_subdivisions+1) *
+ (patch->n_subdivisions+1);
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+ ///////////////////////
+ // preamble
+ if (flags.write_preamble)
+ {
+ 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_nodes << ' '
+ << n_cells << ' '
+ << n_data_sets << ' '
+ << 0 << ' ' // no cell data at present
+ << 0 // no model data
+ << endl;
+
+ ///////////////////////////////
+ // first make up the list of used
+ // nodes along with their
+ // coordinates. number them
+ // consecutively starting with 1
+ //
+ // note that we have to print
+ // d=1..3 dimensions
+ if (true)
+ {
+ unsigned int present_node = 1;
+
+ for (vector<Patch<dim> >::const_iterator patch=patches.begin();
+ patch!=patches.end(); ++patch)
+ {
+ const unsigned int n_subdivisions = patch->n_subdivisions;
+
+ // if we have nonzero values for
+ // this coordinate
+ switch (dim)
+ {
+ case 1:
+ {
+ for (unsigned int i=0; i<n_subdivisions+1; ++i, ++present_node)
+ out << present_node
+ << " "
+ << ((patch->vertices[1](0) * i / n_subdivisions) +
+ (patch->vertices[0](0) * (n_subdivisions-i) / n_subdivisions))
+ << " 0 0\n"; // fill with zeroes
+ break;
+ };
+
+ case 2:
+ {
+ for (unsigned int i=0; i<n_subdivisions+1; ++i)
+ for (unsigned int j=0; j<n_subdivisions+1; ++j)
+ {
+ const double x_frac = i * 1./n_subdivisions,
+ y_frac = j * 1./n_subdivisions;
+
+ // compute coordinates for
+ // this patch point
+ out << present_node
+ << " "
+ << (((patch->vertices[1] * x_frac) +
+ (patch->vertices[0] * (1-x_frac))) * (1-y_frac) +
+ ((patch->vertices[2] * x_frac) +
+ (patch->vertices[3] * (1-x_frac))) * y_frac)
+ << " 0\n"; // fill with zeroes
+
+ ++present_node;
+ };
+
+ break;
+ };
+
+ case 3:
+ {
+ for (unsigned int i=0; i<n_subdivisions+1; ++i)
+ for (unsigned int j=0; j<n_subdivisions+1; ++j)
+ for (unsigned int k=0; k<n_subdivisions+1; ++k)
+ {
+ // note the broken
+ // design of hexahedra
+ // in deal.II, where
+ // first the z-component
+ // is counted up, before
+ // increasing the y-
+ // coordinate
+ const double x_frac = i * 1./n_subdivisions,
+ y_frac = k * 1./n_subdivisions,
+ z_frac = j * 1./n_subdivisions;
+
+ // compute coordinates for
+ // this patch point
+ out << present_node
+ << " "
+ << ((((patch->vertices[1] * x_frac) +
+ (patch->vertices[0] * (1-x_frac))) * (1-y_frac) +
+ ((patch->vertices[2] * x_frac) +
+ (patch->vertices[3] * (1-x_frac))) * y_frac) * (1-z_frac) +
+ (((patch->vertices[5] * x_frac) +
+ (patch->vertices[4] * (1-x_frac))) * (1-y_frac) +
+ ((patch->vertices[6] * x_frac) +
+ (patch->vertices[7] * (1-x_frac))) * y_frac) * z_frac)
+ << endl;
+
+ ++present_node;
+ };
+
+ break;
+ };
+
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
+
+ // note that we number starting with 1!
+ Assert (present_node == n_nodes+1,
+ ExcInternalError());
+ };
+
+ /////////////////////////////////////////
+ // write cell. number them consecutively,
+ // starting with 1
+ if (true)
+ {
+ unsigned int present_cell = 1;
+ unsigned int first_vertex_of_patch = 0;
+
+ for (vector<Patch<dim> >::const_iterator patch=patches.begin();
+ patch!=patches.end(); ++patch)
+ {
+ const unsigned int n_subdivisions = patch->n_subdivisions;
+
+ // write out the cells making
+ // up this patch
+ switch (dim)
+ {
+ case 1:
+ {
+ for (unsigned int i=0; i<n_subdivisions; ++i, ++present_cell)
+ out << present_cell
+ << " 0 line " // set material id to 0
+ << first_vertex_of_patch+i+1 << ' '
+ << first_vertex_of_patch+i+1+1 << endl;
+ break;
+ };
+
+ case 2:
+ {
+ for (unsigned int i=0; i<n_subdivisions; ++i)
+ for (unsigned int j=0; j<n_subdivisions; ++j)
+ {
+ out << present_cell
+ << " 0 quad " // set material id to 0
+ << first_vertex_of_patch+i*(n_subdivisions+1)+j+1 << ' '
+ << first_vertex_of_patch+(i+1)*(n_subdivisions+1)+j+1 << ' '
+ << first_vertex_of_patch+(i+1)*(n_subdivisions+1)+j+1+1 << ' '
+ << first_vertex_of_patch+i*(n_subdivisions+1)+j+1+1
+ << endl;
+ ++present_cell;
+ };
+ break;
+ };
+
+ case 3:
+ {
+ for (unsigned int i=0; i<n_subdivisions; ++i)
+ for (unsigned int j=0; j<n_subdivisions; ++j)
+ for (unsigned int k=0; k<n_subdivisions; ++k)
+ {
+ out << present_cell
+ << " 0 hex " // set material id to 0
+ // note: vertex indices start with 1!
+ << first_vertex_of_patch+(i*(n_subdivisions+1)+j )*(n_subdivisions+1)+k +1 << ' '
+ << first_vertex_of_patch+((i+1)*(n_subdivisions+1)+j )*(n_subdivisions+1)+k +1 << ' '
+ << first_vertex_of_patch+((i+1)*(n_subdivisions+1)+j+1)*(n_subdivisions+1)+k +1 << ' '
+ << first_vertex_of_patch+(i*(n_subdivisions+1)+j+1 )*(n_subdivisions+1)+k +1 << ' '
+ << first_vertex_of_patch+(i*(n_subdivisions+1)+j )*(n_subdivisions+1)+k+1+1 << ' '
+ << first_vertex_of_patch+((i+1)*(n_subdivisions+1)+j )*(n_subdivisions+1)+k+1+1 << ' '
+ << first_vertex_of_patch+((i+1)*(n_subdivisions+1)+j+1)*(n_subdivisions+1)+k+1+1 << ' '
+ << first_vertex_of_patch+(i*(n_subdivisions+1)+j+1 )*(n_subdivisions+1)+k+1+1 << ' '
+ << endl;
+ ++present_cell;
+ };
+ break;
+ };
+
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+
+// finally update the number
+ // of the first vertex of this patch
+ switch (dim)
+ {
+ case 1:
+ first_vertex_of_patch += n_subdivisions+1;
+ break;
+ case 2:
+ first_vertex_of_patch += (n_subdivisions+1) *
+ (n_subdivisions+1);
+ break;
+ case 3:
+ first_vertex_of_patch += (n_subdivisions+1) *
+ (n_subdivisions+1) *
+ (n_subdivisions+1);
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
+ out << endl;
+
+ // note that we number starting with 1!
+ Assert (present_cell == n_cells+1,
+ ExcInternalError());
+ };
+
+
+/////////////////////////////
+ // now write data
+ if (n_data_sets != 0)
+ {
+ out << n_data_sets << " "; // number of vectors
+ for (unsigned int i=0; i<n_data_sets; ++i)
+ out << 1 << ' '; // number of components;
+ // only 1 supported presently
+ out << endl;
+
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << data_names[data_set]
+ << ",dimensionless" // no units supported at present
+ << endl;
+
+
+// loop over all patches
+ unsigned int present_node = 1;
+ for (typename vector<Patch<dim> >::const_iterator patch=patches.begin();
+ patch != patches.end(); ++patch)
+ {
+ const unsigned int n_subdivisions = patch->n_subdivisions;
+
+ Assert (patch->data.m() == n_data_sets,
+ ExcUnexpectedNumberOfDatasets (patch->data.m(), n_data_sets));
+ Assert (patch->data.n() == (dim==1 ?
+ n_subdivisions+1 :
+ (dim==2 ?
+ (n_subdivisions+1)*(n_subdivisions+1) :
+ (dim==3 ?
+ (n_subdivisions+1)*(n_subdivisions+1)*(n_subdivisions+1) :
+ 0))),
+ ExcInvalidDatasetSize (patch->data.n(), n_subdivisions+1));
+
+ switch (dim)
+ {
+ case 1:
+ {
+ for (unsigned int i=0; i<n_subdivisions+1; ++i, ++present_node)
+ {
+ out << present_node
+ << " ";
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << patch->data(data_set,i) << ' ';
+
+ out << endl;
+ };
+
+ break;
+ };
+
+ case 2:
+ {
+ for (unsigned int i=0; i<n_subdivisions+1; ++i)
+ for (unsigned int j=0; j<n_subdivisions+1; ++j)
+ {
+ out << present_node
+ << " ";
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << patch->data(data_set,i*(n_subdivisions+1) + j) << ' ';
+
+ out << endl;
+
+ ++present_node;
+ };
+
+ break;
+ };
+
+ case 3:
+ {
+ for (unsigned int i=0; i<n_subdivisions+1; ++i)
+ for (unsigned int j=0; j<n_subdivisions+1; ++j)
+ for (unsigned int k=0; k<n_subdivisions+1; ++k)
+ {
+ out << present_node
+ << " ";
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << patch->data(data_set,
+ (i*(n_subdivisions+1)+j)*(n_subdivisions+1)+k)
+ << ' ';
+
+ out << endl;
+
+ ++present_node;
+ };
+
+ break;
+ };
+
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
+ };
+
+ // no model data
+
+ // assert the stream is still ok
+ AssertThrow (out, ExcIO());
+};
+
+
+template <int dim>
+void DataOutBase::write_gnuplot (const vector<Patch<dim> > &patches,
+ const vector<string> &data_names,
+ const GnuplotFlags &/*flags*/,
+ ostream &out)
+{
+ AssertThrow (out, ExcIO());
+
+ Assert (patches.size() > 0, ExcNoPatches());
+
+ const unsigned int n_data_sets = data_names.size();
+
+ // write preamble
+ if (true)
+ {
+ // block this to have local
+ // variables destroyed after
+ // use
+ const time_t time1= time (0);
+ const 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<data_names.size(); ++i)
+ out << '<'
+ << data_names[i]
+ << "> ";
+ out << endl;
+ };
+
+
+// loop over all patches
+ for (typename vector<Patch<dim> >::const_iterator patch=patches.begin();
+ patch != patches.end(); ++patch)
+ {
+ const unsigned int n_subdivisions = patch->n_subdivisions;
+
+ Assert (patch->data.m() == n_data_sets,
+ ExcUnexpectedNumberOfDatasets (patch->data.m(), n_data_sets));
+ Assert (patch->data.n() == (dim==1 ?
+ n_subdivisions+1 :
+ (dim==2 ?
+ (n_subdivisions+1)*(n_subdivisions+1) :
+ (dim==3 ?
+ (n_subdivisions+1)*(n_subdivisions+1)*(n_subdivisions+1) :
+ 0))),
+ ExcInvalidDatasetSize (patch->data.n(), n_subdivisions+1));
+
+ switch (dim)
+ {
+ case 1:
+ {
+ for (unsigned int i=0; i<n_subdivisions+1; ++i)
+ {
+ // compute coordinates for
+ // this patch point
+ out << ((patch->vertices[1] * i / n_subdivisions) +
+ (patch->vertices[0] * (n_subdivisions-i) / n_subdivisions))
+ << ' ';
+
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << patch->data(data_set,i) << ' ';
+
+ out << endl;
+ };
+
+ // end of patch
+ out << endl
+ << endl;
+
+ break;
+ };
+
+ case 2:
+ {
+ for (unsigned int i=0; i<n_subdivisions+1; ++i)
+ {
+ for (unsigned int j=0; j<n_subdivisions+1; ++j)
+ {
+ const double x_frac = i * 1./n_subdivisions,
+ y_frac = j * 1./n_subdivisions;
+
+ // compute coordinates for
+ // this patch point
+ out << (((patch->vertices[1] * x_frac) +
+ (patch->vertices[0] * (1-x_frac))) * (1-y_frac) +
+ ((patch->vertices[2] * x_frac) +
+ (patch->vertices[3] * (1-x_frac))) * y_frac)
+ << ' ';
+
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << patch->data(data_set,i*(n_subdivisions+1) + j) << ' ';
+
+ out << endl;
+ };
+
+ // end of row in patch
+ out << endl;
+ };
+
+ // end of patch
+ out << endl;
+
+ break;
+ };
+
+ case 3:
+ {
+ // for all grid points: draw
+ // lines into all positive
+ // coordinate directions if
+ // there is another grid point
+ // there
+ for (unsigned int i=0; i<n_subdivisions+1; ++i)
+ for (unsigned int j=0; j<n_subdivisions+1; ++j)
+ for (unsigned int k=0; k<n_subdivisions+1; ++k)
+ {
+ // note the broken
+ // design of hexahedra
+ // in deal.II, where
+ // first the z-component
+ // is counted up, before
+ // increasing the y-
+ // coordinate
+ const double x_frac = i * 1./n_subdivisions,
+ y_frac = k * 1./n_subdivisions,
+ z_frac = j * 1./n_subdivisions;
+
+ // compute coordinates for
+ // this patch point
+ const Point<dim> this_point
+ = ((((patch->vertices[1] * x_frac) +
+ (patch->vertices[0] * (1-x_frac))) * (1-y_frac) +
+ ((patch->vertices[2] * x_frac) +
+ (patch->vertices[3] * (1-x_frac))) * y_frac) * (1-z_frac) +
+ (((patch->vertices[5] * x_frac) +
+ (patch->vertices[4] * (1-x_frac))) * (1-y_frac) +
+ ((patch->vertices[6] * x_frac) +
+ (patch->vertices[7] * (1-x_frac))) * y_frac) * z_frac);
+
+ // line into positive x-direction
+ // if possible
+ if (i < n_subdivisions)
+ {
+ // write point here
+ // and its data
+ out << this_point;
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << ' '
+ << patch->data(data_set,
+ (i*(n_subdivisions+1) + j)*(n_subdivisions+1)+k);
+ out << endl;
+
+ // write point there
+ // and its data
+ const double x_frac_new = x_frac + 1./n_subdivisions;
+ out << ((((patch->vertices[1] * x_frac_new) +
+ (patch->vertices[0] * (1-x_frac_new))) * (1-y_frac) +
+ ((patch->vertices[2] * x_frac_new) +
+ (patch->vertices[3] * (1-x_frac_new))) * y_frac) * (1-z_frac) +
+ (((patch->vertices[5] * x_frac_new) +
+ (patch->vertices[4] * (1-x_frac_new))) * (1-y_frac) +
+ ((patch->vertices[6] * x_frac_new) +
+ (patch->vertices[7] * (1-x_frac_new))) * y_frac) * z_frac);
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << ' '
+ << patch->data(data_set,
+ ((i+1)*(n_subdivisions+1) + j)*(n_subdivisions+1)+k);
+ out << endl;
+
+ // end of line
+ out << endl
+ << endl;
+ };
+
+ // line into positive y-direction
+ // if possible
+ if (j < n_subdivisions)
+ {
+ // write point here
+ // and its data
+ out << this_point;
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << ' '
+ << patch->data(data_set,
+ (i*(n_subdivisions+1) + j)*(n_subdivisions+1)+k);
+ out << endl;
+
+ // write point there
+ // and its data
+ const double z_frac_new = z_frac + 1./n_subdivisions;
+ out << ((((patch->vertices[1] * x_frac) +
+ (patch->vertices[0] * (1-x_frac))) * (1-y_frac) +
+ ((patch->vertices[2] * x_frac) +
+ (patch->vertices[3] * (1-x_frac))) * y_frac) * (1-z_frac_new) +
+ (((patch->vertices[5] * x_frac) +
+ (patch->vertices[4] * (1-x_frac))) * (1-y_frac) +
+ ((patch->vertices[6] * x_frac) +
+ (patch->vertices[7] * (1-x_frac))) * y_frac) * z_frac_new);
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << ' '
+ << patch->data(data_set,
+ (i*(n_subdivisions+1) + (j+1))*(n_subdivisions+1)+k);
+ out << endl;
+
+ // end of line
+ out << endl
+ << endl;
+ };
+
+ // line into positive z-direction
+ // if possible
+ if (k < n_subdivisions)
+ {
+ // write point here
+ // and its data
+ out << this_point;
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << ' '
+ << patch->data(data_set,
+ (i*(n_subdivisions+1) + j)*(n_subdivisions+1)+k);
+ out << endl;
+
+ // write point there
+ // and its data
+ const double y_frac_new = y_frac + 1./n_subdivisions;
+ out << ((((patch->vertices[1] * x_frac) +
+ (patch->vertices[0] * (1-x_frac))) * (1-y_frac_new) +
+ ((patch->vertices[2] * x_frac) +
+ (patch->vertices[3] * (1-x_frac))) * y_frac_new) * (1-z_frac) +
+ (((patch->vertices[5] * x_frac) +
+ (patch->vertices[4] * (1-x_frac))) * (1-y_frac_new) +
+ ((patch->vertices[6] * x_frac) +
+ (patch->vertices[7] * (1-x_frac))) * y_frac_new) * z_frac);
+ for (unsigned int data_set=0; data_set<n_data_sets; ++data_set)
+ out << ' '
+ << patch->data(data_set,
+ (i*(n_subdivisions+1) + j)*(n_subdivisions+1)+k+1);
+ out << endl;
+
+ // end of line
+ out << endl
+ << endl;
+ };
+
+ };
+
+ break;
+ };
+
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
+
+ AssertThrow (out, ExcIO());
+};
+
+
+template <int dim>
+void DataOutBase::write_povray (const vector<Patch<dim> > &patches,
+ const vector<string> &data_names,
+ const PovrayFlags &flags,
+ ostream &out)
+{
+ AssertThrow (out, ExcIO());
+
+ Assert (patches.size() > 0, ExcNoPatches());
+ Assert (dim==2, ExcNotImplemented()); // only for 2-D
+
+ const unsigned int n_data_sets = data_names.size();
+
+ // write preamble
+ if (true)
+ {
+ // block this to have local
+ // variables destroyed after use
+ const time_t time1= time (0);
+ const 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 POVRAY format see the POVRAY manual."
+ << endl
+ << "*/ " << endl;
+
+ // include files
+ out << "#include \"colors.inc\" " << endl
+ << "#include \"textures.inc\" " << endl;
+
+
+// use external include file for textures,
+ // camera and light
+ if (flags.external_data)
+ out << "#include \"data.inc\" " << endl;
else // all definitions in data file
{
// camera