From: Ralf Hartmann Date: Mon, 20 Sep 2004 15:56:00 +0000 (+0000) Subject: Implementation for GridOut::write_gnuplot for curved boundaries in 3d. X-Git-Tag: v8.0.0~14822 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5b28dbc70eeb636a738e5cf8639fb987a7c60f36;p=dealii.git Implementation for GridOut::write_gnuplot for curved boundaries in 3d. git-svn-id: https://svn.dealii.org/trunk@9647 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/grid/grid_out.cc b/deal.II/deal.II/source/grid/grid_out.cc index 6a131f88a6..04a1e7f964 100644 --- a/deal.II/deal.II/source/grid/grid_out.cc +++ b/deal.II/deal.II/source/grid/grid_out.cc @@ -615,8 +615,8 @@ void GridOut::write_gnuplot (const Triangulation &tria, const unsigned int n_additional_points= gnuplot_flags.n_boundary_face_points; - const unsigned int n_points=GeometryInfo::vertices_per_face+n_additional_points; - + const unsigned int n_points=2+n_additional_points; + typename Triangulation::active_cell_iterator cell=tria.begin_active(); const typename Triangulation::active_cell_iterator endc=tria.end(); @@ -626,15 +626,21 @@ void GridOut::write_gnuplot (const Triangulation &tria, // used to probe boundary points at // curved faces Quadrature *q_projector=0; + std::vector > boundary_points; if (mapping!=0) { - std::vector > boundary_points(n_points); + boundary_points.resize(n_points); boundary_points[0][0]=0; boundary_points[n_points-1][0]=1; for (unsigned int i=1; i quadrature(boundary_points); + std::vector 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)); } @@ -767,22 +773,97 @@ void GridOut::write_gnuplot (const Triangulation &tria, case 3: { -//TODO:[RH] curved boundaries in 3d gnuplot not supported - Assert (mapping == 0, ExcNotImplemented()); - - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - { - const typename Triangulation::face_iterator - face = cell->face(face_no); - - for (unsigned int v=0; v::vertices_per_face; ++v) - out << face->vertex(v) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl; - out << std::endl; - out << std::endl; - } - break; + if (mapping==0 || n_points==2 || !cell->has_boundary_lines()) + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + const typename Triangulation::face_iterator + face = cell->face(face_no); + + out << "#face" << face_no << std::endl; + for (unsigned int v=0; v::vertices_per_face; ++v) + out << face->vertex(v) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl; + out << face->vertex(0) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl; + out << std::endl; + out << std::endl; + } + else + { + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + const typename Triangulation::face_iterator + face = cell->face(face_no); + + if (face->at_boundary()) + { + const unsigned int offset=face_no*n_points*n_points; + for (unsigned int i=0; i p0=mapping->transform_unit_to_real_cell( + cell, q_projector->point(offset+i*n_points+j)); + out << p0 + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl; + out << (mapping->transform_unit_to_real_cell( + cell, q_projector->point(offset+(i+1)*n_points+j))) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl; + 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()) << std::endl; + out << (mapping->transform_unit_to_real_cell( + cell, q_projector->point(offset+i*n_points+j+1))) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl; + out << p0 + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl; + out << std::endl << std::endl; + } + } + else + { + for (unsigned int l=0; l::lines_per_face; ++l) + { + const typename Triangulation::line_iterator + line=face->line(l); + + const Point &v0=line->vertex(0), + &v1=line->vertex(1); + if (line->at_boundary()) + { + // 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()) << std::endl; + } + else + out << v0 + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << v1 + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl; + out << std::endl << std::endl; + } + } + } + } + + break; }; }; };