]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Audit the code for the use of at_boundary, and check whether we should rather be...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 23 Feb 2001 09:41:03 +0000 (09:41 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 23 Feb 2001 09:41:03 +0000 (09:41 +0000)
git-svn-id: https://svn.dealii.org/trunk@4022 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/dofs/dof_handler.cc
deal.II/deal.II/source/dofs/dof_tools.cc
deal.II/deal.II/source/grid/tria.cc
deal.II/deal.II/source/numerics/vectors.cc

index b7d9817f2f4efc01c0e77435b5de0369254355af..635b550ac4a9e30763066b144432019a3a396e29 100644 (file)
@@ -1180,6 +1180,17 @@ unsigned int DoFHandler<dim>::n_boundary_dofs () const
 
   const unsigned int dofs_per_face = selected_fe->dofs_per_face;
   std::vector<unsigned int> 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)
index 7c206cda164b989ad7ea13d2530a85dad54bec65..af734185fa79ed66bae09a2f90e9e4c3fac3e342 100644 (file)
@@ -157,8 +157,19 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
 
   const unsigned int dofs_per_face = dof.get_fe().dofs_per_face;
   std::vector<unsigned int> dofs_on_this_face(dofs_per_face);
-  DoFHandler<dim>::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<dim>::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<dim> &dof,
   DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
                                        endc = dof.end();
 
-  (const_cast<Triangulation<dim>& > (dof.get_tria())).clear_user_flags();
+                                  // clear user flags for further use
+  const_cast<Triangulation<dim>&>(dof.get_tria()).clear_user_flags();
   
   for (; cell!=endc; ++cell)
     {
@@ -850,6 +862,20 @@ DoFTools::extract_boundary_dofs (const DoFHandler<dim>         &dof_handler,
   selected_dofs.clear ();
   selected_dofs.resize (dof_handler.n_dofs(), false);
   std::vector<unsigned int> 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<dim>::active_cell_iterator cell=dof_handler.begin_active();
        cell!=dof_handler.end(); ++cell)
     for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
@@ -1706,6 +1732,19 @@ void DoFTools::map_dof_to_boundary_indices (const DoFHandler<dim>     &dof_handl
   std::vector<unsigned int> 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<dim>::active_face_iterator face = dof_handler.begin_active_face(),
                                                 endf = dof_handler.end_face();
   for (; face!=endf; ++face)
index aada3a51403101517362ce560e803e1699309c9d..a0b5532ed52488e94d0329a792cb6f5cbd915915 100644 (file)
@@ -29,7 +29,7 @@
 
 
 template <int dim>
-const StraightBoundary<dim>& Triangulation<dim>::straight_boundary
+const StraightBoundary<dim> & Triangulation<dim>::straight_boundary
 = StraightBoundary<dim>();
 
 
@@ -1312,8 +1312,24 @@ void Triangulation<3>::create_triangulation (const std::vector<Point<3> >    &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<double>             minimal_length (vertices.size(),
-                                                 almost_infinite_length);
+  std::vector<double> minimal_length (vertices.size(),
+                                     almost_infinite_length);
                                   // also note if a vertex is at
                                   // the boundary
-  std::vector<bool>               at_boundary (vertices.size(), false);
+  std::vector<bool>   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 <int dim>
 void Triangulation<dim>::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<dim>::distort_random (const double factor,
   for (cell_iterator cell=begin(0); cell!=end(0); ++cell)
     almost_infinite_length += cell->diameter();
   
-  std::vector<double>             minimal_length (vertices.size(),
-                                                 almost_infinite_length);
+  std::vector<double> minimal_length (vertices.size(),
+                                     almost_infinite_length);
 
                                   // also note if a vertex is at
                                   // the boundary
-  std::vector<bool>               at_boundary (vertices.size(), false);
+  std::vector<bool>   at_boundary (vertices.size(), false);
   
   for (active_line_iterator line=begin_active_line();
        line != end_line(); ++line)
@@ -1683,8 +1700,10 @@ void Triangulation<dim>::distort_random (const double factor,
 };
 
 
+
 template <int dim>
-void Triangulation<dim>::set_all_refine_flags () {
+void Triangulation<dim>::set_all_refine_flags ()
+{
   active_cell_iterator cell = begin_active(),
                       endc = end();
 
@@ -1696,8 +1715,10 @@ void Triangulation<dim>::set_all_refine_flags () {
 };
 
 
+
 template <int dim>
-void Triangulation<dim>::refine_global (const unsigned int times) {
+void Triangulation<dim>::refine_global (const unsigned int times)
+{
   for (unsigned int i=0; i<times; ++i)
     {
       set_all_refine_flags();
@@ -1706,6 +1727,7 @@ void Triangulation<dim>::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
index 4e4198eaae4e9845826001cc8a781b0a21bb57d2..b99173ebc91da0f9eb1229e50f3ebb5f11009441 100644 (file)
@@ -381,12 +381,34 @@ void VectorTools::project (const DoFHandler<dim>    &dof,
   std::map<unsigned int,double> 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<dim>::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<dim>::active_face_iterator face = dof.begin_active_face(),
+                                                    endf = dof.end_face();
       std::vector<unsigned int> face_dof_indices (fe.dofs_per_face);
       for (; face!=endf; ++face)
        if (face->at_boundary())

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.