From 3fd2e9386d0a7a6018337e0cf1b1b1788e5adb1d Mon Sep 17 00:00:00 2001 From: kayser-herold Date: Sat, 29 Jul 2006 19:30:46 +0000 Subject: [PATCH] Corrected a few errors and inserted a branch into do_make_hanging_node_constraints (3D), which goes into the new hp capable version of the code, when the FiniteElement supports it the new concept. git-svn-id: https://svn.dealii.org/trunk@13496 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/dofs/dof_accessor.templates.h | 20 +++++++- deal.II/deal.II/include/dofs/hp_dof_objects.h | 1 + deal.II/deal.II/source/dofs/dof_tools.cc | 46 +++++++++++++------ deal.II/deal.II/source/fe/fe_q.cc | 33 +++++++++---- deal.II/examples/step-3/step-3.cc | 5 +- 5 files changed, 79 insertions(+), 26 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_accessor.templates.h b/deal.II/deal.II/include/dofs/dof_accessor.templates.h index d65f16ad4b..8b16e784ac 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.templates.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.templates.h @@ -297,11 +297,16 @@ DoFObjectAccessor<1,DH>::get_dof_indices (std::vector &dof_indices // freedom on sub-objects which are // not allocated for this // non-active thing +// TODO: This assertion leads to problems with the do_make_hp_hanging_node_constraints +// method!. Check how this assertion can be stated that it does not create a conflict +// with the above mentioned method. +/* Assert (!this->has_children() || (this->dof_handler->get_fe()[fe_index].dofs_per_cell == 2*this->dof_handler->get_fe()[fe_index].dofs_per_vertex), typename BaseClass::ExcNotActive()); - +*/ + const unsigned int dofs_per_vertex = this->dof_handler->get_fe()[fe_index].dofs_per_vertex, dofs_per_line = this->dof_handler->get_fe()[fe_index].dofs_per_line; std::vector::iterator next = dof_indices.begin(); @@ -373,11 +378,17 @@ DoFObjectAccessor<2,DH>::get_dof_indices (std::vector &dof_indices // freedom on sub-objects which are // not allocated for this // non-active thing + +// TODO: This assertion leads to problems with the do_make_hp_hanging_node_constraints +// method!. Check how this assertion can be stated that it does not create a conflict +// with the above mentioned method. +/* Assert (!this->has_children() || (this->dof_handler->get_fe()[fe_index].dofs_per_cell == 4*this->dof_handler->get_fe()[fe_index].dofs_per_vertex), typename BaseClass::ExcNotActive()); - +*/ + Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), ExcMessage ("DoFHandler not initialized")); @@ -472,10 +483,15 @@ DoFObjectAccessor<3,DH>::get_dof_indices (std::vector &dof_indices // freedom on sub-objects which are // not allocated for this // non-active thing +// TODO: This assertion leads to problems with the do_make_hp_hanging_node_constraints +// method!. Check how this assertion can be stated that it does not create a conflict +// with the above mentioned method. +/* Assert (!this->has_children() || (this->dof_handler->get_fe()[fe_index].dofs_per_cell == 8*this->dof_handler->get_fe()[fe_index].dofs_per_vertex), typename BaseClass::ExcNotActive()); +*/ Assert (static_cast(this->present_level) < this->dof_handler->levels.size(), ExcMessage ("DoFHandler not initialized")); diff --git a/deal.II/deal.II/include/dofs/hp_dof_objects.h b/deal.II/deal.II/include/dofs/hp_dof_objects.h index b4be1e955b..216a9e1426 100644 --- a/deal.II/deal.II/include/dofs/hp_dof_objects.h +++ b/deal.II/deal.II/include/dofs/hp_dof_objects.h @@ -15,6 +15,7 @@ #include #include +#include #include diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 6933b0f9c7..c603fd779d 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -1584,7 +1584,7 @@ namespace internal Assert (n_dofs_mother <= n_dofs_child, ExcInternalError ()); - + for (unsigned int col=0; col!=n_dofs_child; ++col) { bool constrain = true; @@ -1593,11 +1593,10 @@ namespace internal // which is already satisfied. for (unsigned int i=0; i!=n_dofs_mother; ++i) { - // Check for a value of almost exactly + // Check for a value of // 1.0 -//TODO: move this check into the FE classes and let them reset the values to exactly 1.0. similarly, let them make sure that very small values are exactly 0.0 - if (fabs (face_constraints (i,col) - 1.0) < 1.0e-15) - constrain = (dofs_mother[i] != dofs_child[col]); + if (face_constraints (i,col) == 1.0) + constrain = (dofs_mother[i] != dofs_child[col]); } // If this constraint is not @@ -1666,10 +1665,9 @@ namespace internal if (dof_handler.get_fe().hp_constraints_are_implemented ()) { do_make_hp_hanging_node_constraints (dof_handler, - constraints, - ::internal::int2type()); - return; - } + constraints); + return; + } const unsigned int dim = 2; @@ -1806,16 +1804,16 @@ namespace internal ExcInternalError()); } } +#endif - +#if deal_II_dimension > 1 template static void do_make_hp_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints, - ::internal::int2type<2>) + ConstraintMatrix &constraints) { - const unsigned int dim = 2; + const unsigned int dim = DH::dimension; std::vector dofs_on_mother; std::vector dofs_on_children; @@ -1961,7 +1959,12 @@ namespace internal Assert (false, ExcNotImplemented ()); // TODO: That's the difficult one. - } +// Sketch of how this has to be done: +// The coarse element is constrained to a lower order element with +// the degree of the lowest order element. +// Afterwards the two finer elements are constrained to the this +// constrained element. + } } else { @@ -2070,6 +2073,21 @@ namespace internal ConstraintMatrix &constraints, ::internal::int2type<3>) { + // Decide whether to use the + // new or old make_hanging_node_constraints + // function. If all the FiniteElement + // or all elements in a FECollection support + // the new face constraint matrix, the + // new code will be used. + // Otherwise, the old implementation is used + // for the moment. + if (dof_handler.get_fe().hp_constraints_are_implemented ()) + { + do_make_hp_hanging_node_constraints (dof_handler, + constraints); + return; + } + const unsigned int dim = 3; std::vector dofs_on_mother; diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 34bfa21f0d..08c2202a7f 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -400,8 +400,14 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, // which returns the support // points on the face. Quadrature quad_face_support (source_fe.get_unit_face_support_points ()); - - + + // Rule of thumb for FP accuracy, + // that can be expected for a + // given polynomial degree. + // This value is used to cut + // off values close to zero. + double eps = 2e-14*this->degree*(dim-1); + // compute the interpolation // matrix by simply taking the // value at the support points. @@ -413,15 +419,24 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, Point p = QProjector::project_to_face (quad_face_support, 0).point (i); for (unsigned int j=0; jdofs_per_face; ++j) - interpolation_matrix(j,i) = this->shape_value (this->face_to_cell_index(j, 0), p); + { + double matrix_entry = this->shape_value (this->face_to_cell_index(j, 0), p); + + // Correct the interpolated + // value. I.e. if it is close + // to 1 or 0, make it exactly + // 1 or 0. Unfortunately, this + // is required to avoid problems + // with higher order elements. + if (fabs (matrix_entry - 1.0) < eps) + matrix_entry = 1.0; + if (fabs (matrix_entry) < eps) + matrix_entry = 0.0; + + interpolation_matrix(j,i) = matrix_entry; + } } - // cut off very small values - for (unsigned int i=0; idofs_per_face; ++i) - for (unsigned int j=0; j