From: Luca Heltai Date: Fri, 17 Feb 2006 18:09:51 +0000 (+0000) Subject: Added GridOut::write_gmsh. X-Git-Tag: v8.0.0~12283 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=59f752f255709e512345a17ac3b1e303125c2522;p=dealii.git Added GridOut::write_gmsh. git-svn-id: https://svn.dealii.org/trunk@12399 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/grid_in.h b/deal.II/deal.II/include/grid/grid_in.h index b2cff3da8e..333ba455a3 100644 --- a/deal.II/deal.II/include/grid/grid_in.h +++ b/deal.II/deal.II/include/grid/grid_in.h @@ -30,7 +30,7 @@ class SubCellData; /** * This class implements an input mechanism for grid data. It allows * to read a grid structure into a triangulation object. At present, - * only UCD (unstructured cell data) and DB Mesh is supported as input + * UCD (unstructured cell data), DB Mesh and Gmsh are supported as input * format for grid data. Any numerical data after the block of * topological information is ignored. * diff --git a/deal.II/deal.II/include/grid/grid_out.h b/deal.II/deal.II/include/grid/grid_out.h index 0d515a4429..20c33a4a3a 100644 --- a/deal.II/deal.II/include/grid/grid_out.h +++ b/deal.II/deal.II/include/grid/grid_out.h @@ -77,7 +77,41 @@ namespace GridOutFlags const bool write_measure = false, const bool write_all_faces = true); }; - + + /** + * Flags describing the details of + * output in MSH format. + * + * @ingroup output + */ + struct Msh + { + /** + * When writing a mesh, write boundary + * faces explicitly if their boundary + * indicator is not the default + * boundary indicator, which is zero. + * This is necessary if you later + * want to re-read the grid and want + * to get the same boundary indicators + * for the different parts of the + * boundary of the triangulation. + * + * It is not necessary if you only want + * to write the triangulation to + * view or print it. + * + * Default: @p false. + */ + bool write_faces; + + /** + * Constructor. + */ + Msh (const bool write_faces = false); + }; + + /** * Flags describing the details of * output in UCD format. @@ -582,7 +616,7 @@ namespace GridOutFlags * * @ingroup grid * @ingroup output - * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2003; postscript format based on an implementation by Stefan Nauber, 1999 + * @author Wolfgang Bangerth, Guido Kanschat, Luca Heltai, 1999, 2003, 2006; postscript format based on an implementation by Stefan Nauber, 1999 */ class GridOut { @@ -608,7 +642,9 @@ class GridOut /// write() calls write_ucd() ucd, /// write() calls write_xfig() - xfig + xfig, + /// write() calls write_msh() + msh }; /** @@ -703,6 +739,44 @@ class GridOut std::ostream &out, const Mapping<1> *mapping=0); + /** + * Write the triangulation in the + * msh format. + * + * Msh + * is the format used by Gmsh and it is + * described in the gmsh + * user's guide. Besides the + * usual output of the grid + * only, you can decide through + * additional flags (see below, + * and the documentation of the + * GridOutFlags::Msh() class) + * whether boundary faces with + * non-zero boundary indicator + * shall be written to the file + * explicitly. This is useful, + * if you want to re-read the + * grid later on, since + * deal.II sets the boundary + * indicator to zero by default; + * therefore, to obtain the + * same triangulation as before, + * you have to specify faces + * with differing boundary + * indicators explicitly, which + * is done by this flag. + * + * Names and values of further + * flags controlling the output + * can be found in the + * documentation of the + * GridOut::Msh() class. + */ + template + void write_msh (const Triangulation &tria, + std::ostream &out); + /** * Write the triangulation in the * ucd format. @@ -858,6 +932,11 @@ class GridOut */ void set_flags (const GridOutFlags::DX &flags); + /** + * Set flags for GMSH output + */ + void set_flags (const GridOutFlags::Msh &flags); + /** * Set flags for UCD output */ @@ -960,6 +1039,13 @@ class GridOut */ GridOutFlags::DX dx_flags; + /** + * Flags for GMSH output. Can be + * changed by using the + * @p set_flags function. + */ + GridOutFlags::Msh msh_flags; + /** * Flags for UCD output. Can be * changed by using the @@ -1002,7 +1088,45 @@ class GridOut * Flags used for XFig output. */ GridOutFlags::XFig xfig_flags; - + + /** + * Write the grid information about + * faces to @p out. Only those faces + * are printed which are on the boundary + * and which have a boundary indicator + * not equal to zero, since the latter + * is the default for boundary faces. + * + * Since cells and faces are continuously + * numbered, the @p starting_index for + * the numbering of the faces is passed + * also. + * + * This function unfortunately can not + * be included in the regular @p write_msh + * function, since it needs special + * treatment for the case dim==1, in + * which case the face iterators are + * void*'s and lack the member functions + * which are called. We would not actually + * call these functions, but the compiler + * would complain anyway when compiling + * the function for dim==1. Bad luck. + */ + template + void write_msh_faces (const Triangulation &tria, + const unsigned int starting_index, + std::ostream &out) const; + + /** + * Declaration of the specialization + * of above function for 1d. Does + * nothing. + */ + void write_msh_faces (const Triangulation<1> &tria, + const unsigned int starting_index, + std::ostream &out) const; + /** * Write the grid information about * faces to @p out. Only those faces diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc index 0b357817d5..f2e46282be 100644 --- a/deal.II/deal.II/source/grid/grid_generator.cc +++ b/deal.II/deal.II/source/grid/grid_generator.cc @@ -1394,7 +1394,7 @@ GridGenerator::half_hyper_ball (Triangulation& tria, { for (unsigned int i=0;i::faces_per_cell;++i) { - if (cell->face(i)->boundary_indicator() == 255) + if (!cell->at_boundary(i)) continue; // If the center is on the plane x=0, this is a planar @@ -1408,7 +1408,7 @@ GridGenerator::half_hyper_ball (Triangulation& tria, // With this loop we restore back the indicator of the outer lines for (unsigned int i=0;i::faces_per_cell;++i) { - if (cell->face(i)->boundary_indicator() == 255) + if (!cell->at_boundary(i)) continue; // If the center is not on the plane x=0, this is a curvilinear diff --git a/deal.II/deal.II/source/grid/grid_out.all_dimensions.cc b/deal.II/deal.II/source/grid/grid_out.all_dimensions.cc index 4cf57c2731..b2c1d057d4 100644 --- a/deal.II/deal.II/source/grid/grid_out.all_dimensions.cc +++ b/deal.II/deal.II/source/grid/grid_out.all_dimensions.cc @@ -31,6 +31,11 @@ namespace GridOutFlags + Msh::Msh (const bool write_faces) : + write_faces (write_faces) + {} + + Ucd::Ucd (const bool write_preamble, const bool write_faces) : write_preamble (write_preamble), @@ -138,6 +143,12 @@ void GridOut::set_flags (const GridOutFlags::DX &flags) +void GridOut::set_flags (const GridOutFlags::Msh &flags) +{ + msh_flags = flags; +} + + void GridOut::set_flags (const GridOutFlags::Ucd &flags) { ucd_flags = flags; @@ -197,6 +208,8 @@ GridOut::default_suffix (const OutputFormat output_format) return ".eps"; case xfig: return ".fig"; + case msh: + return ".msh"; default: Assert (false, ExcNotImplemented()); return ""; @@ -226,6 +239,9 @@ GridOut::parse_output_format (const std::string &format_name) if (format_name == "xfig") return xfig; + if (format_name == "msh") + return msh; + AssertThrow (false, ExcInvalidState ()); // return something weird return OutputFormat(-1); @@ -235,7 +251,7 @@ GridOut::parse_output_format (const std::string &format_name) std::string GridOut::get_output_format_names () { - return "none|dx|gnuplot|eps|ucd|xfig"; + return "none|dx|gnuplot|eps|ucd|xfig|msh"; } diff --git a/deal.II/deal.II/source/grid/grid_out.cc b/deal.II/deal.II/source/grid/grid_out.cc index e87ec6df0e..7d4952990e 100644 --- a/deal.II/deal.II/source/grid/grid_out.cc +++ b/deal.II/deal.II/source/grid/grid_out.cc @@ -265,7 +265,7 @@ void GridOut::write_dx (const Triangulation &tria, template -void GridOut::write_ucd (const Triangulation &tria, +void GridOut::write_msh (const Triangulation &tria, std::ostream &out) { AssertThrow (out, ExcIO()); @@ -281,242 +281,380 @@ void GridOut::write_ucd (const Triangulation &tria, 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) - << " 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(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 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 +void GridOut::write_ucd (const Triangulation &tria, + std::ostream &out) +{ + 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) { - out << cell_index << ' ' - << static_cast(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_per_cell; - ++vertex) - out << cell->vertex_index(GeometryInfo::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(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 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 + template void GridOut::write_xfig (const Triangulation&, - std::ostream&, - const Mapping*) + std::ostream&, + const Mapping*) { - Assert (false, ExcNotImplemented()); + Assert (false, ExcNotImplemented()); } #else //TODO:[GK] Obey parameters //TODO:[GK] Flip y-axis? -template + template void GridOut::write_xfig (const Triangulation& tria, - std::ostream& out, - const Mapping* /*mapping*/) + std::ostream& out, + const Mapping* /*mapping*/) { - 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) + 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) + // 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& 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::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& tria, unsigned int GridOut::n_boundary_faces (const Triangulation<1> &) const { - return 0; + return 0; } #endif @@ -537,16 +675,63 @@ unsigned int GridOut::n_boundary_faces (const Triangulation<1> &) const template unsigned int GridOut::n_boundary_faces (const Triangulation &tria) const { - typename Triangulation::active_face_iterator face, endf; - unsigned int n_faces = 0; + 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++; + 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 +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()) + << " 1 " << 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; + }; } @@ -554,10 +739,10 @@ unsigned int GridOut::n_boundary_faces (const Triangulation &tria) const #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 @@ -566,34 +751,34 @@ void GridOut::write_ucd_faces (const Triangulation<1> &, template void GridOut::write_ucd_faces (const Triangulation &tria, - const unsigned int starting_index, - std::ostream &out) const + const unsigned int starting_index, + std::ostream &out) const { - typename Triangulation::active_face_iterator face, endf; - unsigned int index=starting_index; + 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'; + 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; - }; + ++index; + }; } @@ -601,33 +786,33 @@ void GridOut::write_ucd_faces (const Triangulation &tria, #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(cell->material_id()) << '\n' - << cell->vertex(1) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << '\n'; - break; + 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'; + 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 @@ -636,124 +821,124 @@ void GridOut::write_gnuplot (const Triangulation<1> &tria, #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::active_cell_iterator cell=tria.begin_active(); - const 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) + 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::active_cell_iterator cell=tria.begin_active(); + const 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> quadrature(boundary_points, dummy_weights); + std::vector dummy_weights(n_points, 1./n_points); + Quadrature<1> quadrature(boundary_points, dummy_weights); - q_projector = new Quadrature (QProjector::project_to_all_faces(quadrature)); + q_projector = new Quadrature (QProjector::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::vertices_per_cell; ++i) - out << cell->vertex(GeometryInfo::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::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'; - out << cell->vertex(0) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << '\n' // double new line for gnuplot 3d plots - << '\n'; + << ' ' << 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 + 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) + for (unsigned int face_no=0; + face_no::faces_per_cell; ++face_no) { - const Triangulation::face_iterator - face = cell->face(face_no); - if (face->at_boundary()) + const Triangulation::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; itransform_unit_to_real_cell - (cell, q_projector->point(offset+i))) - << ' ' << cell->level() - << ' ' << static_cast(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; itransform_unit_to_real_cell + (cell, q_projector->point(offset+i))) + << ' ' << cell->level() + << ' ' << static_cast(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(cell->material_id()) - << '\n' - << face->vertex(1) - << ' ' << cell->level() - << ' ' << static_cast(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(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 (); + if (q_projector != 0) + delete q_projector; - AssertThrow (out, ExcIO()); + // make sure everything now gets to + // disk + out.flush (); + + AssertThrow (out, ExcIO()); } #endif @@ -762,207 +947,207 @@ void GridOut::write_gnuplot (const Triangulation<2> &tria, #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::active_cell_iterator cell=tria.begin_active(); - const 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) + 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::active_cell_iterator cell=tria.begin_active(); + const 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)); + 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) + + 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(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'; + // 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 + else { - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) { - const Triangulation::face_iterator - face = cell->face(face_no); - - if (face->at_boundary()) + const 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'; - } + 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 + else { - for (unsigned int l=0; l::lines_per_face; ++l) + for (unsigned int l=0; l::lines_per_face; ++l) { - const Triangulation::line_iterator - line=face->line(l); - - const Point &v0=line->vertex(0), - &v1=line->vertex(1); - if (line->at_boundary()) + const Triangulation::line_iterator + line=face->line(l); + + const Point &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 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)) + // transform_real_to_unit_cell + // could be + // replaced + // by using + // QProjector<3>::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'; - } - else - out << v0 - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << '\n' - << v1 - << ' ' << cell->level() - << ' ' << static_cast(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 @@ -974,493 +1159,497 @@ struct LineEntry 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 + template void GridOut::write_eps (const Triangulation &tria, - std::ostream &out, - const Mapping *mapping) + std::ostream &out, + const Mapping *mapping) { - - 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) + + 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_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, - // 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(), - 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 > 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) + case 1: + { + Assert(false, ExcInternalError()); + break; + }; + + case 2: + { + typename Triangulation::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, + // 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(), + 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::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 > 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 Point 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 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::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(), + face->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(), + 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::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 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 (; 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::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 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 (; 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(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)) + // 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 << ("/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::active_cell_iterator - cell = tria.begin_active (), - endc = tria.end (); - for (; cell!=endc; ++cell) + 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'; + 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 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 << "(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()); + out << "showpage" << '\n'; + + // make sure everything now gets to + // disk + out.flush (); + + AssertThrow (out, ExcIO()); } #endif -template + template void GridOut::write (const Triangulation &tria, - std::ostream &out, - const OutputFormat output_format, - const Mapping *mapping) + std::ostream &out, + const OutputFormat output_format, + const Mapping *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()); } diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index cc9d12813e..0f2d4d961f 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -426,6 +426,14 @@ inconvenience this causes.

deal.II

    +
  1. + New: There is now a GridOut::write_gmsh function to write in the gmsh + format. It has the same features of the write_ucd one. +
    + (Luca Heltai 2006/02/17) +

    +
  2. Changed: The hp::FECollection::get_fe,