From 60c79f56f966937570fe708a9c402b024c3a5a75 Mon Sep 17 00:00:00 2001 From: wolf Date: Sat, 18 Oct 2003 19:20:41 +0000 Subject: [PATCH] Move computation of offsets in QProjections into the class that does the projections. Document this in various places. Several other updates, including purging the use of GeometryInfo::child_cell_on_face and replacing it by the use of cell->neighbor_child_on_subface. Update news.html. git-svn-id: https://svn.dealii.org/trunk@8108 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/quadrature.h | 180 ++++++++++++++++-- deal.II/base/source/quadrature.cc | 114 +++++++++++ deal.II/deal.II/include/fe/mapping_q1.h | 53 ++---- deal.II/deal.II/include/grid/geometry_info.h | 126 ++++++++++-- deal.II/deal.II/source/dofs/dof_tools.cc | 77 ++++---- deal.II/deal.II/source/fe/fe_dgp.cc | 19 +- deal.II/deal.II/source/fe/fe_dgq.cc | 19 +- deal.II/deal.II/source/fe/fe_nedelec.cc | 19 +- deal.II/deal.II/source/fe/fe_q.cc | 19 +- .../deal.II/source/fe/fe_q_hierarchical.cc | 19 +- .../deal.II/source/fe/fe_raviart_thomas.cc | 19 +- deal.II/deal.II/source/fe/fe_tools.cc | 1 - deal.II/deal.II/source/fe/fe_values.cc | 4 + .../deal.II/source/fe/mapping_cartesian.cc | 84 ++++---- deal.II/deal.II/source/fe/mapping_q.cc | 12 +- deal.II/deal.II/source/fe/mapping_q1.cc | 141 +------------- .../source/numerics/error_estimator.cc | 48 ++--- deal.II/doc/news/c-4-0.html | 46 ++++- deal.II/examples/step-12/step-12.cc | 18 +- deal.II/examples/step-14/step-14.cc | 4 +- 20 files changed, 630 insertions(+), 392 deletions(-) diff --git a/deal.II/base/include/base/quadrature.h b/deal.II/base/include/base/quadrature.h index 1a1d59946b..36a4bf1316 100644 --- a/deal.II/base/include/base/quadrature.h +++ b/deal.II/base/include/base/quadrature.h @@ -267,7 +267,7 @@ class QIterated : public Quadrature /** * This class is a helper class to facilitate the usage of quadrature formulae * on faces or subfaces of cells. It computes the locations of quadrature - * points on the unit cell from a quadrature object for a mannifold of + * points on the unit cell from a quadrature object for a manifold of * one dimension less than that of the cell and the number of the face. * For example, giving the Simpson rule in one dimension and using the * @p{project_to_face} function with face number 1, the returned points will @@ -280,19 +280,34 @@ class QIterated : public Quadrature * unit cell), the same applies as above. Note the order in which the * children of a face are numbered, which in two dimensions coincides * with the orientation of the face. + * + * The second set of functions generates a quadrature formula by + * projecting a given quadrature rule on @em{all} faces and + * subfaces. This is used in the @ref{FEFaceValues} and + * @ref{FESubfaceValues} classes. Since we now have the quadrature + * points of all faces and subfaces in one array, we need to have a + * way to find the starting index of the points and weights + * corresponding to one face or subface within this array. This is + * done through the @ref{DataSetDescriptor} member class. * - * The different functions are grouped into a common class to avoid putting - * them into global namespace (and to make documentation easier, since - * presently the documentation tool can only handle classes, not global - * functions). However, since they have no local data, all functions are - * declared @p{static} and can be called without creating an object of this - * class. + * The different functions are grouped into a common class to avoid + * putting them into global namespace. However, since they have no + * local data, all functions are declared @p{static} and can be + * called without creating an object of this class. + * + * For the 3d case, you should note that the orientation of faces is + * even more intricate than for two dimensions. Quadrature formulae + * are projected upon the faces in their standard orientation, not to + * the inside or outside of the hexahedron. Refer to the + * documentation of the @p{Triangulation} class for a description of + * the orientation of the different faces. To make things more + * complicated, in 3d we allow faces in two orientations (which can + * be identified using @p{cell->face_orientation(face)}), so we have + * to project quadrature formula onto faces and subfaces in two + * orientations. The @ref{DataSetDescriptor} member class is used to + * identify where each dataset starts. * - * For the 3d case, you should note that the orientation of faces is even - * more intricate than for two dimensions. Quadrature formulae are projected - * upon the faces in their standard orientation, not to the inside or outside - * of the hexahedron! Refer to the documentation of the @p{Triangulation} class - * for a description of the orientation of the different faces. + * @author Wolfgang Bangerth, 1998, 1999, 2003 */ template class QProjector @@ -399,6 +414,147 @@ class QProjector project_to_child (const Quadrature &quadrature, const unsigned int child_no); + /** + * Since the + * @p{project_to_all_faces} and + * @p{project_to_all_subfaces} + * functions chain together the + * quadrature points and weights + * of all projections of a face + * quadrature formula to the + * faces or subfaces of a cell, + * we need a way to identify + * where the starting index of + * the points and weights for a + * particular face or subface + * is. This class provides this: + * there are static member + * functions that generate + * objects of this type, given + * face or subface indices, and + * you can then use the generated + * object in place of an integer + * that denotes the offset of a + * given dataset. + * + * @author Wolfgang Bangerth, 2003 + */ + class DataSetDescriptor + { + public: + /** + * Default constructor. This + * doesn't do much except + * generating an invalid + * index, since you didn't + * give a valid descriptor of + * the cell, face, or subface + * you wanted. + */ + DataSetDescriptor (); + + /** + * Static function to + * generate the offset of a + * cell. Since we only have + * one cell per quadrature + * object, this offset is of + * course zero, but we carry + * this function around for + * consistency with the other + * static functions. + */ + static DataSetDescriptor cell (); + + /** + * Static function to + * generate an offset object + * for a given face of a cell + * with the given face + * orientation. This function + * of course is only allowed + * if @p{dim>=2}, and the + * face orientation is + * ignored if the space + * dimension equals 2. + * + * The last argument denotes + * the number of quadrature + * points the + * lower-dimensional face + * quadrature formula (the + * one that has been + * projected onto the faces) + * has. + */ + static + DataSetDescriptor + face (const unsigned int face_no, + const bool face_orientation, + const unsigned int n_quadrature_points); + + /** + * Static function to + * generate an offset object + * for a given subface of a + * cell with the given face + * orientation. This function + * of course is only allowed + * if @p{dim>=2}, and the + * face orientation is + * ignored if the space + * dimension equals 2. + * + * The last argument denotes + * the number of quadrature + * points the + * lower-dimensional face + * quadrature formula (the + * one that has been + * projected onto the faces) + * has. + */ + static + DataSetDescriptor + sub_face (const unsigned int face_no, + const unsigned int subface_no, + const bool face_orientation, + const unsigned int n_quadrature_points); + + /** + * Conversion operator to an + * integer denoting the + * offset of the first + * element of this dataset in + * the set of quadrature + * formulas all projected + * onto faces and + * subfaces. This conversion + * operator allows us to use + * offset descriptor objects + * in place of integer + * offsets. + */ + operator unsigned int () const; + + private: + /** + * Store the integer offset + * for a given cell, face, or + * subface. + */ + const unsigned int dataset_offset; + + /** + * This is the real + * constructor, but it is + * private and thus only + * available to the static + * member functions above. + */ + DataSetDescriptor (const unsigned int dataset_offset); + }; + private: /** * Given a quadrature object in diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc index 182c84d4ac..b591d3c5c5 100644 --- a/deal.II/base/source/quadrature.cc +++ b/deal.II/base/source/quadrature.cc @@ -878,6 +878,120 @@ QProjector::project_to_child (const Quadrature &quadrature, +template +typename QProjector::DataSetDescriptor +QProjector::DataSetDescriptor::cell () +{ + return 0; +} + + +template +typename QProjector::DataSetDescriptor +QProjector::DataSetDescriptor:: +face (const unsigned int face_no, + const bool face_orientation, + const unsigned int n_quadrature_points) +{ + Assert (dim != 1, ExcInternalError()); + Assert (face_no < GeometryInfo::faces_per_cell, + ExcInternalError()); + Assert (n_quadrature_points > 0, ExcInternalError()); + + switch (dim) + { + case 1: + case 2: + return face_no * n_quadrature_points; + + // in 3d, we have to + // account for faces + // that have reverse + // orientation. thus, + // we have to store + // _two_ data sets per + // face or subface + case 3: + return ((face_no + + (face_orientation == true ? + 0 : GeometryInfo::faces_per_cell)) + * n_quadrature_points); + + default: + Assert (false, ExcInternalError()); + } + return static_cast(-1); +} + + + +template +typename QProjector::DataSetDescriptor +QProjector::DataSetDescriptor:: +sub_face (const unsigned int face_no, + const unsigned int subface_no, + const bool face_orientation, + const unsigned int n_quadrature_points) +{ + Assert (dim != 1, ExcInternalError()); + Assert (face_no < GeometryInfo::faces_per_cell, + ExcInternalError()); + // the trick with +1 prevents + // that we get a warning in 1d + Assert (subface_no+1 < GeometryInfo::subfaces_per_face+1, + ExcInternalError()); + Assert (n_quadrature_points > 0, ExcInternalError()); + + switch (dim) + { + case 1: + case 2: + return ((face_no * GeometryInfo::subfaces_per_face + + subface_no) + * n_quadrature_points); + + // for 3d, same as above: + case 3: + return (((face_no * GeometryInfo::subfaces_per_face + + subface_no) + + (face_orientation == true ? + 0 : + GeometryInfo::faces_per_cell * + GeometryInfo::subfaces_per_face) + ) + * n_quadrature_points); + default: + Assert (false, ExcInternalError()); + } + return static_cast(-1); +} + + +template +QProjector::DataSetDescriptor::operator unsigned int () const +{ + return dataset_offset; +} + + + +template +QProjector::DataSetDescriptor:: +DataSetDescriptor (const unsigned int dataset_offset) + : + dataset_offset (dataset_offset) +{} + + +template +QProjector::DataSetDescriptor:: +DataSetDescriptor () + : + dataset_offset (static_cast(-1)) +{} + + + // ------------------------------------------------------------ // diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index fc7080e53a..f3dce621a5 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -16,49 +16,21 @@ #include #include -#include +#include #include #include #include +#include -/*!@addtogroup fe */ -/*@{*/ -namespace internal -{ - template - class DataSetDescriptor - { - public: - DataSetDescriptor (); - - static DataSetDescriptor cell (); - - static - DataSetDescriptor - face (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int n_quadrature_points); - - static - DataSetDescriptor - sub_face (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const unsigned int n_quadrature_points); - - unsigned int offset () const; - - private: - const unsigned int dataset_offset; - - DataSetDescriptor (const unsigned int dataset_offset); - }; -} +/*!@addtogroup fe */ +/*@{*/ + + /** * Mapping of general quadrilateral/hexahedra by d-linear shape * functions. @@ -385,12 +357,17 @@ class MappingQ1 : public Mapping DeclException0 (ExcAccessToUninitializedField); protected: + /** - * Typedef the right data set - * descriptor to a local type to - * make use somewhat simpler. + * Declare a convenience typedef + * for the class that describes + * offsets into quadrature + * formulas projected onto faces + * and subfaces. */ - typedef internal::DataSetDescriptor DataSetDescriptor; + typedef + typename QProjector::DataSetDescriptor + DataSetDescriptor; /** * Implementation of the interface in diff --git a/deal.II/deal.II/include/grid/geometry_info.h b/deal.II/deal.II/include/grid/geometry_info.h index ba01b98725..fe10fd3ca8 100644 --- a/deal.II/deal.II/include/grid/geometry_info.h +++ b/deal.II/deal.II/include/grid/geometry_info.h @@ -159,9 +159,34 @@ struct GeometryInfo * 3 and 2, etc. * * For three spatial dimensions, - * the exact order of the children is - * laid down in the documentation of - * the @ref{Triangulation} class. + * the exact order of the + * children is laid down in the + * documentation of the + * @ref{Triangulation} class. + * However, it must be noted that + * this class and function only + * deals with faces in standard + * orientation. In 3d, faces can + * exist in two orientations, + * though, and if a face is in + * the wrong orientation, then + * this function may not give you + * what you want. You can inquire + * about the face orientation + * using the + * @p{cell->face_orientation} + * function, and the function to + * ask for a neighbor's cell + * behind a given face and + * subface is + * @p{cell->neighbor_child_on_subface}. + * The latter function, in + * contrast to the present one, + * also takes into account the + * actual orientation of the + * faces of a cell and will + * return the correct result in + * all cases. */ static unsigned int child_cell_on_face (const unsigned int face, const unsigned int subface); @@ -510,9 +535,35 @@ struct GeometryInfo<1> * 3 and 2, etc. * * For three spatial dimensions, - * the exact order of the children is - * laid down in the documentation of - * the @ref{Triangulation} class. + * the exact order of the + * children is laid down in the + * documentation of the + * @ref{Triangulation} + * class. However, it must be + * noted that this class and + * function only deals with faces + * in standard orientation. In + * 3d, faces can exist in two + * orientations, though, and if a + * face is in the wrong + * orientation, then this + * function may not give you what + * you want. You can inquire + * about the face orientation + * using the + * @p{cell->face_orientation} + * function, and the function to + * ask for a neighbor's cell + * behind a given face and + * subface is + * @p{cell->neighbor_child_on_subface}. + * The latter function, in + * contrast to the present one, + * also takes into account the + * actual orientation of the + * faces of a cell and will + * return the correct result in + * all cases. */ static unsigned int child_cell_on_face (const unsigned int face, const unsigned int subface); @@ -769,9 +820,35 @@ struct GeometryInfo<2> * 3 and 2, etc. * * For three spatial dimensions, - * the exact order of the children is - * laid down in the documentation of - * the @ref{Triangulation} class. + * the exact order of the + * children is laid down in the + * documentation of the + * @ref{Triangulation} + * class. However, it must be + * noted that this class and + * function only deals with faces + * in standard orientation. In + * 3d, faces can exist in two + * orientations, though, and if a + * face is in the wrong + * orientation, then this + * function may not give you what + * you want. You can inquire + * about the face orientation + * using the + * @p{cell->face_orientation} + * function, and the function to + * ask for a neighbor's cell + * behind a given face and + * subface is + * @p{cell->neighbor_child_on_subface}. + * The latter function, in + * contrast to the present one, + * also takes into account the + * actual orientation of the + * faces of a cell and will + * return the correct result in + * all cases. */ static unsigned int child_cell_on_face (const unsigned int face, const unsigned int subface); @@ -1030,9 +1107,34 @@ struct GeometryInfo<3> * 3 and 2, etc. * * For three spatial dimensions, - * the exact order of the children is - * laid down in the documentation of - * the @ref{Triangulation} class. + * the exact order of the + * children is laid down in the + * documentation of the + * @ref{Triangulation} + * class. However, it must be + * noted that this class and + * function only deals with faces + * in standard orientation. In + * 3d, faces can exist in two + * orientations, though, and if a + * face is in the wrong + * orientation, then this + * function may not give you what + * you want. You can inquire + * about the face orientation + * using the + * @p{cell->face_orientation} + * function, and the function to + * ask for a neighbor's cell + * behind a given face and + * subface is + * @p{cell->neighbor_child_on_subface}. + * The latter function, in contrast + * to the present one, also takes + * into account the actual + * orientation of the faces of a + * cell and will return the + * correct result in all cases. */ static unsigned int child_cell_on_face (const unsigned int face, const unsigned int subface); diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index de6186001e..aaa7d1aaea 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -373,7 +373,8 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler &dof, if (neighbor->level() < cell->level()) continue; - unsigned int neighbor_face = cell->neighbor_of_neighbor(face); + const unsigned int neighbor_face + = cell->neighbor_of_neighbor(face); if (neighbor->has_children()) { @@ -381,35 +382,35 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler &dof, sub_nr != GeometryInfo::subfaces_per_face; ++sub_nr) { - typename DoFHandler::cell_iterator sub_neighbor - = neighbor-> - child(GeometryInfo::child_cell_on_face(neighbor_face, sub_nr)); + const typename DoFHandler::cell_iterator + sub_neighbor + = cell->neighbor_child_on_subface (face, sub_nr); sub_neighbor->get_dof_indices (dofs_on_other_cell); for (unsigned int i=0; iface(neighbor_face)->set_user_flag (); } - } else { + } + else + { neighbor->get_dof_indices (dofs_on_other_cell); for (unsigned int i=0; iface(neighbor_face)->set_user_flag (); } } @@ -575,7 +576,8 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler &dof, face < GeometryInfo::faces_per_cell; ++face) { - typename DoFHandler::face_iterator cell_face = cell->face(face); + const typename DoFHandler::face_iterator + cell_face = cell->face(face); if (cell_face->user_flag_set ()) continue; @@ -589,16 +591,12 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler &dof, const bool j_non_zero_i = support_on_face (j, face); if (flux_dof_mask(i,j) == 'f') - { - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); - } + sparsity.add (dofs_on_this_cell[i], + dofs_on_this_cell[j]); if (flux_dof_mask(i,j) == 'a') - { - if (i_non_zero_i && j_non_zero_i) - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); - } + if (i_non_zero_i && j_non_zero_i) + sparsity.add (dofs_on_this_cell[i], + dofs_on_this_cell[j]); } } } @@ -611,7 +609,8 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler &dof, if (neighbor->level() < cell->level()) continue; - unsigned int neighbor_face = cell->neighbor_of_neighbor(face); + const unsigned int + neighbor_face = cell->neighbor_of_neighbor(face); if (neighbor->has_children()) { @@ -619,9 +618,9 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler &dof, sub_nr != GeometryInfo::subfaces_per_face; ++sub_nr) { - typename DoFHandler::cell_iterator sub_neighbor - = neighbor-> - child(GeometryInfo::child_cell_on_face(neighbor_face, sub_nr)); + const typename DoFHandler::cell_iterator + sub_neighbor + = cell->neighbor_child_on_subface (face, sub_nr); sub_neighbor->get_dof_indices (dofs_on_other_cell); for (unsigned int i=0; i &dof, } sub_neighbor->face(neighbor_face)->set_user_flag (); } - } else { + } + else + { neighbor->get_dof_indices (dofs_on_other_cell); for (unsigned int i=0; i #include #include -#include #include #ifdef HAVE_STD_STRINGSTREAM @@ -352,7 +351,7 @@ FE_DGP::fill_fe_values (const Mapping &mapping, if (flags & update_second_derivatives) this->compute_2nd (mapping, cell, - internal::DataSetDescriptor::cell().offset(), + QProjector::DataSetDescriptor::cell(), mapping_data, fe_data, data); } @@ -377,10 +376,10 @@ FE_DGP::fill_fe_face_values (const Mapping &mapping, // offset determines which data set // to take (all data sets for all // faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - face (cell, face, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + face (face, cell->face_orientation(face), + quadrature.n_quadrature_points)); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); @@ -428,10 +427,10 @@ FE_DGP::fill_fe_subface_values (const Mapping &mappi // offset determines which data set // to take (all data sets for all // sub-faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - sub_face (cell, face, subface, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + sub_face (face, subface, cell->face_orientation(face), + quadrature.n_quadrature_points)); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index 3accfb98c4..79ace133ba 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -18,7 +18,6 @@ #include #include #include -#include #include #include #include @@ -766,7 +765,7 @@ FE_DGQ::fill_fe_values (const Mapping &mapping, if (flags & update_second_derivatives) this->compute_2nd (mapping, cell, - internal::DataSetDescriptor::cell().offset(), + QProjector::DataSetDescriptor::cell(), mapping_data, fe_data, data); } @@ -791,10 +790,10 @@ FE_DGQ::fill_fe_face_values (const Mapping &mapping, // offset determines which data set // to take (all data sets for all // faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - face (cell, face, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + face (face, cell->face_orientation(face), + quadrature.n_quadrature_points)); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); @@ -842,10 +841,10 @@ FE_DGQ::fill_fe_subface_values (const Mapping &mappi // offset determines which data set // to take (all data sets for all // sub-faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - sub_face (cell, face, subface, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + sub_face (face, subface, cell->face_orientation(face), + quadrature.n_quadrature_points)); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index 03a21ead75..de12d5955c 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -18,7 +18,6 @@ #include #include #include -#include #include #include @@ -1098,7 +1097,7 @@ FE_Nedelec::fill_fe_values (const Mapping &mapping, if (flags & update_second_derivatives) this->compute_2nd (mapping, cell, - internal::DataSetDescriptor::cell().offset(), + QProjector::DataSetDescriptor::cell(), mapping_data, fe_data, data); } @@ -1123,10 +1122,10 @@ FE_Nedelec::fill_fe_face_values (const Mapping &mapp // offset determines which data set // to take (all data sets for all // faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - face (cell, face, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + face (face, cell->face_orientation(face), + quadrature.n_quadrature_points)); // get the flags indicating the // fields that have to be filled @@ -1273,10 +1272,10 @@ FE_Nedelec::fill_fe_subface_values (const Mapping &m // offset determines which data set // to take (all data sets for all // faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - sub_face (cell, face, subface, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + sub_face (face, subface, cell->face_orientation(face), + quadrature.n_quadrature_points)); // get the flags indicating the // fields that have to be filled diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 572b02e9a4..f0a3aaa2ea 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -19,7 +19,6 @@ #include #include #include -#include #include #include @@ -1689,7 +1688,7 @@ FE_Q::fill_fe_values (const Mapping &mapping, if (flags & update_second_derivatives) this->compute_2nd (mapping, cell, - internal::DataSetDescriptor::cell().offset(), + QProjector::DataSetDescriptor::cell(), mapping_data, fe_data, data); } @@ -1714,10 +1713,10 @@ FE_Q::fill_fe_face_values (const Mapping &mapping, // offset determines which data set // to take (all data sets for all // faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - face (cell, face, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + face (face, cell->face_orientation(face), + quadrature.n_quadrature_points)); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); @@ -1765,10 +1764,10 @@ FE_Q::fill_fe_subface_values (const Mapping &mapping // offset determines which data set // to take (all data sets for all // sub-faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - sub_face (cell, face, subface, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + sub_face (face, subface, cell->face_orientation(face), + quadrature.n_quadrature_points)); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); diff --git a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc index 39cb5bd1e9..5a4ff17599 100644 --- a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc +++ b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc @@ -19,7 +19,6 @@ #include #include #include -#include #include #include @@ -1061,7 +1060,7 @@ FE_Q_Hierarchical::fill_fe_values ( if (flags & update_second_derivatives) this->compute_2nd (mapping, cell, - internal::DataSetDescriptor::cell().offset(), + QProjector::DataSetDescriptor::cell(), mapping_data, fe_data, data); } @@ -1087,10 +1086,10 @@ FE_Q_Hierarchical::fill_fe_face_values ( // offset determines which data set // to take (all data sets for all // faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - face (cell, face, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + face (face, cell->face_orientation(face), + quadrature.n_quadrature_points)); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); @@ -1139,10 +1138,10 @@ FE_Q_Hierarchical::fill_fe_subface_values ( // offset determines which data set // to take (all data sets for all // sub-faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - sub_face (cell, face, subface, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + sub_face (face, subface, cell->face_orientation(face), + quadrature.n_quadrature_points)); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc index 28d837df14..b63f5b854b 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -18,7 +18,6 @@ #include #include #include -#include #include #include @@ -1611,7 +1610,7 @@ FE_RaviartThomas::fill_fe_values (const Mapping &map if (flags & update_second_derivatives) this->compute_2nd (mapping, cell, - internal::DataSetDescriptor::cell().offset(), + QProjector::DataSetDescriptor::cell(), mapping_data, fe_data, data); } @@ -1636,10 +1635,10 @@ FE_RaviartThomas::fill_fe_face_values (const Mapping // offset determines which data set // to take (all data sets for all // faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - face (cell, face, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + face (face, cell->face_orientation(face), + quadrature.n_quadrature_points)); // get the flags indicating the // fields that have to be filled @@ -1778,10 +1777,10 @@ FE_RaviartThomas::fill_fe_subface_values (const Mapping // offset determines which data set // to take (all data sets for all // faces are stored contiguously) - const unsigned int offset - = (internal::DataSetDescriptor:: - sub_face (cell, face, subface, - quadrature.n_quadrature_points)).offset(); + const typename QProjector::DataSetDescriptor offset + = (QProjector::DataSetDescriptor:: + sub_face (face, subface, cell->face_orientation(face), + quadrature.n_quadrature_points)); // get the flags indicating the // fields that have to be filled diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 461ffb1fde..7b9efaeb96 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -29,7 +29,6 @@ #include #include #include -#include #include #include #include diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 0ea0cdda92..6ae722cd54 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -865,6 +865,10 @@ void FESubfaceValues::reinit (const typename DoFHandler::cell_iterator ExcIndexRange (face_no, 0, GeometryInfo::faces_per_cell)); Assert (subface_no < GeometryInfo::subfaces_per_face, ExcIndexRange (subface_no, 0, GeometryInfo::subfaces_per_face)); + 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.")); this->present_cell = cell; this->my_orientation = this->orientation_table[face_no]; diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index c7ebc88290..5bfcb69184 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -19,7 +19,6 @@ #include #include #include -#include #include #include @@ -135,13 +134,29 @@ MappingCartesian::compute_fill (const typename DoFHandler::cell_iterat std::vector > &quadrature_points, std::vector > &normal_vectors) const { - UpdateFlags update_flags(data.current_update_flags()); + const UpdateFlags update_flags(data.current_update_flags()); const unsigned int npts = quadrature_points.size (); - typedef typename internal::DataSetDescriptor DataSetDescriptor; - unsigned int offset = DataSetDescriptor::cell().offset(); + const typename QProjector::DataSetDescriptor offset + = (face_no == invalid_face_number + ? + QProjector::DataSetDescriptor::cell() + : + (sub_no == invalid_face_number + ? + // called from FEFaceValues + QProjector::DataSetDescriptor::face (face_no, + cell->face_orientation(face_no), + quadrature_points.size()) + : + // called from FESubfaceValues + QProjector::DataSetDescriptor::sub_face (face_no, sub_no, + cell->face_orientation(face_no), + quadrature_points.size()) + )); + // some more sanity checks if (face_no != invalid_face_number) { // Add 1 on both sides of @@ -151,23 +166,11 @@ MappingCartesian::compute_fill (const typename DoFHandler::cell_iterat Assert (face_no+1 < GeometryInfo::faces_per_cell+1, ExcIndexRange (face_no, 0, GeometryInfo::faces_per_cell)); - if (sub_no == invalid_face_number) - // called from FEFaceValues - offset = DataSetDescriptor::face (cell, face_no, - quadrature_points.size()).offset(); - else - { - // called from - // FESubfaceValue (do the - // +1 trick to avoid a - // warning about comparison - // always being false in - // 1d) - Assert (sub_no+1 < GeometryInfo::subfaces_per_face+1, - ExcIndexRange (sub_no, 0, GeometryInfo::subfaces_per_face)); - offset = DataSetDescriptor::sub_face (cell, face_no, sub_no, - quadrature_points.size()).offset(); - } + Assert ((sub_no == invalid_face_number) + || + (sub_no+1 < GeometryInfo::subfaces_per_face+1), + ExcIndexRange (sub_no, 0, + GeometryInfo::subfaces_per_face)); } else // invalid face number, so @@ -365,19 +368,16 @@ MappingCartesian::fill_fe_face_values (const typename DoFHandler::cell J *= data.length[d]; if (data.current_update_flags() & update_JxW_values) - { - for (unsigned int i=0; i void MappingCartesian::fill_fe_subface_values (const typename DoFHandler::cell_iterator &cell, @@ -411,16 +411,12 @@ MappingCartesian::fill_fe_subface_values (const typename DoFHandler::c J *= data.length[d]; if (data.current_update_flags() & update_JxW_values) - { - for (unsigned int i=0; i::subfaces_per_face; - } + for (unsigned int i=0; i::subfaces_per_face; if (data.current_update_flags() & update_boundary_forms) - { - for (unsigned int i=0; i::transform_covariant (Tensor<1,dim> *begin, { for (unsigned int d=0;d::transform_covariant (Tensor<2,dim> *begin, for (unsigned int d=0; d::transform_contravariant (Tensor<2,dim> *begin, for (unsigned int d=0; d::transform_unit_to_real_cell ( Point p_real = cell->vertex(0); for (unsigned int d=0; d #include #include -#include #include #include @@ -339,8 +338,9 @@ MappingQ::fill_fe_face_values (const typename DoFHandler::cell_iterato const unsigned int n_q_points=q.n_quadrature_points; this->compute_fill_face (cell, face_no, false, n_q_points, - MappingQ1::DataSetDescriptor:: - face (cell, face_no, n_q_points), + QProjector::DataSetDescriptor:: + face (face_no, cell->face_orientation(face_no), + n_q_points), q.get_weights(), *p_data, quadrature_points, JxW_values, @@ -395,8 +395,10 @@ MappingQ::fill_fe_subface_values (const typename DoFHandler::cell_iter const unsigned int n_q_points=q.n_quadrature_points; this->compute_fill_face (cell, face_no, true, n_q_points, - MappingQ1::DataSetDescriptor:: - sub_face (cell, face_no, sub_no, n_q_points), + QProjector::DataSetDescriptor:: + sub_face (face_no, sub_no, + cell->face_orientation(face_no), + n_q_points), q.get_weights(), *p_data, quadrature_points, JxW_values, diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index aaae80db02..0c322fe21d 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -18,136 +18,13 @@ #include #include #include -#include #include +#include #include #include -namespace internal -{ - template - DataSetDescriptor - DataSetDescriptor::cell () - { - return 0; - } - - - template - DataSetDescriptor - DataSetDescriptor:: - face (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int n_quadrature_points) - { - Assert (dim != 1, ExcInternalError()); - Assert (face_no < GeometryInfo::faces_per_cell, - ExcInternalError()); - Assert (n_quadrature_points > 0, ExcInternalError()); - - switch (dim) - { - case 1: - case 2: - return face_no * n_quadrature_points; - - // in 3d, we have to - // account for faces - // that have reverse - // orientation. thus, - // we have to store - // _two_ data sets per - // face or subface - case 3: - return ((face_no + - (cell->face_orientation(face_no) == true ? - 0 : GeometryInfo::faces_per_cell)) - * n_quadrature_points); - - default: - Assert (false, ExcInternalError()); - } - return static_cast(-1); - } - - - - template - DataSetDescriptor - DataSetDescriptor:: - sub_face (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const unsigned int n_quadrature_points) - { - Assert (dim != 1, ExcInternalError()); - Assert (face_no < GeometryInfo::faces_per_cell, - ExcInternalError()); - // the trick with +1 prevents - // that we get a warning in 1d - 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) - { - case 1: - case 2: - return ((face_no * GeometryInfo::subfaces_per_face + - subface_no) - * n_quadrature_points); - - // for 3d, same as above: - case 3: - return (((face_no * GeometryInfo::subfaces_per_face + - subface_no) - + (cell->face_orientation(face_no) == true ? - 0 : - GeometryInfo::faces_per_cell * - GeometryInfo::subfaces_per_face) - ) - * n_quadrature_points); - default: - Assert (false, ExcInternalError()); - } - return static_cast(-1); - } - - - template - unsigned int - DataSetDescriptor::offset () const - { - return dataset_offset; - } - - - - template - DataSetDescriptor:: - DataSetDescriptor (const unsigned int dataset_offset) - : - dataset_offset (dataset_offset) - {} - - - template - DataSetDescriptor:: - DataSetDescriptor () - : - dataset_offset (static_cast(-1)) - {} -} - - - template @@ -625,7 +502,7 @@ MappingQ1::compute_fill (const typename DoFHandler::cell_iterator &cel for (unsigned int point=0; point::compute_fill (const typename DoFHandler::cell_iterator &cel for (unsigned int i=0; i::fill_fe_face_values (const typename DoFHandler::cell_iterat compute_fill_face (cell, face_no, false, n_q_points, - DataSetDescriptor::face (cell, face_no, n_q_points), + DataSetDescriptor::face (face_no, + cell->face_orientation(face_no), + n_q_points), q.get_weights(), data, quadrature_points, @@ -841,7 +720,9 @@ MappingQ1::fill_fe_subface_values (const typename DoFHandler::cell_ite compute_fill_face (cell, face_no, true, n_q_points, - DataSetDescriptor::sub_face (cell, face_no, sub_no, n_q_points), + DataSetDescriptor::sub_face (face_no, sub_no, + cell->face_orientation(face_no), + n_q_points), q.get_weights(), data, quadrature_points, @@ -1163,7 +1044,3 @@ MappingQ1::transform_real_to_unit_cell_internal ( //----------------------------------------------------------------------// template class MappingQ1; -namespace internal -{ - template class DataSetDescriptor; -} diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index 907dc6edb8..aee011d7c9 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -802,11 +802,13 @@ integrate_over_regular_face (const DoFHandler &dof_handler, const active_cell_iterator neighbor = cell->neighbor(face_no); - // find which number the current - // face has relative to the neighboring - // cell - const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no); - Assert (neighbor_neighbor::faces_per_cell, ExcInternalError()); + // find which number the + // current face has relative to + // the neighboring cell + const unsigned int neighbor_neighbor + = cell->neighbor_of_neighbor (face_no); + Assert (neighbor_neighbor::faces_per_cell, + ExcInternalError()); // get restriction of finite element // function of @p{neighbor} to the @@ -984,38 +986,24 @@ integrate_over_irregular_face (const DoFHandler &dof_handler, // index the number of the // quadrature point - // store which number @p{cell} has in the - // list of neighbors of @p{neighbor} - const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no); - Assert (neighbor_neighbor::faces_per_cell, ExcInternalError()); + // store which number @p{cell} has + // in the list of neighbors of + // @p{neighbor} + const unsigned int neighbor_neighbor + = cell->neighbor_of_neighbor (face_no); + Assert (neighbor_neighbor::faces_per_cell, + ExcInternalError()); // loop over all subfaces - for (unsigned int subface_no=0; subface_no::subfaces_per_face; + for (unsigned int subface_no=0; + subface_no::subfaces_per_face; ++subface_no) { // get an iterator pointing to the // cell behind the present subface - 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, - (face_orientations_match ? - subface_no : - subface_translation[subface_no]))); const active_cell_iterator neighbor_child - = neighbor->child(neighbor_child_index); - Assert (!neighbor->child(neighbor_child_index)->has_children(), + = cell->neighbor_child_on_subface (face_no, subface_no); + Assert (!neighbor_child->has_children(), ExcInternalError()); // restrict the finite element diff --git a/deal.II/doc/news/c-4-0.html b/deal.II/doc/news/c-4-0.html index 27c1df1674..9699185020 100644 --- a/deal.II/doc/news/c-4-0.html +++ b/deal.II/doc/news/c-4-0.html @@ -32,15 +32,26 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

