From: Wolfgang Bangerth Date: Sun, 24 Apr 2011 04:28:19 +0000 (+0000) Subject: Fix a problem where the hp::DoFHandler got confused with a new triangulation. X-Git-Tag: v8.0.0~4117 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7ac0fe1b34ececb8f04f0864718d5c31454149a7;p=dealii.git Fix a problem where the hp::DoFHandler got confused with a new triangulation. git-svn-id: https://svn.dealii.org/trunk@23640 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index f8844005cb..79ff92045d 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -85,6 +85,15 @@ should be fixed now.

Specific improvements

    +
  1. 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. +
    +(Wolfgang Bangerth, 2011/04/22) +
  2. Fixed: The function FEValuesViews::SymmetricTensor::divergence had a bug. This is now fixed.
    diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc index 2799e9059f..eef01a25d5 100644 --- a/deal.II/examples/step-46/step-46.cc +++ b/deal.II/examples/step-46/step-46.cc @@ -295,7 +295,7 @@ void FluidStructureProblem::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::active_cell_iterator cell = triangulation.begin_active(); cell != triangulation.end(); ++cell) @@ -305,8 +305,6 @@ FluidStructureProblem::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::active_cell_iterator cell = dof_handler.begin_active(); diff --git a/deal.II/include/deal.II/grid/tria.h b/deal.II/include/deal.II/grid/tria.h index 0a9f6d497b..2a9651ddd2 100644 --- a/deal.II/include/deal.II/grid/tria.h +++ b/deal.II/include/deal.II/grid/tria.h @@ -1374,6 +1374,21 @@ class Triangulation : public Subscriptor void copy_notification (const Triangulation &old_tria, const Triangulation &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 &tria); }; /** diff --git a/deal.II/include/deal.II/hp/dof_handler.h b/deal.II/include/deal.II/hp/dof_handler.h index 225fe63489..53b7958e10 100644 --- a/deal.II/include/deal.II/hp/dof_handler.h +++ b/deal.II/include/deal.II/hp/dof_handler.h @@ -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 &tria); virtual void post_refinement_notification (const Triangulation &tria); + virtual void create_notification (const Triangulation &tria); + /** * Compute identities between diff --git a/deal.II/source/grid/tria.cc b/deal.II/source/grid/tria.cc index a58c2ea95d..2c853c7328 100644 --- a/deal.II/source/grid/tria.cc +++ b/deal.II/source/grid/tria.cc @@ -9592,7 +9592,7 @@ create_triangulation (const std::vector > &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 > &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 this_round, next_round; - active_cell_iterator neighbor; + std::list 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::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::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::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 &, +template +void Triangulation:: +RefinementListener::create_notification (const Triangulation &) +{} + + + template void Triangulation::add_refinement_listener (RefinementListener &listener) const diff --git a/deal.II/source/hp/dof_handler.cc b/deal.II/source/hp/dof_handler.cc index 83ed306ef6..dce829ba4d 100644 --- a/deal.II/source/hp/dof_handler.cc +++ b/deal.II/source/hp/dof_handler.cc @@ -3984,7 +3984,9 @@ namespace hp template - void DoFHandler::post_refinement_notification (const Triangulation &tria) + void + DoFHandler:: + post_refinement_notification (const Triangulation &tria) { Assert (has_children.size () == levels.size (), ExcInternalError ()); @@ -4068,6 +4070,18 @@ namespace hp } + template + void + DoFHandler:: + create_notification (const Triangulation &tria) + { + // do the same here as needs to be done + // in the post_refinement hook + post_refinement_notification (tria); + } + + + template void DoFHandler::clear_space () { diff --git a/tests/hp/crash_20.cc b/tests/hp/crash_20.cc new file mode 100644 index 0000000000..c060264599 --- /dev/null +++ b/tests/hp/crash_20.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + +template +void test () +{ + Triangulation tria; + hp::DoFHandler dof_handler(tria); + + GridGenerator::hyper_cube(tria); + + for (typename hp::DoFHandler::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 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/hp/crash_20/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK