template <int dim>
-void GridOut::write_ucd (const Triangulation<dim> &tria,
+void GridOut::write_msh (const Triangulation<dim> &tria,
std::ostream &out)
{
AssertThrow (out, ExcIO());
typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
const typename Triangulation<dim>::active_cell_iterator endc=tria.end();
- // 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)
- << " 0 0 0" // no data
- << '\n';
+ // Write Header
+ // The file format is:
+ /*
- // actually write the vertices.
- // note that we shall number them
- // with first index 1 instead of 0
+
+ $NOD
+ number-of-nodes
+ node-number x-coord y-coord z-coord
+ ...
+ $ENDNOD
+ $ELM
+ number-of-elements
+ elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list
+ ...
+ $ENDELM
+ */
+ out << "$NOD" << std::endl
+ << n_vertices << std::endl;
+
+ // actually write the vertices.
+ // note that we shall number them
+ // with first index 1 instead of 0
for (unsigned int i=0; i<vertices.size(); ++i)
- if (vertex_used[i])
+ if (vertex_used[i])
{
- out << i+1 // vertex index
- << " "
- << vertices[i];
- for (unsigned int d=dim+1; d<=3; ++d)
- out << " 0"; // fill with zeroes
- out << '\n';
+ out << i+1 // vertex index
+ << " "
+ << vertices[i];
+ for (unsigned int d=dim+1; d<=3; ++d)
+ out << " 0"; // fill with zeroes
+ out << std::endl;
};
-
- // write cells. Enumerate cells
- // consecutively, starting with 1
+
+ // Write cells preamble
+ out << "$ENDNOD" << std::endl
+ << "$ELM" << std::endl
+ << tria.n_active_cells() + (msh_flags.write_faces ?
+ n_boundary_faces(tria) : 0) << std::endl;
+
+ /*
+ elm-type
+ defines the geometrical type of the n-th element:
+ 1
+ Line (2 nodes).
+ 2
+ Triangle (3 nodes).
+ 3
+ Quadrangle (4 nodes).
+ 4
+ Tetrahedron (4 nodes).
+ 5
+ Hexahedron (8 nodes).
+ 6
+ Prism (6 nodes).
+ 7
+ Pyramid (5 nodes).
+ 8
+ Second order line (3 nodes: 2 associated with the vertices and 1 with the edge).
+ 9
+ Second order triangle (6 nodes: 3 associated with the vertices and 3 with the edges).
+ 10
+ Second order quadrangle (9 nodes: 4 associated with the vertices, 4 with the edges and 1 with the face).
+ 11
+ Second order tetrahedron (10 nodes: 4 associated with the vertices and 6 with the edges).
+ 12
+ Second order hexahedron (27 nodes: 8 associated with the vertices, 12 with the edges, 6 with the faces and 1 with the volume).
+ 13
+ Second order prism (18 nodes: 6 associated with the vertices, 9 with the edges and 3 with the quadrangular faces).
+ 14
+ Second order pyramid (14 nodes: 5 associated with the vertices, 8 with the edges and 1 with the quadrangular face).
+ 15
+ Point (1 node).
+ */
+ unsigned int elm_type;
+ switch(dim) {
+ case 1:
+ elm_type = 1;
+ break;
+ case 2:
+ elm_type = 3;
+ break;
+ case 3:
+ elm_type = 5;
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+
+ // write cells. Enumerate cells
+ // consecutively, starting with 1
unsigned int cell_index=1;
for (cell=tria.begin_active();
- cell!=endc; ++cell, ++cell_index)
+ cell!=endc; ++cell, ++cell_index)
+ {
+ out << cell_index << ' ' << elm_type << ' '
+ << static_cast<unsigned int>(cell->material_id()) << ' '
+ << cell->subdomain_id() << ' '
+ << GeometryInfo<dim>::vertices_per_cell << ' ';
+
+ // Vertex numbering follows UCD conventions.
+
+ for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex)
+ out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex])+1 << ' ';
+ out << std::endl;
+ };
+
+ // write faces with non-zero boundary
+ // indicator
+ if (msh_flags.write_faces)
+ write_msh_faces (tria, cell_index, out);
+
+ out << "$ENDELM" << std::endl;
+
+ // make sure everything now gets to
+ // disk
+ out.flush ();
+
+ AssertThrow (out, ExcIO());
+}
+
+ template <int dim>
+void GridOut::write_ucd (const Triangulation<dim> &tria,
+ std::ostream &out)
+{
+ AssertThrow (out, ExcIO());
+
+ // get the positions of the
+ // vertices and whether they are
+ // used.
+ const std::vector<Point<dim> > &vertices = tria.get_vertices();
+ const std::vector<bool> &vertex_used = tria.get_used_vertices();
+
+ const unsigned int n_vertices = tria.n_used_vertices();
+
+ typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
+ const typename Triangulation<dim>::active_cell_iterator endc=tria.end();
+
+ // write preamble
+ if (ucd_flags.write_preamble)
{
- out << cell_index << ' '
- << static_cast<unsigned int>(cell->material_id())
- << " ";
- switch (dim)
+ // 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)
+ << " 0 0 0" // no data
+ << '\n';
+
+ // actually write the vertices.
+ // note that we shall number them
+ // with first index 1 instead of 0
+ for (unsigned int i=0; i<vertices.size(); ++i)
+ if (vertex_used[i])
{
- case 1: out << "line "; break;
- case 2: out << "quad "; break;
- case 3: out << "hex "; break;
- default:
- Assert (false, ExcNotImplemented());
+ out << i+1 // vertex index
+ << " "
+ << vertices[i];
+ for (unsigned int d=dim+1; d<=3; ++d)
+ out << " 0"; // fill with zeroes
+ out << '\n';
};
- // it follows a list of the
- // vertices of each cell. in 1d
- // this is simply a list of the
- // two vertices, in 2d its counter
- // clockwise, as usual in this
- // library. in 3d, the same applies
- // (special thanks to AVS for
- // numbering their vertices in a
- // way compatible to deal.II!)
- //
- // technical reference:
- // AVS Developer's Guide, Release 4,
- // May, 1992, p. E6
- //
- // note: vertex numbers are 1-base
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
- ++vertex)
- out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex])+1 << ' ';
- out << '\n';
+ // write cells. Enumerate cells
+ // consecutively, starting with 1
+ unsigned int cell_index=1;
+ for (cell=tria.begin_active();
+ cell!=endc; ++cell, ++cell_index)
+ {
+ out << cell_index << ' '
+ << static_cast<unsigned int>(cell->material_id())
+ << " ";
+ switch (dim)
+ {
+ case 1: out << "line "; break;
+ case 2: out << "quad "; break;
+ case 3: out << "hex "; break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+ // it follows a list of the
+ // vertices of each cell. in 1d
+ // this is simply a list of the
+ // two vertices, in 2d its counter
+ // clockwise, as usual in this
+ // library. in 3d, the same applies
+ // (special thanks to AVS for
+ // numbering their vertices in a
+ // way compatible to deal.II!)
+ //
+ // technical reference:
+ // AVS Developer's Guide, Release 4,
+ // May, 1992, p. E6
+ //
+ // note: vertex numbers are 1-base
+ for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex)
+ out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex])+1 << ' ';
+ out << '\n';
};
- // write faces with non-zero boundary
- // indicator
- if (ucd_flags.write_faces)
- write_ucd_faces (tria, cell_index, out);
-
- // make sure everything now gets to
- // disk
- out.flush ();
+ // write faces with non-zero boundary
+ // indicator
+ if (ucd_flags.write_faces)
+ write_ucd_faces (tria, cell_index, out);
- AssertThrow (out, ExcIO());
+ // make sure everything now gets to
+ // disk
+ out.flush ();
+
+ AssertThrow (out, ExcIO());
}
#if deal_II_dimension != 2
-template <int dim>
+ template <int dim>
void GridOut::write_xfig (const Triangulation<dim>&,
- std::ostream&,
- const Mapping<dim>*)
+ std::ostream&,
+ const Mapping<dim>*)
{
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
}
#else
//TODO:[GK] Obey parameters
//TODO:[GK] Flip y-axis?
-template <int dim>
+ template <int dim>
void GridOut::write_xfig (const Triangulation<dim>& tria,
- std::ostream& out,
- const Mapping<dim>* /*mapping*/)
+ std::ostream& out,
+ const Mapping<dim>* /*mapping*/)
{
- const unsigned int nv = GeometryInfo<dim>::vertices_per_cell;
- const unsigned int nf = GeometryInfo<dim>::faces_per_cell;
- const unsigned int nvf = GeometryInfo<dim>::vertices_per_face;
-
- // The following text was copied
- // from an existing XFig file.
- out << "#FIG 3.2\nLandscape\nCenter\nInches" << '\n'
- << "A4\n100.00\nSingle" << '\n'
- // Background is transparent
- << "-3" << '\n'
- << "# generated by deal.II GridOut class" << '\n'
- << "# reduce first number to scale up image" << '\n'
- << "1200 2" << '\n';
-
- // We write all cells and cells on
- // coarser levels are behind cells
- // on finer levels. Level 0
- // corresponds to a depth of 900,
- // each level subtracting 1
- typename Triangulation<dim>::cell_iterator cell = tria.begin();
- const typename Triangulation<dim>::cell_iterator end = tria.end();
-
- for (;cell != end; ++cell)
+ const unsigned int nv = GeometryInfo<dim>::vertices_per_cell;
+ const unsigned int nf = GeometryInfo<dim>::faces_per_cell;
+ const unsigned int nvf = GeometryInfo<dim>::vertices_per_face;
+
+ // The following text was copied
+ // from an existing XFig file.
+ out << "#FIG 3.2\nLandscape\nCenter\nInches" << '\n'
+ << "A4\n100.00\nSingle" << '\n'
+ // Background is transparent
+ << "-3" << '\n'
+ << "# generated by deal.II GridOut class" << '\n'
+ << "# reduce first number to scale up image" << '\n'
+ << "1200 2" << '\n';
+
+ // We write all cells and cells on
+ // coarser levels are behind cells
+ // on finer levels. Level 0
+ // corresponds to a depth of 900,
+ // each level subtracting 1
+ typename Triangulation<dim>::cell_iterator cell = tria.begin();
+ const typename Triangulation<dim>::cell_iterator end = tria.end();
+
+ for (;cell != end; ++cell)
{
// If depth is not encoded, write finest level only
if (!xfig_flags.level_depth && !cell->active())
continue;
- // Code for polygon
- out << "2 3 "
- << xfig_flags.line_style << ' '
- << xfig_flags.line_thickness
- // with black line
- << " 0 ";
- // Fill color
- if (xfig_flags.level_color)
- out << cell->level() + 8;
- else
- out << cell->material_id() + 1;
- // Depth, unused, fill
- out << ' '
- << (xfig_flags.level_depth
- ? (900-cell->level())
- : (900+cell->material_id()))
- << " 0 "
- << xfig_flags.fill_style << " 0.0 "
- // some style parameters
- << " 0 0 -1 0 0 "
- // number of points
- << nv+1 << '\n';
-
- // For each point, write scaled
- // and shifted coordinates
- // multiplied by 1200
- // (dots/inch)
- for (unsigned int k=0;k<=nv;++k)
+ // Code for polygon
+ out << "2 3 "
+ << xfig_flags.line_style << ' '
+ << xfig_flags.line_thickness
+ // with black line
+ << " 0 ";
+ // Fill color
+ if (xfig_flags.level_color)
+ out << cell->level() + 8;
+ else
+ out << cell->material_id() + 1;
+ // Depth, unused, fill
+ out << ' '
+ << (xfig_flags.level_depth
+ ? (900-cell->level())
+ : (900+cell->material_id()))
+ << " 0 "
+ << xfig_flags.fill_style << " 0.0 "
+ // some style parameters
+ << " 0 0 -1 0 0 "
+ // number of points
+ << nv+1 << '\n';
+
+ // For each point, write scaled
+ // and shifted coordinates
+ // multiplied by 1200
+ // (dots/inch)
+ for (unsigned int k=0;k<=nv;++k)
{
- const Point<dim>& p = cell->vertex(
- GeometryInfo<dim>::ucd_to_deal[k % nv]);
- for (unsigned int d=0;d<dim;++d)
+ const Point<dim>& p = cell->vertex(
+ GeometryInfo<dim>::ucd_to_deal[k % nv]);
+ for (unsigned int d=0;d<dim;++d)
{
- int val = (int)(1200 * xfig_flags.scaling(d) *
- (p(d)-xfig_flags.offset(d)));
- out << '\t' << val;
+ int val = (int)(1200 * xfig_flags.scaling(d) *
+ (p(d)-xfig_flags.offset(d)));
+ out << '\t' << val;
}
- out << '\n';
+ out << '\n';
}
- // Now write boundary edges
- static const unsigned int face_reorder[4]={2,1,3,0};
- if (xfig_flags.draw_boundary)
- for (unsigned int f=0;f<nf;++f)
- {
- typename Triangulation<dim>::face_iterator
- face = cell->face(face_reorder[f]);
- const unsigned char bi = face->boundary_indicator();
- if (bi != 255)
- {
- // Code for polyline
- out << "2 1 "
- // with line style and thickness
- << xfig_flags.boundary_style << ' '
- << xfig_flags.boundary_thickness << ' '
- << (1 + (unsigned int) bi);
- // Fill color
- out << " -1 ";
- // Depth 100 less than cells
- out << (xfig_flags.level_depth
- ? (800-cell->level())
- : 800+bi)
- // unused, no fill
- << " 0 -1 0.0 "
- // some style parameters
- << " 0 0 -1 0 0 "
- // number of points
- << nvf << '\n';
-
- // For each point, write scaled
- // and shifted coordinates
- // multiplied by 1200
- // (dots/inch)
-
- for (unsigned int k=0;k<nvf;++k)
- {
- const Point<dim>& p = face->vertex(k % nv);
- for (unsigned int d=0;d<dim;++d)
- {
- int val = (int)(1200 * xfig_flags.scaling(d) *
- (p(d)-xfig_flags.offset(d)));
- out << '\t' << val;
- }
- out << '\n';
- }
- }
- }
+ // Now write boundary edges
+ static const unsigned int face_reorder[4]={2,1,3,0};
+ if (xfig_flags.draw_boundary)
+ for (unsigned int f=0;f<nf;++f)
+ {
+ typename Triangulation<dim>::face_iterator
+ face = cell->face(face_reorder[f]);
+ const unsigned char bi = face->boundary_indicator();
+ if (bi != 255)
+ {
+ // Code for polyline
+ out << "2 1 "
+ // with line style and thickness
+ << xfig_flags.boundary_style << ' '
+ << xfig_flags.boundary_thickness << ' '
+ << (1 + (unsigned int) bi);
+ // Fill color
+ out << " -1 ";
+ // Depth 100 less than cells
+ out << (xfig_flags.level_depth
+ ? (800-cell->level())
+ : 800+bi)
+ // unused, no fill
+ << " 0 -1 0.0 "
+ // some style parameters
+ << " 0 0 -1 0 0 "
+ // number of points
+ << nvf << '\n';
+
+ // For each point, write scaled
+ // and shifted coordinates
+ // multiplied by 1200
+ // (dots/inch)
+
+ for (unsigned int k=0;k<nvf;++k)
+ {
+ const Point<dim>& p = face->vertex(k % nv);
+ for (unsigned int d=0;d<dim;++d)
+ {
+ int val = (int)(1200 * xfig_flags.scaling(d) *
+ (p(d)-xfig_flags.offset(d)));
+ out << '\t' << val;
+ }
+ out << '\n';
+ }
+ }
+ }
}
- // make sure everything now gets to
- // disk
- out.flush ();
+ // make sure everything now gets to
+ // disk
+ out.flush ();
- AssertThrow (out, ExcIO());
+ AssertThrow (out, ExcIO());
}
#endif
unsigned int GridOut::n_boundary_faces (const Triangulation<1> &) const
{
- return 0;
+ return 0;
}
#endif
template <int dim>
unsigned int GridOut::n_boundary_faces (const Triangulation<dim> &tria) const
{
- typename Triangulation<dim>::active_face_iterator face, endf;
- unsigned int n_faces = 0;
+ typename Triangulation<dim>::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++;
+ 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;
+ return n_faces;
+}
+
+
+
+#if deal_II_dimension == 1
+
+void GridOut::write_msh_faces (const Triangulation<1> &,
+ const unsigned int,
+ std::ostream &) const
+{
+ return;
+}
+
+#endif
+
+
+
+template <int dim>
+void GridOut::write_msh_faces (const Triangulation<dim> &tria,
+ const unsigned int starting_index,
+ std::ostream &out) const
+{
+ typename Triangulation<dim>::active_face_iterator face, endf;
+ unsigned int index=starting_index;
+
+ for (face=tria.begin_active_face(), endf=tria.end_face();
+ face != endf; ++face)
+ if (face->at_boundary() &&
+ (face->boundary_indicator() != 0))
+ {
+ out << index << " ";
+ switch (dim)
+ {
+ case 2: out << 1 << ' '; break;
+ case 3: out << 3 << ' '; break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ out << static_cast<unsigned int>(face->boundary_indicator())
+ << " 1 " << GeometryInfo<dim>::vertices_per_face;
+ // note: vertex numbers are 1-base
+ for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
+ out << ' ' << face->vertex_index(GeometryInfo<dim-1>::ucd_to_deal[vertex])+1;
+ out << '\n';
+
+ ++index;
+ };
}
#if deal_II_dimension == 1
void GridOut::write_ucd_faces (const Triangulation<1> &,
- const unsigned int,
- std::ostream &) const
+ const unsigned int,
+ std::ostream &) const
{
- return;
+ return;
}
#endif
template <int dim>
void GridOut::write_ucd_faces (const Triangulation<dim> &tria,
- const unsigned int starting_index,
- std::ostream &out) const
+ const unsigned int starting_index,
+ std::ostream &out) const
{
- typename Triangulation<dim>::active_face_iterator face, endf;
- unsigned int index=starting_index;
+ typename Triangulation<dim>::active_face_iterator face, endf;
+ unsigned int index=starting_index;
- for (face=tria.begin_active_face(), endf=tria.end_face();
- face != endf; ++face)
- if (face->at_boundary() &&
- (face->boundary_indicator() != 0))
- {
- out << index << " "
- << static_cast<unsigned int>(face->boundary_indicator())
- << " ";
- switch (dim)
- {
- case 2: out << "line "; break;
- case 3: out << "quad "; break;
- default:
- Assert (false, ExcNotImplemented());
- };
- // note: vertex numbers are 1-base
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
- out << face->vertex_index(GeometryInfo<dim-1>::ucd_to_deal[vertex])+1 << ' ';
- out << '\n';
+ for (face=tria.begin_active_face(), endf=tria.end_face();
+ face != endf; ++face)
+ if (face->at_boundary() &&
+ (face->boundary_indicator() != 0))
+ {
+ out << index << " "
+ << static_cast<unsigned int>(face->boundary_indicator())
+ << " ";
+ switch (dim)
+ {
+ case 2: out << "line "; break;
+ case 3: out << "quad "; break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ // note: vertex numbers are 1-base
+ for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
+ out << face->vertex_index(GeometryInfo<dim-1>::ucd_to_deal[vertex])+1 << ' ';
+ out << '\n';
- ++index;
- };
+ ++index;
+ };
}
#if deal_II_dimension==1
void GridOut::write_gnuplot (const Triangulation<1> &tria,
- std::ostream &out,
- const Mapping<1> *)
+ std::ostream &out,
+ const Mapping<1> *)
{
- AssertThrow (out, ExcIO());
+ AssertThrow (out, ExcIO());
- Triangulation<1>::active_cell_iterator cell=tria.begin_active();
- const Triangulation<1>::active_cell_iterator endc=tria.end();
- for (; cell!=endc; ++cell)
+ Triangulation<1>::active_cell_iterator cell=tria.begin_active();
+ const Triangulation<1>::active_cell_iterator endc=tria.end();
+ for (; cell!=endc; ++cell)
{
- if (gnuplot_flags.write_cell_numbers)
- out << "# cell " << cell << '\n';
-
- out << cell->vertex(0)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(1)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << '\n';
- break;
+ if (gnuplot_flags.write_cell_numbers)
+ out << "# cell " << cell << '\n';
+
+ out << cell->vertex(0)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(1)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << '\n';
+ break;
}
- // make sure everything now gets to
- // disk
- out.flush ();
+ // make sure everything now gets to
+ // disk
+ out.flush ();
- AssertThrow (out, ExcIO());
+ AssertThrow (out, ExcIO());
}
#endif
#if deal_II_dimension==2
void GridOut::write_gnuplot (const Triangulation<2> &tria,
- std::ostream &out,
- const Mapping<2> *mapping)
+ std::ostream &out,
+ const Mapping<2> *mapping)
{
- AssertThrow (out, ExcIO());
-
- const unsigned int dim=2;
- const unsigned int n_additional_points=
- gnuplot_flags.n_boundary_face_points;
- const unsigned int n_points=2+n_additional_points;
-
- Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
- const Triangulation<dim>::active_cell_iterator endc=tria.end();
-
- // if we are to treat curved
- // boundaries, then generate a
- // quadrature formula which will be
- // used to probe boundary points at
- // curved faces
- Quadrature<dim> *q_projector=0;
- std::vector<Point<1> > boundary_points;
- if (mapping!=0)
+ AssertThrow (out, ExcIO());
+
+ const unsigned int dim=2;
+ const unsigned int n_additional_points=
+ gnuplot_flags.n_boundary_face_points;
+ const unsigned int n_points=2+n_additional_points;
+
+ Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
+ const Triangulation<dim>::active_cell_iterator endc=tria.end();
+
+ // if we are to treat curved
+ // boundaries, then generate a
+ // quadrature formula which will be
+ // used to probe boundary points at
+ // curved faces
+ Quadrature<dim> *q_projector=0;
+ std::vector<Point<1> > boundary_points;
+ if (mapping!=0)
{
- boundary_points.resize(n_points);
- boundary_points[0][0]=0;
- boundary_points[n_points-1][0]=1;
- for (unsigned int i=1; i<n_points-1; ++i)
- boundary_points[i](0)= 1.*i/(n_points-1);
+ boundary_points.resize(n_points);
+ boundary_points[0][0]=0;
+ boundary_points[n_points-1][0]=1;
+ for (unsigned int i=1; i<n_points-1; ++i)
+ boundary_points[i](0)= 1.*i/(n_points-1);
- std::vector<double> dummy_weights(n_points, 1./n_points);
- Quadrature<1> quadrature(boundary_points, dummy_weights);
+ std::vector<double> dummy_weights(n_points, 1./n_points);
+ Quadrature<1> quadrature(boundary_points, dummy_weights);
- q_projector = new Quadrature<dim> (QProjector<dim>::project_to_all_faces(quadrature));
+ q_projector = new Quadrature<dim> (QProjector<dim>::project_to_all_faces(quadrature));
}
-
- for (; cell!=endc; ++cell)
+
+ for (; cell!=endc; ++cell)
{
- if (gnuplot_flags.write_cell_numbers)
- out << "# cell " << cell << '\n';
+ if (gnuplot_flags.write_cell_numbers)
+ out << "# cell " << cell << '\n';
- if (mapping==0 || !cell->at_boundary())
+ if (mapping==0 || !cell->at_boundary())
{
- // write out the four sides
- // of this cell by putting
- // the four points (+ the
- // initial point again) in
- // a row and lifting the
- // drawing pencil at the
- // end
- for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- out << cell->vertex(GeometryInfo<dim>::ucd_to_deal[i])
+ // write out the four sides
+ // of this cell by putting
+ // the four points (+ the
+ // initial point again) in
+ // a row and lifting the
+ // drawing pencil at the
+ // end
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ out << cell->vertex(GeometryInfo<dim>::ucd_to_deal[i])
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+ out << cell->vertex(0)
<< ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
- out << cell->vertex(0)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << '\n' // double new line for gnuplot 3d plots
- << '\n';
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << '\n' // double new line for gnuplot 3d plots
+ << '\n';
}
- else
- // cell is at boundary and we
- // are to treat curved
- // boundaries. so loop over
- // all faces and draw them as
- // small pieces of lines
+ else
+ // cell is at boundary and we
+ // are to treat curved
+ // boundaries. so loop over
+ // all faces and draw them as
+ // small pieces of lines
{
- for (unsigned int face_no=0;
- face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ for (unsigned int face_no=0;
+ face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
- const Triangulation<dim>::face_iterator
- face = cell->face(face_no);
- if (face->at_boundary())
+ const Triangulation<dim>::face_iterator
+ face = cell->face(face_no);
+ if (face->at_boundary())
{
- // compute offset
- // of quadrature
- // points within
- // set of projected
- // points
- const unsigned int offset=face_no*n_points;
- for (unsigned int i=0; i<n_points; ++i)
- out << (mapping->transform_unit_to_real_cell
- (cell, q_projector->point(offset+i)))
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id())
+ // compute offset
+ // of quadrature
+ // points within
+ // set of projected
+ // points
+ const unsigned int offset=face_no*n_points;
+ for (unsigned int i=0; i<n_points; ++i)
+ out << (mapping->transform_unit_to_real_cell
+ (cell, q_projector->point(offset+i)))
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id())
+ << '\n';
+
+ out << '\n'
<< '\n';
-
- out << '\n'
- << '\n';
}
- else
+ else
{
- // if, however, the
- // face is not at
- // the boundary,
- // then draw it as
- // usual
- out << face->vertex(0)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id())
- << '\n'
- << face->vertex(1)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id())
- << '\n'
- << '\n'
- << '\n';
+ // if, however, the
+ // face is not at
+ // the boundary,
+ // then draw it as
+ // usual
+ out << face->vertex(0)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id())
+ << '\n'
+ << face->vertex(1)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id())
+ << '\n'
+ << '\n'
+ << '\n';
}
}
}
}
- if (q_projector != 0)
- delete q_projector;
-
- // make sure everything now gets to
- // disk
- out.flush ();
+ if (q_projector != 0)
+ delete q_projector;
- AssertThrow (out, ExcIO());
+ // make sure everything now gets to
+ // disk
+ out.flush ();
+
+ AssertThrow (out, ExcIO());
}
#endif
#if deal_II_dimension==3
void GridOut::write_gnuplot (const Triangulation<3> &tria,
- std::ostream &out,
- const Mapping<3> *mapping)
+ std::ostream &out,
+ const Mapping<3> *mapping)
{
- AssertThrow (out, ExcIO());
-
- const unsigned int dim=3;
- const unsigned int n_additional_points=
- gnuplot_flags.n_boundary_face_points;
- const unsigned int n_points=2+n_additional_points;
-
- Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
- const Triangulation<dim>::active_cell_iterator endc=tria.end();
-
- // if we are to treat curved
- // boundaries, then generate a
- // quadrature formula which will be
- // used to probe boundary points at
- // curved faces
- Quadrature<dim> *q_projector=0;
- std::vector<Point<1> > boundary_points;
- if (mapping!=0)
+ AssertThrow (out, ExcIO());
+
+ const unsigned int dim=3;
+ const unsigned int n_additional_points=
+ gnuplot_flags.n_boundary_face_points;
+ const unsigned int n_points=2+n_additional_points;
+
+ Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
+ const Triangulation<dim>::active_cell_iterator endc=tria.end();
+
+ // if we are to treat curved
+ // boundaries, then generate a
+ // quadrature formula which will be
+ // used to probe boundary points at
+ // curved faces
+ Quadrature<dim> *q_projector=0;
+ std::vector<Point<1> > boundary_points;
+ if (mapping!=0)
{
- boundary_points.resize(n_points);
- boundary_points[0][0]=0;
- boundary_points[n_points-1][0]=1;
- for (unsigned int i=1; i<n_points-1; ++i)
- boundary_points[i](0)= 1.*i/(n_points-1);
-
- std::vector<double> dummy_weights(n_points, 1./n_points);
- Quadrature<1> quadrature1d(boundary_points, dummy_weights);
-
- // tensor product of points,
- // only one copy
- QIterated<dim-1> quadrature(quadrature1d, 1);
- q_projector = new Quadrature<dim> (QProjector<dim>::project_to_all_faces(quadrature));
+ boundary_points.resize(n_points);
+ boundary_points[0][0]=0;
+ boundary_points[n_points-1][0]=1;
+ for (unsigned int i=1; i<n_points-1; ++i)
+ boundary_points[i](0)= 1.*i/(n_points-1);
+
+ std::vector<double> dummy_weights(n_points, 1./n_points);
+ Quadrature<1> quadrature1d(boundary_points, dummy_weights);
+
+ // tensor product of points,
+ // only one copy
+ QIterated<dim-1> quadrature(quadrature1d, 1);
+ q_projector = new Quadrature<dim> (QProjector<dim>::project_to_all_faces(quadrature));
}
-
- for (; cell!=endc; ++cell)
+
+ for (; cell!=endc; ++cell)
{
- if (gnuplot_flags.write_cell_numbers)
- out << "# cell " << cell << '\n';
+ if (gnuplot_flags.write_cell_numbers)
+ out << "# cell " << cell << '\n';
- if (mapping==0 || n_points==2 || !cell->has_boundary_lines())
+ if (mapping==0 || n_points==2 || !cell->has_boundary_lines())
{
- // front face
- out << cell->vertex(0)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(1)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(5)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(4)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(0)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << '\n';
- // back face
- out << cell->vertex(2)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(3)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(7)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(6)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(2)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << '\n';
-
- // now for the four connecting lines
- out << cell->vertex(0)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(2)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << '\n';
- out << cell->vertex(1)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(3)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << '\n';
- out << cell->vertex(5)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(7)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << '\n';
- out << cell->vertex(4)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << cell->vertex(6)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << '\n';
+ // front face
+ out << cell->vertex(0)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(1)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(5)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(4)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(0)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << '\n';
+ // back face
+ out << cell->vertex(2)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(3)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(7)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(6)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(2)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << '\n';
+
+ // now for the four connecting lines
+ out << cell->vertex(0)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(2)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << '\n';
+ out << cell->vertex(1)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(3)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << '\n';
+ out << cell->vertex(5)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(7)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << '\n';
+ out << cell->vertex(4)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << cell->vertex(6)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << '\n';
}
- else
+ else
{
- for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
- const Triangulation<dim>::face_iterator
- face = cell->face(face_no);
-
- if (face->at_boundary())
+ const Triangulation<dim>::face_iterator
+ face = cell->face(face_no);
+
+ if (face->at_boundary())
{
- const unsigned int offset=face_no*n_points*n_points;
- for (unsigned int i=0; i<n_points-1; ++i)
- for (unsigned int j=0; j<n_points-1; ++j)
- {
- const Point<dim> p0=mapping->transform_unit_to_real_cell(
- cell, q_projector->point(offset+i*n_points+j));
- out << p0
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
- out << (mapping->transform_unit_to_real_cell(
- cell, q_projector->point(offset+(i+1)*n_points+j)))
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
- out << (mapping->transform_unit_to_real_cell(
- cell, q_projector->point(offset+(i+1)*n_points+j+1)))
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
- out << (mapping->transform_unit_to_real_cell(
- cell, q_projector->point(offset+i*n_points+j+1)))
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
- // and the
- // first
- // point
- // again
- out << p0
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
- out << '\n' << '\n';
- }
+ const unsigned int offset=face_no*n_points*n_points;
+ for (unsigned int i=0; i<n_points-1; ++i)
+ for (unsigned int j=0; j<n_points-1; ++j)
+ {
+ const Point<dim> p0=mapping->transform_unit_to_real_cell(
+ cell, q_projector->point(offset+i*n_points+j));
+ out << p0
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+ out << (mapping->transform_unit_to_real_cell(
+ cell, q_projector->point(offset+(i+1)*n_points+j)))
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+ out << (mapping->transform_unit_to_real_cell(
+ cell, q_projector->point(offset+(i+1)*n_points+j+1)))
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+ out << (mapping->transform_unit_to_real_cell(
+ cell, q_projector->point(offset+i*n_points+j+1)))
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+ // and the
+ // first
+ // point
+ // again
+ out << p0
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+ out << '\n' << '\n';
+ }
}
- else
+ else
{
- for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
+ for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
{
- const Triangulation<dim>::line_iterator
- line=face->line(l);
-
- const Point<dim> &v0=line->vertex(0),
- &v1=line->vertex(1);
- if (line->at_boundary())
+ const Triangulation<dim>::line_iterator
+ line=face->line(l);
+
+ const Point<dim> &v0=line->vertex(0),
+ &v1=line->vertex(1);
+ if (line->at_boundary())
{
- // transform_real_to_unit_cell
- // could be
- // replaced
- // by using
- // QProjector<3>::project_to_line
- // which is
- // not yet
- // implemented
- const Point<dim> u0=mapping->transform_real_to_unit_cell(cell, v0),
- u1=mapping->transform_real_to_unit_cell(cell, v1);
-
- for (unsigned int i=0; i<n_points; ++i)
- out << (mapping->transform_unit_to_real_cell
- (cell, (1-boundary_points[i][0])*u0+boundary_points[i][0]*u1))
+ // transform_real_to_unit_cell
+ // could be
+ // replaced
+ // by using
+ // QProjector<3>::project_to_line
+ // which is
+ // not yet
+ // implemented
+ const Point<dim> u0=mapping->transform_real_to_unit_cell(cell, v0),
+ u1=mapping->transform_real_to_unit_cell(cell, v1);
+
+ for (unsigned int i=0; i<n_points; ++i)
+ out << (mapping->transform_unit_to_real_cell
+ (cell, (1-boundary_points[i][0])*u0+boundary_points[i][0]*u1))
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
+ }
+ else
+ out << v0
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
+ << v1
<< ' ' << cell->level()
<< ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
- }
- else
- out << v0
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n'
- << v1
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
- out << '\n' << '\n';
+ out << '\n' << '\n';
}
}
}
}
}
- if (q_projector != 0)
- delete q_projector;
-
-
- // make sure everything now gets to
- // disk
- out.flush ();
+ if (q_projector != 0)
+ delete q_projector;
- AssertThrow (out, ExcIO());
+
+ // make sure everything now gets to
+ // disk
+ out.flush ();
+
+ AssertThrow (out, ExcIO());
}
#endif
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)
- {}
+ const Point<2> &s,
+ const bool c,
+ const unsigned int l)
+ :
+ first(f), second(s),
+ colorize(c), level(l)
+ {}
};
#if deal_II_dimension==1
void GridOut::write_eps (const Triangulation<1> &,
- std::ostream &,
- const Mapping<1> *)
+ std::ostream &,
+ const Mapping<1> *)
{
- Assert(false, ExcNotImplemented());
+ Assert(false, ExcNotImplemented());
}
#else
-template <int dim>
+ template <int dim>
void GridOut::write_eps (const Triangulation<dim> &tria,
- std::ostream &out,
- const Mapping<dim> *mapping)
+ std::ostream &out,
+ const Mapping<dim> *mapping)
{
-
- typedef std::list<LineEntry> LineList;
-
- // get a pointer to the flags
- // common to all dimensions,
- // in order to avoid the recurring
- // distinctions between
- // eps_flags_1, eps_flags_2, ...
- const GridOutFlags::EpsFlagsBase
- &eps_flags_base = (dim==2 ?
- static_cast<GridOutFlags::EpsFlagsBase&>(eps_flags_2) :
- (dim==3 ?
- static_cast<GridOutFlags::EpsFlagsBase&>(eps_flags_3) :
- *static_cast<GridOutFlags::EpsFlagsBase*>(0)));
-
- AssertThrow (out, ExcIO());
- const unsigned int n_points = eps_flags_base.n_boundary_face_points;
-
- // make up a list of lines by which
- // we will construct the triangulation
- //
- // this part unfortunately is a bit
- // dimension dependent, so we have to
- // treat every dimension different.
- // however, by directly producing
- // the lines to be printed, i.e. their
- // 2d images, we can later do the
- // actual output dimension independent
- // again
- LineList line_list;
-
- switch (dim)
+
+ typedef std::list<LineEntry> LineList;
+
+ // get a pointer to the flags
+ // common to all dimensions,
+ // in order to avoid the recurring
+ // distinctions between
+ // eps_flags_1, eps_flags_2, ...
+ const GridOutFlags::EpsFlagsBase
+ &eps_flags_base = (dim==2 ?
+ static_cast<GridOutFlags::EpsFlagsBase&>(eps_flags_2) :
+ (dim==3 ?
+ static_cast<GridOutFlags::EpsFlagsBase&>(eps_flags_3) :
+ *static_cast<GridOutFlags::EpsFlagsBase*>(0)));
+
+ AssertThrow (out, ExcIO());
+ const unsigned int n_points = eps_flags_base.n_boundary_face_points;
+
+ // make up a list of lines by which
+ // we will construct the triangulation
+ //
+ // this part unfortunately is a bit
+ // dimension dependent, so we have to
+ // treat every dimension different.
+ // however, by directly producing
+ // the lines to be printed, i.e. their
+ // 2d images, we can later do the
+ // actual output dimension independent
+ // again
+ LineList line_list;
+
+ switch (dim)
{
- case 1:
- {
- Assert(false, ExcInternalError());
- break;
- };
-
- case 2:
- {
- typename Triangulation<dim>::active_line_iterator
- line =tria.begin_active_line (),
- endline=tria.end_line ();
-
- // 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
- for (; line!=endline; ++line)
- if (mapping==0 || !line->at_boundary())
- // one would expect
- // make_pair(line->vertex(0),
- // line->vertex(1))
- // here, but that is not
- // dimension independent, since
- // vertex(i) is Point<dim>,
- // but we want a Point<2>.
- // in fact, whenever we're here,
- // the vertex is a Point<dim>,
- // but the compiler does not
- // know this. hopefully, the
- // compiler will optimize away
- // this little kludge
- line_list.push_back (LineEntry(Point<2>(line->vertex(0)(0),
- line->vertex(0)(1)),
- Point<2>(line->vertex(1)(0),
- line->vertex(1)(1)),
- line->user_flag_set(),
- line->level()));
-
- // next if we are to treat
- // curved boundaries
- // specially, then add lines
- // to the list consisting of
- // pieces of the boundary
- // lines
- if (mapping!=0)
- {
- // to do so, first
- // generate a sequence of
- // points on a face and
- // project them onto the
- // faces of a unit cell
- std::vector<Point<dim-1> > boundary_points (n_points);
-
- for (unsigned int i=0; i<n_points; ++i)
- boundary_points[i](0) = 1.*(i+1)/(n_points+1);
-
- Quadrature<dim-1> quadrature (boundary_points);
- Quadrature<dim> q_projector (QProjector<dim>::project_to_all_faces(quadrature));
-
- // next loop over all
- // boundary faces and
- // generate the info from
- // them
- typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active ();
- const typename Triangulation<dim>::active_cell_iterator end=tria.end ();
- for (; cell!=end; ++cell)
- for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ case 1:
+ {
+ Assert(false, ExcInternalError());
+ break;
+ };
+
+ case 2:
+ {
+ typename Triangulation<dim>::active_line_iterator
+ line =tria.begin_active_line (),
+ endline=tria.end_line ();
+
+ // 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
+ for (; line!=endline; ++line)
+ if (mapping==0 || !line->at_boundary())
+ // one would expect
+ // make_pair(line->vertex(0),
+ // line->vertex(1))
+ // here, but that is not
+ // dimension independent, since
+ // vertex(i) is Point<dim>,
+ // but we want a Point<2>.
+ // in fact, whenever we're here,
+ // the vertex is a Point<dim>,
+ // but the compiler does not
+ // know this. hopefully, the
+ // compiler will optimize away
+ // this little kludge
+ line_list.push_back (LineEntry(Point<2>(line->vertex(0)(0),
+ line->vertex(0)(1)),
+ Point<2>(line->vertex(1)(0),
+ line->vertex(1)(1)),
+ line->user_flag_set(),
+ line->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)
{
- const typename Triangulation<dim>::face_iterator
- face = cell->face(face_no);
-
- if (face->at_boundary())
- {
- Point<dim> p0_dim(face->vertex(0));
- Point<2> p0 (p0_dim(0), p0_dim(1));
-
- // loop over
- // all pieces
- // of the line
- // and generate
- // line-lets
- const unsigned int offset=face_no*n_points;
- for (unsigned int i=0; i<n_points; ++i)
+ // to do so, first
+ // generate a sequence of
+ // points on a face and
+ // project them onto the
+ // faces of a unit cell
+ std::vector<Point<dim-1> > boundary_points (n_points);
+
+ for (unsigned int i=0; i<n_points; ++i)
+ boundary_points[i](0) = 1.*(i+1)/(n_points+1);
+
+ Quadrature<dim-1> quadrature (boundary_points);
+ Quadrature<dim> q_projector (QProjector<dim>::project_to_all_faces(quadrature));
+
+ // next loop over all
+ // boundary faces and
+ // generate the info from
+ // them
+ typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active ();
+ const typename Triangulation<dim>::active_cell_iterator end=tria.end ();
+ for (; cell!=end; ++cell)
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
- const Point<dim> p1_dim (mapping->transform_unit_to_real_cell
- (cell, q_projector.point(offset+i)));
- const Point<2> p1 (p1_dim(0), p1_dim(1));
-
- line_list.push_back (LineEntry(p0, p1,
- face->user_flag_set(),
- face->level() ));
- p0=p1;
- }
-
- // generate last piece
- const Point<dim> p1_dim (face->vertex(1));
- const Point<2> p1 (p1_dim(0), p1_dim(1));
- line_list.push_back (LineEntry(p0, p1,
- face->user_flag_set(),
- face->level()));
- };
+ const typename Triangulation<dim>::face_iterator
+ face = cell->face(face_no);
+
+ if (face->at_boundary())
+ {
+ Point<dim> p0_dim(face->vertex(0));
+ Point<2> p0 (p0_dim(0), p0_dim(1));
+
+ // loop over
+ // all pieces
+ // of the line
+ // and generate
+ // line-lets
+ const unsigned int offset=face_no*n_points;
+ for (unsigned int i=0; i<n_points; ++i)
+ {
+ const Point<dim> p1_dim (mapping->transform_unit_to_real_cell
+ (cell, q_projector.point(offset+i)));
+ const Point<2> p1 (p1_dim(0), p1_dim(1));
+
+ line_list.push_back (LineEntry(p0, p1,
+ face->user_flag_set(),
+ face->level() ));
+ p0=p1;
+ }
+
+ // generate last piece
+ const Point<dim> p1_dim (face->vertex(1));
+ const Point<2> p1 (p1_dim(0), p1_dim(1));
+ line_list.push_back (LineEntry(p0, p1,
+ face->user_flag_set(),
+ face->level()));
+ };
+ };
};
- };
-
- break;
- };
-
- case 3:
- {
- // curved boundary output
- // presently not supported
-//TODO:[RH] curved boundaries in eps for 3d
- Assert (mapping == 0, ExcNotImplemented());
-
- typename Triangulation<dim>::active_line_iterator
- line =tria.begin_active_line (),
- endline=tria.end_line ();
-
- // 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 = deal_II_numbers::PI;
- const double z_angle = eps_flags_3.azimut_angle;
- const double turn_angle = eps_flags_3.turn_angle;
- const Point<dim> view_direction(-std::sin(z_angle * 2.*pi / 360.) * std::sin(turn_angle * 2.*pi / 360.),
- +std::sin(z_angle * 2.*pi / 360.) * std::cos(turn_angle * 2.*pi / 360.),
- -std::cos(z_angle * 2.*pi / 360.));
-
- // decide about the two unit vectors
- // in this plane. we chose the first one
- // to be the projection of the z-axis
- // to this plane
- const Point<dim> vector1
- = Point<dim>(0,0,1) - ((Point<dim>(0,0,1) * view_direction) * view_direction);
- const Point<dim> unit_vector1 = vector1 / std::sqrt(vector1.square());
-
- // now the third vector is fixed. we
- // chose the projection of a more or
- // less arbitrary vector to the plane
- // perpendicular to the first one
- const Point<dim> vector2
- = (Point<dim>(1,0,0)
- - ((Point<dim>(1,0,0) * view_direction) * view_direction)
- - ((Point<dim>(1,0,0) * unit_vector1) * unit_vector1));
- const Point<dim> unit_vector2 = vector2 / std::sqrt(vector2.square());
-
- for (; line!=endline; ++line)
- 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(),
- line->level()));
- break;
- };
+ break;
+ };
- default:
+ case 3:
+ {
+ // curved boundary output
+ // presently not supported
+ //TODO:[RH] curved boundaries in eps for 3d
+ Assert (mapping == 0, ExcNotImplemented());
+
+ typename Triangulation<dim>::active_line_iterator
+ line =tria.begin_active_line (),
+ endline=tria.end_line ();
+
+ // 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 = deal_II_numbers::PI;
+ const double z_angle = eps_flags_3.azimut_angle;
+ const double turn_angle = eps_flags_3.turn_angle;
+ const Point<dim> view_direction(-std::sin(z_angle * 2.*pi / 360.) * std::sin(turn_angle * 2.*pi / 360.),
+ +std::sin(z_angle * 2.*pi / 360.) * std::cos(turn_angle * 2.*pi / 360.),
+ -std::cos(z_angle * 2.*pi / 360.));
+
+ // decide about the two unit vectors
+ // in this plane. we chose the first one
+ // to be the projection of the z-axis
+ // to this plane
+ const Point<dim> vector1
+ = Point<dim>(0,0,1) - ((Point<dim>(0,0,1) * view_direction) * view_direction);
+ const Point<dim> unit_vector1 = vector1 / std::sqrt(vector1.square());
+
+ // now the third vector is fixed. we
+ // chose the projection of a more or
+ // less arbitrary vector to the plane
+ // perpendicular to the first one
+ const Point<dim> vector2
+ = (Point<dim>(1,0,0)
+ - ((Point<dim>(1,0,0) * view_direction) * view_direction)
+ - ((Point<dim>(1,0,0) * unit_vector1) * unit_vector1));
+ const Point<dim> unit_vector2 = vector2 / std::sqrt(vector2.square());
+
+ for (; line!=endline; ++line)
+ 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(),
+ line->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;
+ // 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)
+ 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_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));
+ 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_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);
+ 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));
+ // 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)
+ // now write preamble
+ if (true)
{
- // block this to have local
- // variables destroyed after
- // use
- std::time_t time1= std::time (0);
- std::tm *time = std::localtime(&time1);
- out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
- << "%%Title: deal.II Output" << '\n'
- << "%%Creator: the deal.II library" << '\n'
- << "%%Creation Date: "
- << time->tm_year+1900 << "/"
- << time->tm_mon+1 << "/"
- << time->tm_mday << " - "
- << time->tm_hour << ":"
- << std::setw(2) << time->tm_min << ":"
- << std::setw(2) << time->tm_sec << '\n'
- << "%%BoundingBox: "
- // lower left corner
- << "0 0 "
- // upper right corner
- << static_cast<unsigned int>(std::floor(( (x_max-x_min) * scale )+1))
- << ' '
- << static_cast<unsigned int>(std::floor(( (y_max-y_min) * scale )+1))
- << '\n';
-
- // define some abbreviations to keep
- // the output small:
- // m=move turtle to
- // x=execute line stroke
- // b=black pen
- // r=red pen
- out << "/m {moveto} bind def" << '\n'
- << "/x {lineto stroke} bind def" << '\n'
- << "/b {0 0 0 setrgbcolor} def" << '\n'
- << "/r {1 0 0 setrgbcolor} def" << '\n';
-
- // calculate colors for level
- // coloring; level 0 is black,
- // other levels are blue
- // ... red
- if (eps_flags_base.color_lines_level)
- out << "/l { neg "
- << (max_level)
- << " add "
- << (0.66666/std::max(1U,(max_level-1)))
- << " mul 1 0.8 sethsbcolor} def" << '\n';
-
- // in 2d, we can also plot cell
- // and vertex numbers, but this
- // requires a somewhat more
- // lengthy preamble. please
- // don't ask me what most of
- // this means, it is reverse
- // engineered from what GNUPLOT
- // uses in its output
- if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
- eps_flags_2.write_vertex_numbers))
+ // block this to have local
+ // variables destroyed after
+ // use
+ std::time_t time1= std::time (0);
+ std::tm *time = std::localtime(&time1);
+ out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
+ << "%%Title: deal.II Output" << '\n'
+ << "%%Creator: the deal.II library" << '\n'
+ << "%%Creation Date: "
+ << time->tm_year+1900 << "/"
+ << time->tm_mon+1 << "/"
+ << time->tm_mday << " - "
+ << time->tm_hour << ":"
+ << std::setw(2) << time->tm_min << ":"
+ << std::setw(2) << time->tm_sec << '\n'
+ << "%%BoundingBox: "
+ // lower left corner
+ << "0 0 "
+ // upper right corner
+ << static_cast<unsigned int>(std::floor(( (x_max-x_min) * scale )+1))
+ << ' '
+ << static_cast<unsigned int>(std::floor(( (y_max-y_min) * scale )+1))
+ << '\n';
+
+ // define some abbreviations to keep
+ // the output small:
+ // m=move turtle to
+ // x=execute line stroke
+ // b=black pen
+ // r=red pen
+ out << "/m {moveto} bind def" << '\n'
+ << "/x {lineto stroke} bind def" << '\n'
+ << "/b {0 0 0 setrgbcolor} def" << '\n'
+ << "/r {1 0 0 setrgbcolor} def" << '\n';
+
+ // calculate colors for level
+ // coloring; level 0 is black,
+ // other levels are blue
+ // ... red
+ if (eps_flags_base.color_lines_level)
+ out << "/l { neg "
+ << (max_level)
+ << " add "
+ << (0.66666/std::max(1U,(max_level-1)))
+ << " mul 1 0.8 sethsbcolor} def" << '\n';
+
+ // in 2d, we can also plot cell
+ // and vertex numbers, but this
+ // requires a somewhat more
+ // lengthy preamble. please
+ // don't ask me what most of
+ // this means, it is reverse
+ // engineered from what GNUPLOT
+ // uses in its output
+ if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
+ eps_flags_2.write_vertex_numbers))
{
- out << ("/R {rmoveto} bind def\n"
- "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
- "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
- "currentdict end definefont\n"
- "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
- "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
- "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
- "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
- "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
- "5 get stringwidth pop add}\n"
- "{pop} ifelse} forall} bind def\n"
- "/MCshow { currentpoint stroke m\n"
- "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
- << '\n';
+ out << ("/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';
+ 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))
+ // now write the lines
+ const Point<2> offset(x_min, y_min);
+
+ for (LineList::const_iterator line=line_list.begin();
+ line!=line_list.end(); ++line)
+ if (eps_flags_base.color_lines_level && (line->level > 0))
+ // lines colored according to
+ // refinement level,
+ // contributed by Jörg
+ // R. Weimar
+ out << line->level
+ << " l "
+ << (line->first - offset) * scale << " m "
+ << (line->second - offset) * scale << " x" << '\n';
+ else
+ out << ((line->colorize && eps_flags_base.color_lines_on_user_flag) ? "r " : "b ")
+ << (line->first - offset) * scale << " m "
+ << (line->second - offset) * scale << " x" << '\n';
+
+ // finally write the cell numbers
+ // in 2d, if that is desired
+ if ((dim == 2) && (eps_flags_2.write_cell_numbers == true))
{
- out << "(Helvetica) findfont 140 scalefont setfont"
- << '\n';
-
- typename Triangulation<dim>::active_cell_iterator
- cell = tria.begin_active (),
- endc = tria.end ();
- for (; cell!=endc; ++cell)
+ out << "(Helvetica) findfont 140 scalefont setfont"
+ << '\n';
+
+ typename Triangulation<dim>::active_cell_iterator
+ cell = tria.begin_active (),
+ endc = tria.end ();
+ for (; cell!=endc; ++cell)
{
- out << (cell->center()(0)-offset(0))*scale << ' '
- << (cell->center()(1)-offset(1))*scale
- << " m" << '\n'
- << "[ [(Helvetica) 12.0 0.0 true true (";
- if (eps_flags_2.write_cell_number_level)
- out << cell;
- else
- out << cell->index();
-
- out << ")] "
- << "] -6 MCshow"
- << '\n';
+ 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))
+ // and the vertex numbers
+ if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
{
- out << "(Helvetica) findfont 140 scalefont setfont"
- << '\n';
-
- // have a list of those
- // vertices which we have
- // already tracked, to avoid
- // doing this multiply
- std::set<unsigned int> treated_vertices;
- typename Triangulation<dim>::active_cell_iterator
- cell = tria.begin_active (),
- endc = tria.end ();
- for (; cell!=endc; ++cell)
- for (unsigned int vertex=0;
- vertex<GeometryInfo<dim>::vertices_per_cell;
- ++vertex)
- if (treated_vertices.find(cell->vertex_index(vertex))
- ==
- treated_vertices.end())
- {
- treated_vertices.insert (cell->vertex_index(vertex));
-
- out << (cell->vertex(vertex)(0)-offset(0))*scale << ' '
- << (cell->vertex(vertex)(1)-offset(1))*scale
- << " m" << '\n'
- << "[ [(Helvetica) 10.0 0.0 true true ("
- << cell->vertex_index(vertex)
- << ")] "
- << "] -6 MCshow"
- << '\n';
- };
+ out << "(Helvetica) findfont 140 scalefont setfont"
+ << '\n';
+
+ // have a list of those
+ // vertices which we have
+ // already tracked, to avoid
+ // doing this multiply
+ std::set<unsigned int> treated_vertices;
+ typename Triangulation<dim>::active_cell_iterator
+ cell = tria.begin_active (),
+ endc = tria.end ();
+ for (; cell!=endc; ++cell)
+ for (unsigned int vertex=0;
+ vertex<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex)
+ if (treated_vertices.find(cell->vertex_index(vertex))
+ ==
+ treated_vertices.end())
+ {
+ treated_vertices.insert (cell->vertex_index(vertex));
+
+ out << (cell->vertex(vertex)(0)-offset(0))*scale << ' '
+ << (cell->vertex(vertex)(1)-offset(1))*scale
+ << " m" << '\n'
+ << "[ [(Helvetica) 10.0 0.0 true true ("
+ << cell->vertex_index(vertex)
+ << ")] "
+ << "] -6 MCshow"
+ << '\n';
+ };
};
-
- out << "showpage" << '\n';
-
- // make sure everything now gets to
- // disk
- out.flush ();
- AssertThrow (out, ExcIO());
+ out << "showpage" << '\n';
+
+ // make sure everything now gets to
+ // disk
+ out.flush ();
+
+ AssertThrow (out, ExcIO());
}
#endif
-template <int dim>
+ template <int dim>
void GridOut::write (const Triangulation<dim> &tria,
- std::ostream &out,
- const OutputFormat output_format,
- const Mapping<dim> *mapping)
+ std::ostream &out,
+ const OutputFormat output_format,
+ const Mapping<dim> *mapping)
{
- switch (output_format)
+ 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 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());
+
+ Assert (false, ExcInternalError());
}