From 363426e8d93613a836f2bc4936f8bd338d598eab Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 30 Sep 2003 14:48:29 +0000 Subject: [PATCH] Overhaul the way we compute offsets. This was a big mess and unintelligible. Now have a descriptor for the data set we want to access. git-svn-id: https://svn.dealii.org/trunk@8072 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/mapping_q1.h | 50 ++++- .../deal.II/source/fe/mapping_cartesian.cc | 19 +- deal.II/deal.II/source/fe/mapping_q.cc | 23 +- deal.II/deal.II/source/fe/mapping_q1.cc | 199 +++++++++++++----- 4 files changed, 218 insertions(+), 73 deletions(-) diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index c6e24aa129..fc7080e53a 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -25,6 +25,40 @@ /*!@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); + }; +} + + + /** * Mapping of general quadrilateral/hexahedra by d-linear shape * functions. @@ -351,6 +385,12 @@ class MappingQ1 : public Mapping DeclException0 (ExcAccessToUninitializedField); protected: + /** + * Typedef the right data set + * descriptor to a local type to + * make use somewhat simpler. + */ + typedef internal::DataSetDescriptor DataSetDescriptor; /** * Implementation of the interface in @@ -442,9 +482,9 @@ class MappingQ1 : public Mapping * @p{fill_*} functions. */ void compute_fill (const typename DoFHandler::cell_iterator &cell, - const unsigned int npts, - const unsigned int offset, - InternalData &data, + const unsigned int npts, + const DataSetDescriptor data_set, + InternalData &data, std::vector > &quadrature_points) const; /** @@ -455,7 +495,7 @@ class MappingQ1 : public Mapping const unsigned int face_no, const bool is_subface, const unsigned int npts, - const unsigned int offset, + const DataSetDescriptor data_set, const std::vector &weights, InternalData &mapping_data, std::vector > &quadrature_points, @@ -683,7 +723,7 @@ template <> void MappingQ1<1>::compute_fill_face ( const unsigned int, const bool, const unsigned int, - const unsigned int, + const DataSetDescriptor, const std::vector &, InternalData &, std::vector > &, diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index fcfc88ae3a..c7ebc88290 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -137,7 +138,9 @@ MappingCartesian::compute_fill (const typename DoFHandler::cell_iterat UpdateFlags update_flags(data.current_update_flags()); const unsigned int npts = quadrature_points.size (); - unsigned int offset = 0; + + typedef typename internal::DataSetDescriptor DataSetDescriptor; + unsigned int offset = DataSetDescriptor::cell().offset(); if (face_no != invalid_face_number) { @@ -150,14 +153,20 @@ MappingCartesian::compute_fill (const typename DoFHandler::cell_iterat if (sub_no == invalid_face_number) // called from FEFaceValues - offset = face_no * quadrature_points.size(); + offset = DataSetDescriptor::face (cell, face_no, + quadrature_points.size()).offset(); else { - // called from FESubfaceValue + // 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 = (face_no * GeometryInfo::subfaces_per_face + sub_no) - * quadrature_points.size(); + offset = DataSetDescriptor::sub_face (cell, face_no, sub_no, + quadrature_points.size()).offset(); } } else diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index 47543ddd12..0b8542ac44 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -226,7 +226,7 @@ MappingQ::get_face_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { InternalData *data = new InternalData(n_shape_functions); - Quadrature q (QProjector::project_to_all_faces(quadrature)); + const Quadrature q (QProjector::project_to_all_faces(quadrature)); this->compute_face_data (update_flags, q, quadrature.n_quadrature_points, *data); if (!use_mapping_q_on_all_cells) @@ -244,7 +244,7 @@ MappingQ::get_subface_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { InternalData *data = new InternalData(n_shape_functions); - Quadrature q (QProjector::project_to_all_subfaces(quadrature)); + const Quadrature q (QProjector::project_to_all_subfaces(quadrature)); this->compute_face_data (update_flags, q, quadrature.n_quadrature_points, *data); if (!use_mapping_q_on_all_cells) @@ -327,9 +327,6 @@ MappingQ::fill_fe_face_values (const typename DoFHandler::cell_iterato data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - const unsigned int npts=q.n_quadrature_points; - const unsigned int offset=face_no*npts; - // depending on this result, use // this or the other data object // for the mapping @@ -339,8 +336,12 @@ MappingQ::fill_fe_face_values (const typename DoFHandler::cell_iterato else p_data=&data; + const unsigned int n_q_points=q.n_quadrature_points; this->compute_fill_face (cell, face_no, false, - npts, offset, q.get_weights(), + n_q_points, + MappingQ1::DataSetDescriptor:: + face (cell, face_no, n_q_points), + q.get_weights(), *p_data, quadrature_points, JxW_values, exterior_forms, normal_vectors); @@ -382,10 +383,6 @@ MappingQ::fill_fe_subface_values (const typename DoFHandler::cell_iter data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - const unsigned int npts=q.n_quadrature_points; - const unsigned int offset= - (face_no*GeometryInfo::subfaces_per_face + sub_no)*npts; - // depending on this result, use // this or the other data object // for the mapping @@ -395,8 +392,12 @@ MappingQ::fill_fe_subface_values (const typename DoFHandler::cell_iter else p_data=&data; + const unsigned int n_q_points=q.n_quadrature_points; this->compute_fill_face (cell, face_no, true, - npts, offset, q.get_weights(), + n_q_points, + MappingQ1::DataSetDescriptor:: + sub_face (cell, face_no, sub_no, n_q_points), + q.get_weights(), *p_data, quadrature_points, JxW_values, exterior_forms, normal_vectors); diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index c342fd7f8a..0140474a97 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -25,6 +25,110 @@ #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) + { + 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->get_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) + { + switch (dim) + { + case 1: + case 2: + return ((face_no * GeometryInfo::faces_per_cell + + subface_no) + * n_quadrature_points); + + // for 3d, same as above: + case 3: + return (((face_no * GeometryInfo::faces_per_cell + + subface_no) + + (cell->get_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 const unsigned int MappingQ1::n_shape_functions; @@ -291,7 +395,7 @@ MappingQ1::compute_data (const UpdateFlags update_flags, const unsigned int n_original_q_points, InternalData &data) const { - const unsigned int npts = q.n_quadrature_points; + const unsigned int n_q_points = q.n_quadrature_points; data.update_once = update_once(update_flags); data.update_each = update_each(update_flags); @@ -300,10 +404,10 @@ MappingQ1::compute_data (const UpdateFlags update_flags, const UpdateFlags flags(data.update_flags); if (flags & update_transformation_values) - data.shape_values.resize(data.n_shape_functions * npts); + data.shape_values.resize(data.n_shape_functions * n_q_points); if (flags & update_transformation_gradients) - data.shape_derivatives.resize(data.n_shape_functions * npts); + data.shape_derivatives.resize(data.n_shape_functions * n_q_points); if (flags & update_covariant_transformation) data.covariant.resize(n_original_q_points); @@ -438,17 +542,17 @@ MappingQ1::get_subface_data (const UpdateFlags update_flags, template void MappingQ1::compute_fill (const typename DoFHandler::cell_iterator &cell, - const unsigned int npts, - const unsigned int offset, - InternalData &data, + const unsigned int n_q_points, + const DataSetDescriptor data_set, + InternalData &data, std::vector > &quadrature_points) const { const UpdateFlags update_flags(data.current_update_flags()); if (update_flags & update_q_points) { - Assert (quadrature_points.size() == npts, - ExcDimensionMismatch(quadrature_points.size(), npts)); + Assert (quadrature_points.size() == n_q_points, + ExcDimensionMismatch(quadrature_points.size(), n_q_points)); std::fill(quadrature_points.begin(), quadrature_points.end(), Point()); @@ -456,14 +560,14 @@ MappingQ1::compute_fill (const typename DoFHandler::cell_iterator &cel if (update_flags & update_covariant_transformation) { - Assert (data.covariant.size() == npts, - ExcDimensionMismatch(data.covariant.size(), npts)); + Assert (data.covariant.size() == n_q_points, + ExcDimensionMismatch(data.covariant.size(), n_q_points)); } if (update_flags & update_contravariant_transformation) { - Assert (data.contravariant.size() == npts, - ExcDimensionMismatch(data.contravariant.size(), npts)); + Assert (data.contravariant.size() == n_q_points, + ExcDimensionMismatch(data.contravariant.size(), n_q_points)); std::fill(data.contravariant.begin(), data.contravariant.end(), Tensor<2,dim>()); @@ -472,8 +576,8 @@ MappingQ1::compute_fill (const typename DoFHandler::cell_iterator &cel if (update_flags & update_jacobian_grads) { Assert(false, ExcNotImplemented()); -// Assert (covariant_grads.size () == npts, -// ExcDimensionMismatch(covariant_grads.size(), npts)); +// Assert (covariant_grads.size () == n_q_points, +// ExcDimensionMismatch(covariant_grads.size(), n_q_points)); } // if necessary, recompute the @@ -498,19 +602,20 @@ MappingQ1::compute_fill (const typename DoFHandler::cell_iterator &cel // first compute quadrature points if (update_flags & update_q_points) - for (unsigned int point=0; point::compute_fill (const typename DoFHandler::cell_iterator &cel // covariant transformation // matrices if (update_flags & update_covariant_transformation) - for (unsigned int point=0; point::fill_fe_values (const typename DoFHandler::cell_iterator &c Assert(data_ptr!=0, ExcInternalError()); InternalData &data=*data_ptr; - const unsigned int npts=q.n_quadrature_points; + const unsigned int n_q_points=q.n_quadrature_points; - compute_fill (cell, - npts, - 0, - data, - quadrature_points); + compute_fill (cell, n_q_points, DataSetDescriptor::cell (), + data, quadrature_points); const UpdateFlags update_flags(data.current_update_flags()); @@ -565,9 +667,9 @@ MappingQ1::fill_fe_values (const typename DoFHandler::cell_iterator &c // by Jacobian determinants if (update_flags & update_JxW_values) { - Assert (JxW_values.size() == npts, - ExcDimensionMismatch(JxW_values.size(), npts)); - for (unsigned int point=0; point::compute_fill_face (const typename DoFHandler::cell_iterator &cell, const unsigned int face_no, const bool is_subface, - const unsigned int npts, - const unsigned int offset, + const unsigned int n_q_points, + const DataSetDescriptor data_set, const std::vector &weights, InternalData &data, std::vector > &quadrature_points, @@ -589,24 +691,20 @@ MappingQ1::compute_fill_face (const typename DoFHandler::cell_iterator std::vector > &boundary_forms, std::vector > &normal_vectors) const { - compute_fill (cell, - npts, - offset, - data, - quadrature_points); + compute_fill (cell, n_q_points, data_set, data, quadrature_points); const UpdateFlags update_flags(data.current_update_flags()); if (update_flags & update_boundary_forms) { - Assert (boundary_forms.size()==npts, - ExcDimensionMismatch(boundary_forms.size(), npts)); + Assert (boundary_forms.size()==n_q_points, + ExcDimensionMismatch(boundary_forms.size(), n_q_points)); if (update_flags & update_normal_vectors) - Assert (normal_vectors.size()==npts, - ExcDimensionMismatch(normal_vectors.size(), npts)); + Assert (normal_vectors.size()==n_q_points, + ExcDimensionMismatch(normal_vectors.size(), n_q_points)); if (update_flags & update_JxW_values) - Assert (JxW_values.size() == npts, - ExcDimensionMismatch(JxW_values.size(), npts)); + Assert (JxW_values.size() == n_q_points, + ExcDimensionMismatch(JxW_values.size(), n_q_points)); Assert (data.aux[0].size() <= data.unit_tangentials[face_no].size(), @@ -665,7 +763,7 @@ MappingQ1::compute_fill_face (const typename DoFHandler::cell_iterator { JxW_values[i] = f * weights[i]; if (is_subface) - JxW_values[i]/=GeometryInfo::subfaces_per_face; + JxW_values[i] /= GeometryInfo::subfaces_per_face; } if (update_flags & update_normal_vectors) normal_vectors[i] = boundary_forms[i] / f; @@ -689,12 +787,11 @@ MappingQ1::fill_fe_face_values (const typename DoFHandler::cell_iterat Assert(data_ptr!=0, ExcInternalError()); InternalData &data=*data_ptr; - const unsigned int npts=q.n_quadrature_points; - const unsigned int offset=face_no*npts; + const unsigned int n_q_points = q.n_quadrature_points; compute_fill_face (cell, face_no, false, - npts, - offset, + n_q_points, + DataSetDescriptor::face (cell, face_no, n_q_points), q.get_weights(), data, quadrature_points, @@ -720,13 +817,11 @@ MappingQ1::fill_fe_subface_values (const typename DoFHandler::cell_ite Assert(data_ptr!=0, ExcInternalError()); InternalData &data=*data_ptr; - const unsigned int npts=q.n_quadrature_points; - const unsigned int offset= - (face_no*GeometryInfo::subfaces_per_face + sub_no)*npts; + const unsigned int n_q_points = q.n_quadrature_points; compute_fill_face (cell, face_no, true, - npts, - offset, + n_q_points, + DataSetDescriptor::sub_face (cell, face_no, sub_no, n_q_points), q.get_weights(), data, quadrature_points, @@ -744,7 +839,7 @@ MappingQ1<1>::compute_fill_face (const DoFHandler<1>::cell_iterator &, const unsigned int, const bool, const unsigned int, - const unsigned int, + const DataSetDescriptor, const std::vector &, InternalData &, std::vector > &, -- 2.39.5