From ee5992e3c2c1550cad39660642475171a49a5fe0 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 23 Feb 2001 09:41:03 +0000 Subject: [PATCH] Audit the code for the use of at_boundary, and check whether we should rather be using has_boundary_lines. Add comments at various places. git-svn-id: https://svn.dealii.org/trunk@4022 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/dofs/dof_handler.cc | 11 +++ deal.II/deal.II/source/dofs/dof_tools.cc | 45 ++++++++++- deal.II/deal.II/source/grid/tria.cc | 86 ++++++++++++++++++---- deal.II/deal.II/source/numerics/vectors.cc | 30 +++++++- 4 files changed, 151 insertions(+), 21 deletions(-) diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index b7d9817f2f..635b550ac4 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -1180,6 +1180,17 @@ unsigned int DoFHandler::n_boundary_dofs () const const unsigned int dofs_per_face = selected_fe->dofs_per_face; std::vector dofs_on_face(dofs_per_face); + + // loop over all faces to check + // whether they are at a + // boundary. note that we need not + // take special care of single + // lines (using + // @p{cell->has_boundary_lines}), + // since we do not support + // boundaries of dimension dim-2, + // and so every boundary line is + // also part of a boundary face. active_face_iterator face = begin_active_face (), endf = end_face(); for (; face!=endf; ++face) diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 7c206cda16..af734185fa 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -157,8 +157,19 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, const unsigned int dofs_per_face = dof.get_fe().dofs_per_face; std::vector dofs_on_this_face(dofs_per_face); - DoFHandler::active_face_iterator face = dof.begin_active_face(), - endf = dof.end_face(); + + // loop over all faces to check + // whether they are at a + // boundary. note that we need not + // take special care of single + // lines (using + // @p{cell->has_boundary_lines}), + // since we do not support + // boundaries of dimension dim-2, + // and so every boundary line is + // also part of a boundary face. + typename DoFHandler::active_face_iterator face = dof.begin_active_face(), + endf = dof.end_face(); for (; face!=endf; ++face) if (face->at_boundary()) { @@ -250,7 +261,8 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler &dof, DoFHandler::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); - (const_cast& > (dof.get_tria())).clear_user_flags(); + // clear user flags for further use + const_cast&>(dof.get_tria()).clear_user_flags(); for (; cell!=endc; ++cell) { @@ -850,6 +862,20 @@ DoFTools::extract_boundary_dofs (const DoFHandler &dof_handler, selected_dofs.clear (); selected_dofs.resize (dof_handler.n_dofs(), false); std::vector face_dof_indices (dof_handler.get_fe().dofs_per_face); + + // now loop over all cells and + // check whether their faces are at + // the boundary. note that we need + // not take special care of single + // lines being at the boundary + // (using + // @p{cell->has_boundary_lines}), + // since we do not support + // boundaries of dimension dim-2, + // and so every isolated boundary + // line is also part of a boundary + // face which we will be visiting + // sooner or later for (DoFHandler::active_cell_iterator cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell) for (unsigned int face=0; face::faces_per_cell; ++face) @@ -1706,6 +1732,19 @@ void DoFTools::map_dof_to_boundary_indices (const DoFHandler &dof_handl std::vector dofs_on_face(dofs_per_face); unsigned int next_boundary_index = 0; + // now loop over all cells and + // check whether their faces are at + // the boundary. note that we need + // not take special care of single + // lines being at the boundary + // (using + // @p{cell->has_boundary_lines}), + // since we do not support + // boundaries of dimension dim-2, + // and so every isolated boundary + // line is also part of a boundary + // face which we will be visiting + // sooner or later typename DoFHandler::active_face_iterator face = dof_handler.begin_active_face(), endf = dof_handler.end_face(); for (; face!=endf; ++face) diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index aada3a5140..a0b5532ed5 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -29,7 +29,7 @@ template -const StraightBoundary& Triangulation::straight_boundary +const StraightBoundary & Triangulation::straight_boundary = StraightBoundary(); @@ -1312,8 +1312,24 @@ void Triangulation<3>::create_triangulation (const std::vector > &v, // interior for (line_iterator line=begin_line(); line!=end_line(); ++line) line->set_boundary_indicator (255); - // next reset all lines bounding boundary - // quads as on the boundary also + // next reset all lines bounding + // boundary quads as on the + // boundary also. note that since + // we are in 3d, there are cases + // where one or more lines of a + // quad that is not on the + // boundary, are actually boundary + // lines. they will not be marked + // when visiting this + // face. however, since we do not + // support dim-2 dimensional + // boundaries (i.e. internal lines + // constituting boundaries), every + // such line is also part of a face + // that is actually on the + // boundary, so sooner or later we + // get to mark that line for being + // on the boundary for (quad_iterator quad=begin_quad(); quad!=end_quad(); ++quad) if (quad->at_boundary()) for (unsigned int l=0; l<4; ++l) @@ -1533,11 +1549,11 @@ void Triangulation<1>::distort_random (const double factor, for (cell_iterator cell=begin(0); cell!=end(0); ++cell) almost_infinite_length += cell->diameter(); - std::vector minimal_length (vertices.size(), - almost_infinite_length); + std::vector minimal_length (vertices.size(), + almost_infinite_length); // also note if a vertex is at // the boundary - std::vector at_boundary (vertices.size(), false); + std::vector at_boundary (vertices.size(), false); for (active_line_iterator line=begin_active_line(); line != end_line(); ++line) @@ -1585,7 +1601,8 @@ void Triangulation<1>::distort_random (const double factor, template void Triangulation::distort_random (const double factor, - const bool keep_boundary) { + const bool keep_boundary) +{ // this function is mostly equivalent to // that for the general dimensional case // the only difference being the correction @@ -1607,12 +1624,12 @@ void Triangulation::distort_random (const double factor, for (cell_iterator cell=begin(0); cell!=end(0); ++cell) almost_infinite_length += cell->diameter(); - std::vector minimal_length (vertices.size(), - almost_infinite_length); + std::vector minimal_length (vertices.size(), + almost_infinite_length); // also note if a vertex is at // the boundary - std::vector at_boundary (vertices.size(), false); + std::vector at_boundary (vertices.size(), false); for (active_line_iterator line=begin_active_line(); line != end_line(); ++line) @@ -1683,8 +1700,10 @@ void Triangulation::distort_random (const double factor, }; + template -void Triangulation::set_all_refine_flags () { +void Triangulation::set_all_refine_flags () +{ active_cell_iterator cell = begin_active(), endc = end(); @@ -1696,8 +1715,10 @@ void Triangulation::set_all_refine_flags () { }; + template -void Triangulation::refine_global (const unsigned int times) { +void Triangulation::refine_global (const unsigned int times) +{ for (unsigned int i=0; i::refine_global (const unsigned int times) { }; + // if necessary try to work around a bug in the IBM xlC compiler #ifdef XLC_WORK_AROUND_STD_BUG using namespace std; @@ -4515,7 +4537,8 @@ void Triangulation<2>::execute_refinement () { #if deal_II_dimension == 3 template <> -void Triangulation<3>::execute_refinement () { +void Triangulation<3>::execute_refinement () +{ const unsigned int dim = 3; // check whether a new level is needed @@ -4818,9 +4841,44 @@ void Triangulation<3>::execute_refinement () { vertices[next_unused_vertex] = boundary[quad->boundary_indicator()]->get_new_point_on_quad (quad); else + // it might be that the + // quad itself is not + // at the boundary, but + // that one of its line + // actually is. in this + // case, the newly + // created vertices at + // the centers of the + // lines are not + // necessarily the mean + // values of the + // adjacent vertices, + // so do not compute + // the new vertex as + // the mean value of + // the 4 vertices of + // the face, but rather + // as the mean value of + // the 8 vertices which + // we already have (the + // four old ones, and + // the four ones + // inserted as middle + // points for the four + // lines). summing up + // some more points is + // generally cheaper + // than first asking + // whether one of the + // lines is at the + // boundary vertices[next_unused_vertex] = (quad->vertex(0) + quad->vertex(1) + - quad->vertex(2) + quad->vertex(3)) / 4; + quad->vertex(2) + quad->vertex(3) + + quad->line(0)->child(0)->vertex(1) + + quad->line(1)->child(0)->vertex(1) + + quad->line(2)->child(0)->vertex(1) + + quad->line(3)->child(0)->vertex(1) ) / 8; // now that we created the right // point, make up the four diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 4e4198eaae..b99173ebc9 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -381,12 +381,34 @@ void VectorTools::project (const DoFHandler &dof, std::map boundary_values; if (enforce_zero_boundary == true) - // no need to project boundary values, but - // enforce homogeneous boundary values + // no need to project boundary + // values, but enforce + // homogeneous boundary values // anyway { - DoFHandler::active_face_iterator face = dof.begin_active_face(), - endf = dof.end_face(); + // loop over all boundary faces + // to get all dof indices of + // dofs on the boundary. note + // that in 3d there are cases + // where a face is not at the + // boundary, yet one of its + // lines is, and we should + // consider the degrees of + // freedom on it as boundary + // nodes. likewise, in 2d and + // 3d there are cases where a + // cell is only at the boundary + // by one vertex. nevertheless, + // since we do not support + // boundaries with dimension + // less or equal to dim-2, each + // such boundary dof is also + // found from some other face + // that is actually wholly on + // the boundary, not only by + // one line or one vertex + typename DoFHandler::active_face_iterator face = dof.begin_active_face(), + endf = dof.end_face(); std::vector face_dof_indices (fe.dofs_per_face); for (; face!=endf; ++face) if (face->at_boundary()) -- 2.39.5