Incompatibilities

    +
  1. + Changed: The QProjector has functions that + project a lower-dimensional quadrature formula onto all faces or + subfaces of a cell. In 3d, it now does this but also adds projections of + these quadrature formula onto the faces from the other side. We need + this due to the fact that we now support faces in 3d that have a normal + vector opposite to the standard direction. +
    + (WB 2003/10/17) +

    +
  2. Moved and changed: The header file - include/numerics/dof_renumbering.h has been moved to the - directory include/dofs, where it logically - belongs. Furthermore, the sorting parameter of the function DoFRenumbering::component_wise has changed its meaning. See - the reference documentation for details. -
    - (GK 2003/07/03) + include/numerics/dof_renumbering.h has been moved to the + directory include/dofs, where it logically + belongs. Furthermore, the sorting parameter of the function DoFRenumbering::component_wise has changed its meaning. See + the reference documentation for details. +
    + (GK 2003/07/03)

@@ -48,6 +59,16 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

General

    +
  1. + Augmented: The GeometryInfo::child_cell_on_face + results in a result that might not be what you expect in 3d in some + cases. A warning has been added to the documentation, as well as a + reference to the function to call instead. +
    + (WB 2003/10/17) +

    +
  2. Fixed: The step-14 program had a bug in the rare case that there are more CPUs in a machine than there are cells. This is now fixed. @@ -235,6 +256,15 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

    deal.II

      +
    1. + New: There is now a function CellAccessor::neighbor_child_on_subface + that should be called instead of using GeometryInfo::child_cell_on_face in most cases. +
      + (WB 2003/10/17) +

      +
    2. New: GridGenerator has a new function subdivided_hyper_rectangle diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc index f6e8c1916a..20c6f5fe0e 100644 --- a/deal.II/examples/step-12/step-12.cc +++ b/deal.II/examples/step-12/step-12.cc @@ -959,9 +959,9 @@ void DGMethod::assemble_system1 () // `behind' the // current // subface. - typename DoFHandler::active_cell_iterator neighbor_child= - neighbor->child(GeometryInfo:: - child_cell_on_face(neighbor2,subface_no)); + typename DoFHandler::active_cell_iterator + neighbor_child + = cell->neighbor_child_on_subface (face_no, subface_no); // As these are // quite @@ -1123,9 +1123,10 @@ void DGMethod::assemble_system1 () const unsigned int neighbor_face_no=faceno_subfaceno.first, neighbor_subface_no=faceno_subfaceno.second; - Assert (neighbor->neighbor(neighbor_face_no) - ->child(GeometryInfo::child_cell_on_face( - face_no,neighbor_subface_no)) == cell, ExcInternalError()); + Assert (neighbor->neighbor_child_on_subface (neighbor_face_no, + neighbor_subface_no) + == cell, + ExcInternalError()); // Reinit the // appropriate @@ -1310,9 +1311,8 @@ void DGMethod::assemble_system2 () subface_no::subfaces_per_face; ++subface_no) { - typename DoFHandler::cell_iterator neighbor_child= - neighbor->child(GeometryInfo::child_cell_on_face( - neighbor2,subface_no)); + typename DoFHandler::cell_iterator neighbor_child + = cell->neighbor_child_on_subface (face_no, subface_no); Assert (neighbor_child->face(neighbor2) == face->child(subface_no), ExcInternalError()); Assert (!neighbor_child->has_children(), ExcInternalError()); diff --git a/deal.II/examples/step-14/step-14.cc b/deal.II/examples/step-14/step-14.cc index 2ca8e13af3..11a158a714 100644 --- a/deal.II/examples/step-14/step-14.cc +++ b/deal.II/examples/step-14/step-14.cc @@ -3557,9 +3557,7 @@ namespace LaplaceSolver // assertion will be removed // anyway. const active_cell_iterator neighbor_child - = neighbor->child(GeometryInfo:: - child_cell_on_face(neighbor_neighbor, - subface_no)); + = cell->neighbor_child_on_subface (face_no, subface_no); Assert (neighbor_child->face(neighbor_neighbor) == cell->face(face_no)->child(subface_no), ExcInternalError()); -- 2.39.5