]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix a problem where the hp::DoFHandler got confused with a new triangulation.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 24 Apr 2011 04:28:19 +0000 (04:28 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 24 Apr 2011 04:28:19 +0000 (04:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@23640 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/examples/step-46/step-46.cc
deal.II/include/deal.II/grid/tria.h
deal.II/include/deal.II/hp/dof_handler.h
deal.II/source/grid/tria.cc
deal.II/source/hp/dof_handler.cc
tests/hp/crash_20.cc [new file with mode: 0644]
tests/hp/crash_20/cmp/generic [new file with mode: 0644]

index f8844005cbe9ff3335260bd497f685b4e9632dbb..79ff92045d10b6fd85b3196701c87d9c4b0f9146 100644 (file)
@@ -85,6 +85,15 @@ should be fixed now.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: If Triangulation::create_triangulation is called after an
+hp::DoFHandler object is attached to the triangulation object, setting active
+FE indices leads to a crash. The problem did not happen if the mesh was
+refined before setting the FE indices. This is now fixed. In the process, the
+Triangulation::RefinementListener::create_notification function was
+introduced.
+<br>
+(Wolfgang Bangerth, 2011/04/22)
+
 <li> Fixed: The function FEValuesViews::SymmetricTensor::divergence had a bug.
 This is now fixed.
 <br>
index 2799e9059fe1bee6991e0f477b330a76aaad2e71..eef01a25d515c8a2580772df707ea8a074d331a8 100644 (file)
@@ -295,7 +295,7 @@ void
 FluidStructureProblem<dim>::make_grid ()
 {
 // not quite what we want...  
-  GridGenerator::subdivided_hyper_cube (triangulation, 4, -1, 1);
+  GridGenerator::subdivided_hyper_cube (triangulation, 8, -1, 1);
   for (typename Triangulation<dim>::active_cell_iterator
         cell = triangulation.begin_active();
        cell != triangulation.end(); ++cell)
@@ -305,8 +305,6 @@ FluidStructureProblem<dim>::make_grid ()
          (cell->face(f)->center()[dim-1] == 1))
        cell->face(f)->set_all_boundary_indicators(1);
 
-  triangulation.refine_global (3-dim);
-
 
   for (typename Triangulation<dim>::active_cell_iterator
          cell = dof_handler.begin_active();
index 0a9f6d497bd4f26eae5e1fbb4066d8033e3d61ce..2a9651ddd2f7a6a0295a19780740d7bc34d67a8c 100644 (file)
@@ -1374,6 +1374,21 @@ class Triangulation : public Subscriptor
         void
         copy_notification (const Triangulation<dim, spacedim> &old_tria,
                           const Triangulation<dim, spacedim> &new_tria);
+
+                                         /**
+                                          * At the end of a call to
+                                          * create_triangulation() the
+                                          * Triangulation class calls this
+                                          * method on all objects derived from
+                                          * this class and registered with the
+                                          * current Triangulation object. By
+                                          * default this method does nothing,
+                                          * a different behavior has to be
+                                          * implemented in derived classes.
+                                          */
+        virtual
+        void
+        create_notification (const Triangulation<dim, spacedim> &tria);
     };
 
                                     /**
index 225fe63489ef5e70bed182c6afb9d95ac20ffe96..53b7958e10cb1ae90f13430662994515633e7eaf 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
+//    Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -1135,6 +1135,8 @@ namespace hp
                                         */
       virtual void pre_refinement_notification (const Triangulation<dim,spacedim> &tria);
       virtual void post_refinement_notification (const Triangulation<dim,spacedim> &tria);
+      virtual void create_notification (const Triangulation<dim,spacedim> &tria);
+      
 
                                       /**
                                        * Compute identities between
index a58c2ea95dfbe20a27724f09dd173507bd3b8ba7..2c853c7328e7d5f207ab1f602a50caa0cd7e64e7 100644 (file)
@@ -9592,7 +9592,7 @@ create_triangulation (const std::vector<Point<spacedim> >    &v,
     }
 
   compute_number_cache (*this, levels.size(), number_cache);
-
+  
                                   // now verify that there are indeed
                                   // no distorted cells. as per the
                                   // documentation of this class, we
@@ -9640,98 +9640,106 @@ create_triangulation (const std::vector<Point<spacedim> >    &v,
     much bigger than the minimal data required, but it makes the code more readable.
 
   */
-  if (dim < spacedim) {
-
-    Table<2,bool> correct(GeometryInfo< dim >::faces_per_cell,
-                         GeometryInfo< dim >::faces_per_cell);
-    switch (dim)
-      {
-       case 1:
-       {
-         bool values [][2] = {{false,true},
-                              {true,false} };
-         for (unsigned int i=0; i< GeometryInfo< dim >::faces_per_cell; ++i)
-           for (unsigned int j=0; j< GeometryInfo< dim >::faces_per_cell; ++j)
-             correct(i,j) = ( values[i][j]);
-         break;
-       }
-       case 2:
+  if (dim < spacedim)
+    {
+      Table<2,bool> correct(GeometryInfo< dim >::faces_per_cell,
+                           GeometryInfo< dim >::faces_per_cell);
+      switch (dim)
        {
-         bool values [][4]= {{false,true ,true , false},
-                             {true ,false,false, true },
-                             {true ,false,false, true },
-                             {false,true ,true , false} };
-         for (unsigned int i=0; i< GeometryInfo< dim >::faces_per_cell; ++i)
-           for (unsigned int j=0; j< GeometryInfo< dim >::faces_per_cell; ++j)
-             correct(i,j) = ( values[i][j]);
-         break;
+         case 1:
+         {
+           bool values [][2] = {{false,true},
+                                {true,false} };
+           for (unsigned int i=0; i< GeometryInfo< dim >::faces_per_cell; ++i)
+             for (unsigned int j=0; j< GeometryInfo< dim >::faces_per_cell; ++j)
+               correct(i,j) = ( values[i][j]);
+           break;
+         }
+         case 2:
+         {
+           bool values [][4]= {{false,true ,true , false},
+                               {true ,false,false, true },
+                               {true ,false,false, true },
+                               {false,true ,true , false} };
+           for (unsigned int i=0; i< GeometryInfo< dim >::faces_per_cell; ++i)
+             for (unsigned int j=0; j< GeometryInfo< dim >::faces_per_cell; ++j)
+               correct(i,j) = ( values[i][j]);
+           break;
+         }
+         default:
+               Assert (false, ExcNotImplemented());
        }
-       default:
-             Assert (false, ExcNotImplemented());
-      }
 
 
-    std::list<active_cell_iterator> this_round, next_round;
-    active_cell_iterator neighbor;
+      std::list<active_cell_iterator> this_round, next_round;
+      active_cell_iterator neighbor;
 
-    this_round.push_back (begin_active());
-    begin_active()->set_direction_flag (true);
-    begin_active()->set_user_flag ();
+      this_round.push_back (begin_active());
+      begin_active()->set_direction_flag (true);
+      begin_active()->set_user_flag ();
 
-    while (this_round.size() > 0)
-      {
-       for ( typename std::list<active_cell_iterator>::iterator cell = this_round.begin();
-             cell != this_round.end(); ++cell)
-         {
-           for (unsigned int i = 0; i < GeometryInfo< dim >::faces_per_cell; ++i)
-             {
-               if ( !((*cell)->face(i)->at_boundary()) )
-                 {
-                   neighbor = (*cell)->neighbor(i);
+      while (this_round.size() > 0)
+       {
+         for ( typename std::list<active_cell_iterator>::iterator cell = this_round.begin();
+               cell != this_round.end(); ++cell)
+           {
+             for (unsigned int i = 0; i < GeometryInfo< dim >::faces_per_cell; ++i)
+               {
+                 if ( !((*cell)->face(i)->at_boundary()) )
+                   {
+                     neighbor = (*cell)->neighbor(i);
 
-                   unsigned int cf = (*cell)->face_index(i);
-                   unsigned int j = 0;
-                   while(neighbor->face_index(j) != cf)
-                     {++j;}
+                     unsigned int cf = (*cell)->face_index(i);
+                     unsigned int j = 0;
+                     while(neighbor->face_index(j) != cf)
+                       {++j;}
 
-                   if ( (correct(i,j) && !(*cell)->direction_flag())
-                        ||
-                        (!correct(i,j) && (*cell)->direction_flag()) )
-                     {
-                       if (neighbor->user_flag_set() == false)
-                         {
-                           neighbor->set_direction_flag (false);
-                           neighbor->set_user_flag ();
-                           next_round.push_back (neighbor);
-                         }
-                       else
-                         Assert (neighbor->direction_flag() == false,
-                                 ExcNonOrientableTriangulation());
-
-                     }
-                 }
-             }
-         }
+                     if ( (correct(i,j) && !(*cell)->direction_flag())
+                          ||
+                          (!correct(i,j) && (*cell)->direction_flag()) )
+                       {
+                         if (neighbor->user_flag_set() == false)
+                           {
+                             neighbor->set_direction_flag (false);
+                             neighbor->set_user_flag ();
+                             next_round.push_back (neighbor);
+                           }
+                         else
+                           Assert (neighbor->direction_flag() == false,
+                                   ExcNonOrientableTriangulation());
 
-                                        // Before we quit let's check
-                                        // that if the triangulation
-                                        // is disconnected that we
-                                        // still get all cells
-      if (next_round.size() == 0)
-       for (active_cell_iterator cell = begin_active();
-            cell != end(); ++cell)
-         if (cell->user_flag_set() == false)
-           {
-             next_round.push_back (cell);
-             cell->set_direction_flag (true);
-             cell->set_user_flag ();
-             break;
+                       }
+                   }
+               }
            }
 
-      this_round = next_round;
-      next_round.clear();
+                                          // Before we quit let's check
+                                          // that if the triangulation
+                                          // is disconnected that we
+                                          // still get all cells
+         if (next_round.size() == 0)
+           for (active_cell_iterator cell = begin_active();
+                cell != end(); ++cell)
+             if (cell->user_flag_set() == false)
+               {
+                 next_round.push_back (cell);
+                 cell->set_direction_flag (true);
+                 cell->set_user_flag ();
+                 break;
+               }
+
+         this_round = next_round;
+         next_round.clear();
+       }
     }
-  }
+
+                                  // inform all listeners that the
+                                  // triangulation has been created
+  typename std::list<RefinementListener *>::iterator
+    ref_listener = refinement_listeners.begin (),
+    end_listener = refinement_listeners.end ();
+  for (; ref_listener != end_listener; ++ref_listener)
+    (*ref_listener)->create_notification (*this);
 }
 
 
@@ -14496,6 +14504,13 @@ RefinementListener::copy_notification (const Triangulation<dim, spacedim> &,
 
 
 
+template<int dim, int spacedim>
+void Triangulation<dim, spacedim>::
+RefinementListener::create_notification (const Triangulation<dim, spacedim> &)
+{}
+
+
+
 template<int dim, int spacedim>
 void
 Triangulation<dim, spacedim>::add_refinement_listener (RefinementListener &listener) const
index 83ed306ef6d97024d27305cc0784e19fb239b1d2..dce829ba4dbc1cf978d2cc58bbf2aa78ad6b6a40 100644 (file)
@@ -3984,7 +3984,9 @@ namespace hp
 
 
   template<int dim, int spacedim>
-  void DoFHandler<dim,spacedim>::post_refinement_notification (const Triangulation<dim,spacedim> &tria)
+  void
+  DoFHandler<dim,spacedim>::
+  post_refinement_notification (const Triangulation<dim,spacedim> &tria)
   {
     Assert (has_children.size () == levels.size (), ExcInternalError ());
 
@@ -4068,6 +4070,18 @@ namespace hp
   }
 
 
+  template<int dim, int spacedim>
+  void
+  DoFHandler<dim,spacedim>::
+  create_notification (const Triangulation<dim,spacedim> &tria)
+  {
+                                    // do the same here as needs to be done
+                                    // in the post_refinement hook
+    post_refinement_notification (tria);
+  }
+
+  
+
   template<int dim, int spacedim>
   void DoFHandler<dim,spacedim>::clear_space ()
   {
diff --git a/tests/hp/crash_20.cc b/tests/hp/crash_20.cc
new file mode 100644 (file)
index 0000000..c060264
--- /dev/null
@@ -0,0 +1,71 @@
+//----------------------------  crash_20.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 2006, 2011 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  crash_20.cc  ---------------------------
+
+
+// if the mesh is generated after the hp::DoFHandler is attached to the
+// triangulation object, then we can't set active fe indices -- which is
+// somewhat tragic since we have to assign active fe indices before we can
+// call distribute_dofs
+//
+// originally, this problem was avoided because the hp::DoFHandler listens to
+// the refinement listener signal to rebuild its data structures; so if you
+// create a triangulation object, attach the hp::DoFHandler, create a coarse
+// mesh, then refine the mesh, everything is ok again. the solution is to also
+// listen to the creation of triangulations.
+
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <hp/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_dgq.h>
+
+#include <fstream>
+
+
+
+template <int dim>
+void test ()
+{
+  Triangulation<dim> tria;
+  hp::DoFHandler<dim> dof_handler(tria);
+
+  GridGenerator::hyper_cube(tria);
+
+  for (typename hp::DoFHandler<dim>::active_cell_iterator
+         cell=dof_handler.begin_active();
+       cell!=dof_handler.end(); ++cell)
+    cell->set_active_fe_index (0);
+}
+
+
+int main ()
+{
+  std::ofstream logfile("crash_20/output");
+  logfile.precision(2);
+  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);  
+
+  test<1> ();
+  test<2> ();
+  test<3> ();
+  
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/hp/crash_20/cmp/generic b/tests/hp/crash_20/cmp/generic
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

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.