From 4f8944c3c084e719c185d5c2abd0e408c65fd967 Mon Sep 17 00:00:00 2001 From: bangerth Date: Tue, 9 Nov 2010 17:12:56 +0000 Subject: [PATCH] One more file moved to compile-once mode. git-svn-id: https://svn.dealii.org/trunk@22651 0785d39b-7218-0410-832d-ea1e28bc413d --- .../source/grid/grid_out.all_dimensions.cc | 1907 ++++++++++++++++ deal.II/source/grid/grid_out.cc | 1957 ----------------- deal.II/source/grid/grid_out.inst.in | 60 + 3 files changed, 1967 insertions(+), 1957 deletions(-) delete mode 100644 deal.II/source/grid/grid_out.cc create mode 100644 deal.II/source/grid/grid_out.inst.in diff --git a/deal.II/source/grid/grid_out.all_dimensions.cc b/deal.II/source/grid/grid_out.all_dimensions.cc index 880c5d361f..d95514c166 100644 --- a/deal.II/source/grid/grid_out.all_dimensions.cc +++ b/deal.II/source/grid/grid_out.all_dimensions.cc @@ -14,6 +14,22 @@ #include #include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + 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 +void GridOut::write_dx (const Triangulation &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 > &vertices = tria.get_vertices(); + const std::vector &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 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::active_cell_iterator cell; + const typename Triangulation::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; i1) ? dx_flags.write_faces : false; + + const unsigned int n_cells = tria.n_active_cells(); + const unsigned int n_faces = tria.n_active_cells() + * GeometryInfo::faces_per_cell; + + const unsigned int n_vertices_per_cell = GeometryInfo::vertices_per_cell; + const unsigned int n_vertices_per_face = GeometryInfo::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::vertices_per_cell; ++v) + out << '\t' << renumber[cell->vertex_index(GeometryInfo::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::faces_per_cell;++f) + { + typename Triangulation::face_iterator face = cell->face(f); + + for (unsigned int v=0; v::vertices_per_face; ++v) + out << '\t' << renumber[face->vertex_index(GeometryInfo::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::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::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::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 +void GridOut::write_msh (const Triangulation &tria, + std::ostream &out) const +{ + AssertThrow (out, ExcIO()); + + // get the positions of the + // vertices and whether they are + // used. + const std::vector > &vertices = tria.get_vertices(); + const std::vector &vertex_used = tria.get_used_vertices(); + + const unsigned int n_vertices = tria.n_used_vertices(); + + typename Triangulation::active_cell_iterator cell=tria.begin_active(); + const typename Triangulation::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(cell->material_id()) << ' ' + << cell->subdomain_id() << ' ' + << GeometryInfo::vertices_per_cell << ' '; + + // Vertex numbering follows UCD conventions. + + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + out << cell->vertex_index(GeometryInfo::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 +void GridOut::write_ucd (const Triangulation &tria, + std::ostream &out) const +{ + AssertThrow (out, ExcIO()); + + // get the positions of the + // vertices and whether they are + // used. + const std::vector > &vertices = tria.get_vertices(); + const std::vector &vertex_used = tria.get_used_vertices(); + + const unsigned int n_vertices = tria.n_used_vertices(); + + typename Triangulation::active_cell_iterator cell=tria.begin_active(); + const typename Triangulation::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(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::vertices_per_cell; + ++vertex) + out << cell->vertex_index(GeometryInfo::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 +void GridOut::write_xfig ( + const Triangulation&, + std::ostream&, + const Mapping*) 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::vertices_per_cell; + const unsigned int nf = GeometryInfo::faces_per_cell; + const unsigned int nvf = GeometryInfo::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::cell_iterator cell = tria.begin(); + const typename Triangulation::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& p = cell->vertex( + GeometryInfo::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::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& 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()); +} + + + + + +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 +unsigned int GridOut::n_boundary_faces (const Triangulation &tria) const +{ + typename Triangulation::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 +unsigned int GridOut::n_boundary_lines (const Triangulation &tria) const +{ + typename Triangulation::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 +void GridOut::write_msh_faces (const Triangulation &tria, + const unsigned int starting_index, + std::ostream &out) const +{ + typename Triangulation::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(face->boundary_indicator()) + << ' ' + << static_cast(face->boundary_indicator()) + << ' ' << GeometryInfo::vertices_per_face; + // note: vertex numbers are 1-base + for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) + out << ' ' + << face->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1; + out << '\n'; + + ++index; + }; +} + + +template +void GridOut::write_msh_lines (const Triangulation &tria, + const unsigned int starting_index, + std::ostream &out) const +{ + typename Triangulation::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(line->boundary_indicator()) + << ' ' + << static_cast(line->boundary_indicator()) + << " 2 "; + // note: vertex numbers are 1-base + for (unsigned int vertex=0; vertex<2; ++vertex) + out << ' ' + << line->vertex_index(GeometryInfo::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 +void GridOut::write_ucd_faces (const Triangulation &tria, + const unsigned int starting_index, + std::ostream &out) const +{ + typename Triangulation::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(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::vertices_per_face; ++vertex) + out << face->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1 << ' '; + out << '\n'; + + ++index; + }; +} + + + +template +void GridOut::write_ucd_lines (const Triangulation &tria, + const unsigned int starting_index, + std::ostream &out) const +{ + typename Triangulation::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(line->boundary_indicator()) + << " line "; + // note: vertex numbers are 1-base + for (unsigned int vertex=0; vertex<2; ++vertex) + out << line->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1 << ' '; + out << '\n'; + + ++index; + }; +} + + +namespace internal +{ + namespace + { + template + 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::active_cell_iterator + cell=tria.begin_active(); + const typename dealii::Triangulation::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(cell->material_id()) << '\n' + << cell->vertex(1) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << '\n'; + } + + // make sure everything now gets to + // disk + out.flush (); + + AssertThrow (out, ExcIO()); + } + + + + template + 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::active_cell_iterator + cell=tria.begin_active(); + const typename dealii::Triangulation::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 *q_projector=0; + std::vector > 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 dummy_weights(n_points, 1./n_points); + Quadrature quadrature(boundary_points, dummy_weights); + + q_projector = new Quadrature (QProjector::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::vertices_per_cell; ++i) + out << cell->vertex(GeometryInfo::ucd_to_deal[i]) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n'; + out << cell->vertex(0) + << ' ' << cell->level() + << ' ' << static_cast(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::faces_per_cell; ++face_no) + { + const typename dealii::Triangulation::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; itransform_unit_to_real_cell + (cell, q_projector->point(offset+i))) + << ' ' << cell->level() + << ' ' << static_cast(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(cell->material_id()) + << '\n' + << face->vertex(1) + << ' ' << cell->level() + << ' ' << static_cast(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 + 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::active_cell_iterator + cell=tria.begin_active(); + const typename dealii::Triangulation::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 *q_projector=0; + std::vector > 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 dummy_weights(n_points, 1./n_points); + Quadrature<1> quadrature1d(boundary_points, dummy_weights); + + // tensor product of points, + // only one copy + QIterated quadrature(quadrature1d, 1); + q_projector = new Quadrature (QProjector::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(cell->material_id()) << '\n' + << cell->vertex(1) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << cell->vertex(5) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << cell->vertex(4) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << cell->vertex(0) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << '\n'; + // back face + out << cell->vertex(2) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << cell->vertex(3) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << cell->vertex(7) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << cell->vertex(6) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << cell->vertex(2) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << '\n'; + + // now for the four connecting lines + out << cell->vertex(0) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << cell->vertex(2) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << '\n'; + out << cell->vertex(1) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << cell->vertex(3) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << '\n'; + out << cell->vertex(5) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << cell->vertex(7) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << '\n'; + out << cell->vertex(4) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << cell->vertex(6) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << '\n'; + } + else + { + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + const typename dealii::Triangulation::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 p0=mapping->transform_unit_to_real_cell( + cell, q_projector->point(offset+i*n_points+j)); + out << p0 + << ' ' << cell->level() + << ' ' << static_cast(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(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(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(cell->material_id()) << '\n'; + // and the + // first + // point + // again + out << p0 + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n'; + out << '\n' << '\n'; + } + } + else + { + for (unsigned int l=0; l::lines_per_face; ++l) + { + const typename dealii::Triangulation::line_iterator + line=face->line(l); + + const Point &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::project_to_line + // which is + // not yet + // implemented + const Point u0=mapping->transform_real_to_unit_cell(cell, v0), + u1=mapping->transform_real_to_unit_cell(cell, v1); + + for (unsigned int i=0; itransform_unit_to_real_cell + (cell, (1-boundary_points[i][0])*u0+boundary_points[i][0]*u1)) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n'; + } + else + out << v0 + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << '\n' + << v1 + << ' ' << cell->level() + << ' ' << static_cast(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 +void GridOut::write_gnuplot ( + const Triangulation &tria, + std::ostream &out, + const Mapping *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 + void write_eps (const dealii::Triangulation &tria, + std::ostream &out, + const Mapping *mapping, + const GridOutFlags::Eps<2> &eps_flags_2, + const GridOutFlags::Eps<3> &eps_flags_3) + { + typedef std::list 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(eps_flags_2) : + (dim==3 ? + static_cast(eps_flags_3) : + *static_cast(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::active_cell_iterator + cell=tria.begin_active(), + endc=tria.end(); + + for (; cell!=endc; ++cell) + for (unsigned int line_no=0; + line_no::lines_per_cell; ++line_no) + { + typename dealii::Triangulation::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, but we + // want a Point<2>. + // in fact, whenever + // we're here, the + // vertex is a + // Point, 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 > boundary_points (n_points); + + for (unsigned int i=0; i quadrature (boundary_points); + Quadrature q_projector (QProjector::project_to_all_faces(quadrature)); + + // next loop over all + // boundary faces and + // generate the info from + // them + typename dealii::Triangulation::active_cell_iterator + cell=tria.begin_active (); + const typename dealii::Triangulation::active_cell_iterator + end=tria.end (); + for (; cell!=end; ++cell) + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + const typename dealii::Triangulation::face_iterator + face = cell->face(face_no); + + if (face->at_boundary()) + { + Point 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 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 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::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 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 vector1 + = Point(0,0,1) - ((Point(0,0,1) * view_direction) * view_direction); + const Point 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 vector2 + = (Point(1,0,0) + - ((Point(1,0,0) * view_direction) * view_direction) + - ((Point(1,0,0) * unit_vector1) * unit_vector1)); + const Point unit_vector2 = vector2 / std::sqrt(vector2.square()); + + + for (; cell!=endc; ++cell) + for (unsigned int line_no=0; + line_no::lines_per_cell; ++line_no) + { + typename dealii::Triangulation::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(std::floor(( (x_max-x_min) * scale )+1)) + << ' ' + << static_cast(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::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 treated_vertices; + typename dealii::Triangulation::active_cell_iterator + cell = tria.begin_active (), + endc = tria.end (); + for (; cell!=endc; ++cell) + for (unsigned int vertex=0; + vertex::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 +void GridOut::write_eps (const Triangulation &tria, + std::ostream &out, + const Mapping *mapping) const +{ + internal::write_eps (tria, out, mapping, + eps_flags_2, eps_flags_3); +} + + +template +void GridOut::write (const Triangulation &tria, + std::ostream &out, + const OutputFormat output_format, + const Mapping *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 +void GridOut::write ( + const Triangulation &tria, + std::ostream &out, + const Mapping *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 index e843f634c6..0000000000 --- a/deal.II/source/grid/grid_out.cc +++ /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 -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -DEAL_II_NAMESPACE_OPEN - - -#if deal_II_dimension == 1 - - -template -void GridOut::write_dx (const Triangulation &, - std::ostream &) const -{ - Assert (false, ExcNotImplemented()); -} - -#else - - -template -void GridOut::write_dx (const Triangulation &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 > &vertices = tria.get_vertices(); - const std::vector &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 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::active_cell_iterator cell; - const typename Triangulation::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; i1) ? dx_flags.write_faces : false; - - const unsigned int n_cells = tria.n_active_cells(); - const unsigned int n_faces = tria.n_active_cells() - * GeometryInfo::faces_per_cell; - - const unsigned int n_vertices_per_cell = GeometryInfo::vertices_per_cell; - const unsigned int n_vertices_per_face = GeometryInfo::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::vertices_per_cell; ++v) - out << '\t' << renumber[cell->vertex_index(GeometryInfo::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::faces_per_cell;++f) - { - typename Triangulation::face_iterator face = cell->face(f); - - for (unsigned int v=0; v::vertices_per_face; ++v) - out << '\t' << renumber[face->vertex_index(GeometryInfo::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::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::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::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 -void GridOut::write_msh (const Triangulation &tria, - std::ostream &out) const -{ - AssertThrow (out, ExcIO()); - - // get the positions of the - // vertices and whether they are - // used. - const std::vector > &vertices = tria.get_vertices(); - const std::vector &vertex_used = tria.get_used_vertices(); - - const unsigned int n_vertices = tria.n_used_vertices(); - - typename Triangulation::active_cell_iterator cell=tria.begin_active(); - const typename Triangulation::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(cell->material_id()) << ' ' - << cell->subdomain_id() << ' ' - << GeometryInfo::vertices_per_cell << ' '; - - // Vertex numbering follows UCD conventions. - - for (unsigned int vertex=0; vertex::vertices_per_cell; - ++vertex) - out << cell->vertex_index(GeometryInfo::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 -void GridOut::write_ucd (const Triangulation &tria, - std::ostream &out) const -{ - AssertThrow (out, ExcIO()); - - // get the positions of the - // vertices and whether they are - // used. - const std::vector > &vertices = tria.get_vertices(); - const std::vector &vertex_used = tria.get_used_vertices(); - - const unsigned int n_vertices = tria.n_used_vertices(); - - typename Triangulation::active_cell_iterator cell=tria.begin_active(); - const typename Triangulation::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(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::vertices_per_cell; - ++vertex) - out << cell->vertex_index(GeometryInfo::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 -void GridOut::write_xfig ( - const Triangulation&, - std::ostream&, - const Mapping*) const -{ - Assert (false, ExcNotImplemented()); -} - -#else - -//TODO:[GK] Obey parameters -//TODO:[GK] Flip y-axis? -template -void GridOut::write_xfig ( - const Triangulation& tria, - std::ostream& out, - const Mapping* /*mapping*/) const -{ - const unsigned int nv = GeometryInfo::vertices_per_cell; - const unsigned int nf = GeometryInfo::faces_per_cell; - const unsigned int nvf = GeometryInfo::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::cell_iterator cell = tria.begin(); - const typename Triangulation::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& p = cell->vertex( - GeometryInfo::ucd_to_deal[k % nv]); - for (unsigned int d=0;d::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& p = face->vertex(k % nv); - for (unsigned int d=0;d &) 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 -unsigned int GridOut::n_boundary_faces (const Triangulation &tria) const -{ - typename Triangulation::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 -unsigned int GridOut::n_boundary_lines (const Triangulation &tria) const -{ - typename Triangulation::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 -void GridOut::write_msh_faces (const Triangulation &tria, - const unsigned int starting_index, - std::ostream &out) const -{ - typename Triangulation::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(face->boundary_indicator()) - << ' ' - << static_cast(face->boundary_indicator()) - << ' ' << GeometryInfo::vertices_per_face; - // note: vertex numbers are 1-base - for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) - out << ' ' - << face->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1; - out << '\n'; - - ++index; - }; -} - - -template -void GridOut::write_msh_lines (const Triangulation &tria, - const unsigned int starting_index, - std::ostream &out) const -{ - typename Triangulation::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(line->boundary_indicator()) - << ' ' - << static_cast(line->boundary_indicator()) - << " 2 "; - // note: vertex numbers are 1-base - for (unsigned int vertex=0; vertex<2; ++vertex) - out << ' ' - << line->vertex_index(GeometryInfo::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 -void GridOut::write_ucd_faces (const Triangulation &tria, - const unsigned int starting_index, - std::ostream &out) const -{ - typename Triangulation::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(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::vertices_per_face; ++vertex) - out << face->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1 << ' '; - out << '\n'; - - ++index; - }; -} - - - -template -void GridOut::write_ucd_lines (const Triangulation &tria, - const unsigned int starting_index, - std::ostream &out) const -{ - typename Triangulation::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(line->boundary_indicator()) - << " line "; - // note: vertex numbers are 1-base - for (unsigned int vertex=0; vertex<2; ++vertex) - out << line->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1 << ' '; - out << '\n'; - - ++index; - }; -} - - -#if deal_II_dimension==1 - -template -void GridOut::write_gnuplot ( - const Triangulation &tria, - std::ostream &out, - const Mapping *) const -{ - AssertThrow (out, ExcIO()); - - typename Triangulation::active_cell_iterator cell=tria.begin_active(); - const typename Triangulation::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(cell->material_id()) << '\n' - << cell->vertex(1) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << '\n'; - } - - // make sure everything now gets to - // disk - out.flush (); - - AssertThrow (out, ExcIO()); -} - -#endif - - -#if deal_II_dimension==2 - -template -void GridOut::write_gnuplot ( - const Triangulation &tria, - std::ostream &out, - const Mapping *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::active_cell_iterator cell=tria.begin_active(); - const typename Triangulation::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 *q_projector=0; - std::vector > 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 dummy_weights(n_points, 1./n_points); - Quadrature quadrature(boundary_points, dummy_weights); - - q_projector = new Quadrature (QProjector::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::vertices_per_cell; ++i) - out << cell->vertex(GeometryInfo::ucd_to_deal[i]) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n'; - out << cell->vertex(0) - << ' ' << cell->level() - << ' ' << static_cast(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::faces_per_cell; ++face_no) - { - const typename Triangulation::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; itransform_unit_to_real_cell - (cell, q_projector->point(offset+i))) - << ' ' << cell->level() - << ' ' << static_cast(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(cell->material_id()) - << '\n' - << face->vertex(1) - << ' ' << cell->level() - << ' ' << static_cast(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 -void GridOut::write_gnuplot (const Triangulation &tria, - std::ostream &out, - const Mapping *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::active_cell_iterator cell=tria.begin_active(); - const typename Triangulation::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 *q_projector=0; - std::vector > 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 dummy_weights(n_points, 1./n_points); - Quadrature<1> quadrature1d(boundary_points, dummy_weights); - - // tensor product of points, - // only one copy - QIterated quadrature(quadrature1d, 1); - q_projector = new Quadrature (QProjector::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(cell->material_id()) << '\n' - << cell->vertex(1) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << cell->vertex(5) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << cell->vertex(4) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << cell->vertex(0) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << '\n'; - // back face - out << cell->vertex(2) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << cell->vertex(3) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << cell->vertex(7) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << cell->vertex(6) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << cell->vertex(2) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << '\n'; - - // now for the four connecting lines - out << cell->vertex(0) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << cell->vertex(2) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << '\n'; - out << cell->vertex(1) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << cell->vertex(3) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << '\n'; - out << cell->vertex(5) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << cell->vertex(7) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << '\n'; - out << cell->vertex(4) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << cell->vertex(6) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << '\n'; - } - else - { - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - { - const typename Triangulation::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 p0=mapping->transform_unit_to_real_cell( - cell, q_projector->point(offset+i*n_points+j)); - out << p0 - << ' ' << cell->level() - << ' ' << static_cast(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(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(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(cell->material_id()) << '\n'; - // and the - // first - // point - // again - out << p0 - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n'; - out << '\n' << '\n'; - } - } - else - { - for (unsigned int l=0; l::lines_per_face; ++l) - { - const typename Triangulation::line_iterator - line=face->line(l); - - const Point &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::project_to_line - // which is - // not yet - // implemented - const Point u0=mapping->transform_real_to_unit_cell(cell, v0), - u1=mapping->transform_real_to_unit_cell(cell, v1); - - for (unsigned int i=0; itransform_unit_to_real_cell - (cell, (1-boundary_points[i][0])*u0+boundary_points[i][0]*u1)) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n'; - } - else - out << v0 - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << v1 - << ' ' << cell->level() - << ' ' << static_cast(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 -void GridOut::write_eps ( - const Triangulation &, - std::ostream &, - const Mapping *) const -{ - Assert(false, ExcNotImplemented()); -} - - -#else - -template -void GridOut::write_eps ( - const Triangulation &tria, - std::ostream &out, - const Mapping *mapping) const -{ - - typedef std::list 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(eps_flags_2) : - (dim==3 ? - static_cast(eps_flags_3) : - *static_cast(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::active_cell_iterator - cell=tria.begin_active(), - endc=tria.end(); - - for (; cell!=endc; ++cell) - for (unsigned int line_no=0; - line_no::lines_per_cell; ++line_no) - { - typename Triangulation::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, but we - // want a Point<2>. - // in fact, whenever - // we're here, the - // vertex is a - // Point, 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 > boundary_points (n_points); - - for (unsigned int i=0; i quadrature (boundary_points); - Quadrature q_projector (QProjector::project_to_all_faces(quadrature)); - - // next loop over all - // boundary faces and - // generate the info from - // them - typename Triangulation::active_cell_iterator cell=tria.begin_active (); - const typename Triangulation::active_cell_iterator end=tria.end (); - for (; cell!=end; ++cell) - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - { - const typename Triangulation::face_iterator - face = cell->face(face_no); - - if (face->at_boundary()) - { - Point 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 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 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::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 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 vector1 - = Point(0,0,1) - ((Point(0,0,1) * view_direction) * view_direction); - const Point 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 vector2 - = (Point(1,0,0) - - ((Point(1,0,0) * view_direction) * view_direction) - - ((Point(1,0,0) * unit_vector1) * unit_vector1)); - const Point unit_vector2 = vector2 / std::sqrt(vector2.square()); - - - for (; cell!=endc; ++cell) - for (unsigned int line_no=0; - line_no::lines_per_cell; ++line_no) - { - typename Triangulation::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(std::floor(( (x_max-x_min) * scale )+1)) - << ' ' - << static_cast(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::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 treated_vertices; - typename Triangulation::active_cell_iterator - cell = tria.begin_active (), - endc = tria.end (); - for (; cell!=endc; ++cell) - for (unsigned int vertex=0; - vertex::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 -void GridOut::write ( - const Triangulation &tria, - std::ostream &out, - const OutputFormat output_format, - const Mapping *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 -void GridOut::write ( - const Triangulation &tria, - std::ostream &out, - const Mapping *mapping) const -{ - write(tria, out, default_format, mapping); -} - - -//TODO: Should all write functions be instantiated? - -// explicit instantiations -template void GridOut::write_dx -(const Triangulation&, - std::ostream&) const; -template void GridOut::write_msh -(const Triangulation&, - std::ostream&) const; -template void GridOut::write_xfig -(const Triangulation&, - std::ostream&, - const Mapping*) const; -template void GridOut::write_gnuplot -(const Triangulation&, - std::ostream&, - const Mapping*) const; -template void GridOut::write_ucd -(const Triangulation &, - std::ostream &) const; -template void GridOut::write_eps -(const Triangulation &, - std::ostream &, - const Mapping *) const; -template void GridOut::write -(const Triangulation &, - std::ostream &, const OutputFormat, - const Mapping *) const; -template void GridOut::write -(const Triangulation &, - std::ostream &, const Mapping *) const; - -#if deal_II_dimension == 1 || deal_II_dimension ==2 -// explicit instantiations for codimension one -template void GridOut::write_msh -(const Triangulation&, - std::ostream&) const; -template void GridOut::write_ucd -(const Triangulation &, - std::ostream &) const; -template void GridOut::write_gnuplot -(const Triangulation&, - std::ostream&, - const Mapping*) 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 index 0000000000..70af179131 --- /dev/null +++ b/deal.II/source/grid/grid_out.inst.in @@ -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&, + std::ostream&) const; + template void GridOut::write_msh + (const Triangulation&, + std::ostream&) const; + template void GridOut::write_xfig + (const Triangulation&, + std::ostream&, + const Mapping*) const; + template void GridOut::write_gnuplot + (const Triangulation&, + std::ostream&, + const Mapping*) const; + template void GridOut::write_ucd + (const Triangulation &, + std::ostream &) const; + template void GridOut::write_eps + (const Triangulation &, + std::ostream &, + const Mapping *) const; + template void GridOut::write + (const Triangulation &, + std::ostream &, const OutputFormat, + const Mapping *) const; + template void GridOut::write + (const Triangulation &, + std::ostream &, const Mapping *) const; + +#if deal_II_dimension == 1 || deal_II_dimension ==2 +// explicit instantiations for codimension one + template void GridOut::write_msh + (const Triangulation&, + std::ostream&) const; + template void GridOut::write_ucd + (const Triangulation &, + std::ostream &) const; + template void GridOut::write_gnuplot + (const Triangulation&, + std::ostream&, + const Mapping*) const; +#endif + } + -- 2.39.5