From 77f54e9142b7d5f451d17935574803ec55a3e281 Mon Sep 17 00:00:00 2001 From: wolf Date: Thu, 16 Oct 2003 14:27:33 +0000 Subject: [PATCH] Fix various things for 3d with misoriented faces. git-svn-id: https://svn.dealii.org/trunk@8095 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/tria_accessor.h | 106 +++++++++++------- .../include/grid/tria_accessor.templates.h | 49 +++++++- deal.II/deal.II/source/fe/mapping_q1.cc | 12 +- deal.II/deal.II/source/grid/tria.cc | 18 +-- deal.II/deal.II/source/grid/tria_accessor.cc | 44 +------- .../source/numerics/error_estimator.cc | 104 +++++++++++------ 6 files changed, 202 insertions(+), 131 deletions(-) diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h index 9249aec2ac..38fa731872 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.h +++ b/deal.II/deal.II/include/grid/tria_accessor.h @@ -676,20 +676,28 @@ class TriaObjectAccessor : public TriaAccessor unsigned int number_of_children () const; /** - * Enquire some information about - * whether the face with the - * given number is in standard - * orientation or not. This - * information is only useful - * (and available) in 3d, see the + * Return whether the face with + * index @p{face} has its normal + * pointing in the standard + * direction (@p{true}) or + * whether it is the opposite + * (@p{false}). Which is the + * standard direction is + * documented with the + * @ref{Triangulation} class. In + * 1d and 2d, this is always + * @p{true}, but in 3d it may be + * different, see the respective + * discussion in the * documentation of the - * respective class. For all - * other dimensions, this - * function should not be called - * and triggers and assertion if - * done so. + * @ref{Triangulation} class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. */ - bool get_face_orientation (const unsigned int face) const; + bool face_orientation (const unsigned int face) const; private: /** @@ -1135,20 +1143,28 @@ class TriaObjectAccessor<1, dim> : public TriaAccessor unsigned int number_of_children () const; /** - * Enquire some information about - * whether the face with the - * given number is in standard - * orientation or not. This - * information is only useful - * (and available) in 3d, see the + * Return whether the face with + * index @p{face} has its normal + * pointing in the standard + * direction (@p{true}) or + * whether it is the opposite + * (@p{false}). Which is the + * standard direction is + * documented with the + * @ref{Triangulation} class. In + * 1d and 2d, this is always + * @p{true}, but in 3d it may be + * different, see the respective + * discussion in the * documentation of the - * respective class. For all - * other dimensions, this - * function should not be called - * and triggers and assertion if - * done so. + * @ref{Triangulation} class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. */ - bool get_face_orientation (const unsigned int face) const; + bool face_orientation (const unsigned int face) const; private: /** @@ -1616,20 +1632,28 @@ class TriaObjectAccessor<2, dim> : public TriaAccessor unsigned int number_of_children () const; /** - * Enquire some information about - * whether the face with the - * given number is in standard - * orientation or not. This - * information is only useful - * (and available) in 3d, see the + * Return whether the face with + * index @p{face} has its normal + * pointing in the standard + * direction (@p{true}) or + * whether it is the opposite + * (@p{false}). Which is the + * standard direction is + * documented with the + * @ref{Triangulation} class. In + * 1d and 2d, this is always + * @p{true}, but in 3d it may be + * different, see the respective + * discussion in the * documentation of the - * respective class. For all - * other dimensions, this - * function should not be called - * and triggers and assertion if - * done so. + * @ref{Triangulation} class. + * + * This function is really only + * for internal use in the + * library unless you absolutely + * know what this is all about. */ - bool get_face_orientation (const unsigned int face) const; + bool face_orientation (const unsigned int face) const; private: /** @@ -2098,7 +2122,7 @@ class TriaObjectAccessor<3, dim> : public TriaAccessor unsigned int number_of_children () const; /** - * Return whether the quad with + * Return whether the face with * index @p{face} has its normal * pointing in the standard * direction (@p{true}) or @@ -2106,6 +2130,12 @@ class TriaObjectAccessor<3, dim> : public TriaAccessor * (@p{false}). Which is the * standard direction is * documented with the + * @ref{Triangulation} class. In + * 1d and 2d, this is always + * @p{true}, but in 3d it may be + * different, see the respective + * discussion in the + * documentation of the * @ref{Triangulation} class. * * This function is really only @@ -2113,7 +2143,7 @@ class TriaObjectAccessor<3, dim> : public TriaAccessor * library unless you absolutely * know what this is all about. */ - bool get_face_orientation (const unsigned int face) const; + bool face_orientation (const unsigned int face) const; /** * Set whether the quad with diff --git a/deal.II/deal.II/include/grid/tria_accessor.templates.h b/deal.II/deal.II/include/grid/tria_accessor.templates.h index 98e5833172..7eb16fed5c 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.templates.h +++ b/deal.II/deal.II/include/grid/tria_accessor.templates.h @@ -25,6 +25,8 @@ #include +//TODO[WB]: TriaObjectAccessor should never be used since we have specializations. Can we remove the implementations of functions, or what are they there for? + /*------------------------ Functions: TriaAccessor ---------------------------*/ @@ -213,6 +215,16 @@ TriaObjectAccessor<1,dim>::max_refinement_depth () const +template +inline +bool +TriaObjectAccessor<1,dim>::face_orientation (const unsigned int) const +{ + return true; +} + + + template inline void @@ -400,6 +412,16 @@ TriaObjectAccessor<2,dim>::max_refinement_depth () const +template +inline +bool +TriaObjectAccessor<2,dim>::face_orientation (const unsigned int) const +{ + return true; +} + + + template inline void @@ -534,7 +556,7 @@ TriaObjectAccessor<3,dim>::line (const unsigned int i) const { 4, 3, 0 }}; return (this->quad(lookup_table[i][0]) - ->line(get_face_orientation(lookup_table[i][0]) ? + ->line(face_orientation(lookup_table[i][0]) ? lookup_table[i][1] : lookup_table[i][2])); } @@ -592,7 +614,7 @@ TriaObjectAccessor<3,dim>::line_index (const unsigned int i) const { 4, 3, 0 }}; return (this->quad(lookup_table[i][0]) - ->line_index(get_face_orientation(lookup_table[i][0]) ? + ->line_index(face_orientation(lookup_table[i][0]) ? lookup_table[i][1] : lookup_table[i][2])); } @@ -672,6 +694,29 @@ TriaObjectAccessor<3,dim>::max_refinement_depth () const +template +inline +bool +TriaObjectAccessor<3, dim>:: +face_orientation (const unsigned int face) const +{ + Assert (used(), typename TriaAccessor::ExcCellNotUsed()); + Assert (face::faces_per_cell, + ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell)); + Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face + < this->tria->levels[this->present_level] + ->hexes.face_orientations.size(), + ExcInternalError()); + + return (this->tria->levels[this->present_level] + ->hexes.face_orientations[this->present_index * + GeometryInfo<3>::faces_per_cell + + + face]); +} + + + template inline void diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index 92fb98e3b3..aaae80db02 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -62,7 +62,7 @@ namespace internal // face or subface case 3: return ((face_no + - (cell->get_face_orientation(face_no) == true ? + (cell->face_orientation(face_no) == true ? 0 : GeometryInfo::faces_per_cell)) * n_quadrature_points); @@ -90,6 +90,11 @@ namespace internal Assert (subface_no+1 < GeometryInfo::subfaces_per_face+1, ExcInternalError()); Assert (n_quadrature_points > 0, ExcInternalError()); + + Assert (cell->has_children() == false, + ExcMessage ("You can't use subface data for cells that are " + "already refined. Iterate over their children " + "instead in these cases.")); switch (dim) { @@ -103,10 +108,11 @@ namespace internal case 3: return (((face_no * GeometryInfo::subfaces_per_face + subface_no) - + (cell->get_face_orientation(face_no) == true ? + + (cell->face_orientation(face_no) == true ? 0 : GeometryInfo::faces_per_cell * - GeometryInfo::subfaces_per_face)) + GeometryInfo::subfaces_per_face) + ) * n_quadrature_points); default: Assert (false, ExcInternalError()); diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index d3663b6e9b..9afb26c52e 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -6124,12 +6124,12 @@ Triangulation<3>::execute_refinement () // opposite direction, so // store that up front const bool face_orientation[6] - = { hex->get_face_orientation (0), - hex->get_face_orientation (1), - hex->get_face_orientation (2), - hex->get_face_orientation (3), - hex->get_face_orientation (4), - hex->get_face_orientation (5) }; + = { hex->face_orientation (0), + hex->face_orientation (1), + hex->face_orientation (2), + hex->face_orientation (3), + hex->face_orientation (4), + hex->face_orientation (5) }; const unsigned int line_indices[30] = { @@ -6535,7 +6535,7 @@ Triangulation<3>::execute_refinement () for (unsigned int f=0; f::faces_per_cell; ++f) for (unsigned int s=0; s::subfaces_per_face; ++s) new_hexes[GeometryInfo::child_cell_on_face(f,s)] - ->set_face_orientation(f, hex->get_face_orientation(f)); + ->set_face_orientation(f, hex->face_orientation(f)); ///////////////////////////////// @@ -6684,9 +6684,9 @@ Triangulation<3>::execute_refinement () // all // wrong...) const bool orient - = (neighbor->get_face_orientation(nb_nb) + = (neighbor->face_orientation(nb_nb) && - hex->get_face_orientation(face)); + hex->face_orientation(face)); static const unsigned int child_switch_table[GeometryInfo::subfaces_per_face] diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc index fdd4a0e436..0b69a294e8 100644 --- a/deal.II/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/deal.II/source/grid/tria_accessor.cc @@ -250,15 +250,6 @@ unsigned int TriaObjectAccessor<1, dim>::number_of_children () const -template -bool -TriaObjectAccessor<1,dim>::get_face_orientation (const unsigned int) const -{ - Assert (false, ExcInternalError()); - return true; -} - - /*------------------------ Functions: QuadAccessor ---------------------------*/ template @@ -678,15 +669,6 @@ unsigned int TriaObjectAccessor<2, dim>::number_of_children () const -template -bool -TriaObjectAccessor<2,dim>::get_face_orientation (const unsigned int) const -{ - Assert (false, ExcInternalError()); - return true; -} - - /*------------------------ Functions: TriaObjectAccessor ---------------------------*/ @@ -715,14 +697,14 @@ int TriaObjectAccessor<3, dim>::vertex_index (const unsigned int corner) const { 0, 3, 2, 1 }; if (corner<4) { - if (get_face_orientation(0) == true) + if (face_orientation(0) == true) return quad(0)->vertex_index(corner); else return quad(0)->vertex_index(vertex_translation[corner]); } else { - if (get_face_orientation(1) == true) + if (face_orientation(1) == true) return quad(1)->vertex_index(corner-4); else return quad(1)->vertex_index(vertex_translation[corner-4]); @@ -1666,28 +1648,6 @@ unsigned int TriaObjectAccessor<3, dim>::number_of_children () const -template -bool -TriaObjectAccessor<3, dim>:: -get_face_orientation (const unsigned int face) const -{ - Assert (used(), typename TriaAccessor::ExcCellNotUsed()); - Assert (face::faces_per_cell, - ExcIndexRange (face, 0, GeometryInfo<3>::faces_per_cell)); - Assert (this->present_index * GeometryInfo<3>::faces_per_cell + face - < this->tria->levels[this->present_level] - ->hexes.face_orientations.size(), - ExcInternalError()); - - return (this->tria->levels[this->present_level] - ->hexes.face_orientations[this->present_index * - GeometryInfo<3>::faces_per_cell - + - face]); -} - - - template void TriaObjectAccessor<3, dim>:: diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index a9c4bf9884..0e0d61db9a 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -373,25 +373,36 @@ estimate_some (const Mapping &mapping, { const unsigned int n_solution_vectors = solutions.size(); - // make up a fe face values object for the - // restriction of the finite element function - // to a face, for the present cell and its - // neighbor. In principle we would only need - // one at a time, but this way we can - // have more fine grained access to what - // values really need to be computed (we - // need not compute all values on the - // neighbor cells, so using two objects + // make up a fe face values object + // for the restriction of the + // finite element function to a + // face, for the present cell and + // its neighbor. In principle we + // would only need one at a time, + // but this way we can have more + // fine grained access to what + // values really need to be + // computed (we need not compute + // all values on the neighbor + // cells, so using two objects // gives us a performance gain). + // + // in debug mode, make sure that + // some data matches, so compute + // quadrature points always FEFaceValues fe_face_values_cell (mapping, dof_handler.get_fe(), quadrature, - UpdateFlags(update_gradients | - update_JxW_values | - ((!neumann_bc.empty() || - (coefficients != 0)) ? - update_q_points : 0) | - update_normal_vectors)); + UpdateFlags( + update_gradients | + update_JxW_values | + ((!neumann_bc.empty() || + (coefficients != 0)) ? + update_q_points : 0) | +#ifdef DEBUG + update_q_points | +#endif + update_normal_vectors)); FEFaceValues fe_face_values_neighbor (mapping, dof_handler.get_fe(), quadrature, @@ -399,7 +410,12 @@ estimate_some (const Mapping &mapping, FESubfaceValues fe_subface_values (mapping, dof_handler.get_fe(), quadrature, - update_gradients); + UpdateFlags( + update_gradients +#ifdef DEBUG + | update_q_points +#endif + )); active_cell_iterator cell = dof_handler.begin_active(); @@ -987,43 +1003,57 @@ integrate_over_irregular_face (const DoFHandler &dof_handler, { // get an iterator pointing to the // cell behind the present subface - const bool orientation_flag - = (dim != 3 ? - true : - cell->get_face_orientation(face_no) == true - && - neighbor->get_face_orientation(neighbor_neighbor) == true); static const unsigned int subface_translation[4] = { 0, 3, 2, 1 }; + // see whether face and + // the neighbor's + // counterface share the + // same indexing of + // children. if not so, + // translate child + // indices + const bool face_orientations_match + = (neighbor->face_orientation(neighbor_neighbor) == + cell->face_orientation(face_no)); const unsigned int neighbor_child_index = (GeometryInfo:: child_cell_on_face(neighbor_neighbor, - (orientation_flag ? + (face_orientations_match ? subface_no : subface_translation[subface_no]))); const active_cell_iterator neighbor_child = neighbor->child(neighbor_child_index); - Assert (neighbor_child->face(neighbor_neighbor) == - cell->face(face_no)->child(subface_no), - ExcInternalError()); Assert (!neighbor->child(neighbor_child_index)->has_children(), ExcInternalError()); - // restrict the finite element on the - // present cell to the subface and - // store the gradient of the solution - // in psi + // restrict the finite element + // on the present cell to the + // subface fe_subface_values.reinit (cell, face_no, subface_no); + // restrict the finite element + // on the neighbor cell to the + // common @p{subface}. + fe_face_values.reinit (neighbor_child, neighbor_neighbor); + + // make sure that quadrature + // points match + for (unsigned int q=0; q &dof_handler, if (coefficients->n_components == 1) { coefficients->value_list (fe_face_values.get_quadrature_points(), - per_thread_data.coefficient_values1); + per_thread_data.coefficient_values1); for (unsigned int n=0; n