From 256a8c0a7e4f54a4699435b0942ea3dbd4ebb9f6 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth <bangerth@math.tamu.edu> Date: Mon, 27 Jul 2015 17:06:51 -0500 Subject: [PATCH] Remove MappingQ1::compute_fill and simplify MappingQ1::compute_fill_face. Also move around a few functions in the .cc file to make the flow easier to read. The old style used rather generic function names that don't really give away what concretely a function does. --- include/deal.II/fe/mapping_q1.h | 13 +- source/fe/mapping_q.cc | 6 +- source/fe/mapping_q1.cc | 290 ++++++++++++++++++-------------- 3 files changed, 169 insertions(+), 140 deletions(-) diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h index f2b6f05242..81071b09a7 100644 --- a/include/deal.II/fe/mapping_q1.h +++ b/include/deal.II/fe/mapping_q1.h @@ -386,25 +386,14 @@ public: const unsigned int n_orig_q_points, InternalData &data) const; - /** - * Do the computation for the <tt>fill_*</tt> functions. - */ - void compute_fill (const typename Triangulation<dim,spacedim>::cell_iterator &cell, - const unsigned int npts, - const DataSetDescriptor data_set, - const CellSimilarity::Similarity cell_similarity, - const InternalData &data, - std::vector<Point<spacedim> > &quadrature_points) const; - /** * Do the computation for the <tt>fill_*</tt> functions. */ void compute_fill_face (const typename Triangulation<dim,spacedim>::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, - const unsigned int npts, const DataSetDescriptor data_set, - const std::vector<double> &weights, + const Quadrature<dim-1> &quadrature, const InternalData &internal_data, FEValuesData<dim,spacedim> &output_data) const; diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 034c0540c8..edb1827ae4 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -334,14 +334,13 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator & const unsigned int n_q_points = quadrature.size(); this->compute_fill_face (cell, face_no, numbers::invalid_unsigned_int, - n_q_points, QProjector<dim>::DataSetDescriptor:: face (face_no, cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), n_q_points), - quadrature.get_weights(), + quadrature, *p_data, output_data); } @@ -383,7 +382,6 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato const unsigned int n_q_points = quadrature.size(); this->compute_fill_face (cell, face_no, subface_no, - n_q_points, QProjector<dim>::DataSetDescriptor:: subface (face_no, subface_no, cell->face_orientation(face_no), @@ -391,7 +389,7 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato cell->face_rotation(face_no), n_q_points, cell->subface_case(face_no)), - quadrature.get_weights(), + quadrature, *p_data, output_data); } diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index b5bc80917d..7e090cfb08 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -603,121 +603,122 @@ MappingQ1<dim,spacedim>::get_subface_data (const UpdateFlags update_flags, -template<int dim, int spacedim> -void -MappingQ1<dim,spacedim>::compute_fill (const typename Triangulation<dim,spacedim>::cell_iterator &cell, - const unsigned int n_q_points, - const DataSetDescriptor data_set, - const CellSimilarity::Similarity cell_similarity, - const InternalData &data, - std::vector<Point<spacedim> > &quadrature_points) const +namespace internal { - const UpdateFlags update_flags(data.current_update_flags()); - - // if necessary, recompute the - // support points of the - // transformation of this cell - // (note that we need to first - // check the triangulation pointer, - // since otherwise the second test - // might trigger an exception if - // the triangulations are not the - // same) - if ((data.mapping_support_points.size() == 0) - || - (&cell->get_triangulation() != - &data.cell_of_current_support_points->get_triangulation()) - || - (cell != data.cell_of_current_support_points)) - { - compute_mapping_support_points(cell, data.mapping_support_points); - data.cell_of_current_support_points = cell; - } - - // first compute quadrature points - if (update_flags & update_quadrature_points) - { - AssertDimension (quadrature_points.size(), n_q_points); - - for (unsigned int point=0; point<n_q_points; ++point) - { - const double *shape = &data.shape(point+data_set,0); - Point<spacedim> result = (shape[0] * - data.mapping_support_points[0]); - for (unsigned int k=1; k<data.n_shape_functions; ++k) - for (unsigned int i=0; i<spacedim; ++i) - result[i] += shape[k] * data.mapping_support_points[k][i]; - quadrature_points[point] = result; - } - } - - // then Jacobians - if (update_flags & update_contravariant_transformation) + namespace + { + /** + * Compute the locations of quadrature points on the object described by + * the first argument (and the cell for which the mapping support points + * have already been set), but only if the update_flags of the @p data + * argument indicate so. + */ + template <int dim, int spacedim> + void + maybe_compute_q_points (const typename dealii::MappingQ1<dim,spacedim>::DataSetDescriptor data_set, + const typename dealii::MappingQ1<dim,spacedim>::InternalData &data, + std::vector<Point<spacedim> > &quadrature_points) { - AssertDimension (data.contravariant.size(), n_q_points); + const UpdateFlags update_flags(data.current_update_flags()); - // if the current cell is just a - // translation of the previous one, no - // need to recompute jacobians... - if (cell_similarity != CellSimilarity::translation) + if (update_flags & update_quadrature_points) { - std::fill(data.contravariant.begin(), data.contravariant.end(), - DerivativeForm<1,dim,spacedim>()); - - Assert (data.n_shape_functions > 0, ExcInternalError()); - const Tensor<1,spacedim> *supp_pts = - &data.mapping_support_points[0]; - - for (unsigned int point=0; point<n_q_points; ++point) + for (unsigned int point=0; point<quadrature_points.size(); ++point) { - const Tensor<1,dim> *data_derv = - &data.derivative(point+data_set, 0); - - double result [spacedim][dim]; - - // peel away part of sum to avoid zeroing the - // entries and adding for the first time - for (unsigned int i=0; i<spacedim; ++i) - for (unsigned int j=0; j<dim; ++j) - result[i][j] = data_derv[0][j] * supp_pts[0][i]; + const double *shape = &data.shape(point+data_set,0); + Point<spacedim> result = (shape[0] * + data.mapping_support_points[0]); for (unsigned int k=1; k<data.n_shape_functions; ++k) for (unsigned int i=0; i<spacedim; ++i) - for (unsigned int j=0; j<dim; ++j) - result[i][j] += data_derv[k][j] * supp_pts[k][i]; - - // write result into contravariant data. for - // j=dim in the case dim<spacedim, there will - // never be any nonzero data that arrives in - // here, so it is ok anyway because it was - // initialized to zero at the initialization - for (unsigned int i=0; i<spacedim; ++i) - for (unsigned int j=0; j<dim; ++j) - data.contravariant[point][i][j] = result[i][j]; + result[i] += shape[k] * data.mapping_support_points[k][i]; + quadrature_points[point] = result; } } } - if (update_flags & update_covariant_transformation) - { - AssertDimension (data.covariant.size(), n_q_points); - if (cell_similarity != CellSimilarity::translation) - for (unsigned int point=0; point<n_q_points; ++point) - { - data.covariant[point] = (data.contravariant[point]).covariant_form(); - } - } - if (update_flags & update_volume_elements) - if (cell_similarity != CellSimilarity::translation) - for (unsigned int point=0; point<n_q_points; ++point) - data.volume_elements[point] = data.contravariant[point].determinant(); + /** + * Update the co- and contravariant matrices as well as their determinant, for the cell + * described stored in the data object, but only if the update_flags of the @p data + * argument indicate so. + * + * Skip the computation if possible as indicated by the first argument. + */ + template <int dim, int spacedim> + void + maybe_update_Jacobians (const CellSimilarity::Similarity cell_similarity, + const typename dealii::MappingQ1<dim,spacedim>::DataSetDescriptor data_set, + const typename dealii::MappingQ1<dim,spacedim>::InternalData &data) + { + const UpdateFlags update_flags(data.current_update_flags()); + if (update_flags & update_contravariant_transformation) + // if the current cell is just a + // translation of the previous one, no + // need to recompute jacobians... + if (cell_similarity != CellSimilarity::translation) + { + const unsigned int n_q_points = data.contravariant.size(); + + std::fill(data.contravariant.begin(), data.contravariant.end(), + DerivativeForm<1,dim,spacedim>()); + + Assert (data.n_shape_functions > 0, ExcInternalError()); + const Tensor<1,spacedim> *supp_pts = + &data.mapping_support_points[0]; + + for (unsigned int point=0; point<n_q_points; ++point) + { + const Tensor<1,dim> *data_derv = + &data.derivative(point+data_set, 0); + + double result [spacedim][dim]; + + // peel away part of sum to avoid zeroing the + // entries and adding for the first time + for (unsigned int i=0; i<spacedim; ++i) + for (unsigned int j=0; j<dim; ++j) + result[i][j] = data_derv[0][j] * supp_pts[0][i]; + for (unsigned int k=1; k<data.n_shape_functions; ++k) + for (unsigned int i=0; i<spacedim; ++i) + for (unsigned int j=0; j<dim; ++j) + result[i][j] += data_derv[k][j] * supp_pts[k][i]; + + // write result into contravariant data. for + // j=dim in the case dim<spacedim, there will + // never be any nonzero data that arrives in + // here, so it is ok anyway because it was + // initialized to zero at the initialization + for (unsigned int i=0; i<spacedim; ++i) + for (unsigned int j=0; j<dim; ++j) + data.contravariant[point][i][j] = result[i][j]; + } + } + + if (update_flags & update_covariant_transformation) + if (cell_similarity != CellSimilarity::translation) + { + const unsigned int n_q_points = data.contravariant.size(); + for (unsigned int point=0; point<n_q_points; ++point) + { + data.covariant[point] = (data.contravariant[point]).covariant_form(); + } + } + if (update_flags & update_volume_elements) + if (cell_similarity != CellSimilarity::translation) + { + const unsigned int n_q_points = data.contravariant.size(); + for (unsigned int point=0; point<n_q_points; ++point) + data.volume_elements[point] = data.contravariant[point].determinant(); + } + + } + } } - template<int dim, int spacedim> void MappingQ1<dim,spacedim>::compute_mapping_support_points( @@ -748,10 +749,27 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell, const unsigned int n_q_points=quadrature.size(); - compute_fill (cell, n_q_points, DataSetDescriptor::cell (), cell_similarity, - data, - output_data.quadrature_points); + // if necessary, recompute the support points of the transformation of this cell + // (note that we need to first check the triangulation pointer, since otherwise + // the second test might trigger an exception if the triangulations are not the + // same) + if ((data.mapping_support_points.size() == 0) + || + (&cell->get_triangulation() != + &data.cell_of_current_support_points->get_triangulation()) + || + (cell != data.cell_of_current_support_points)) + { + compute_mapping_support_points(cell, data.mapping_support_points); + data.cell_of_current_support_points = cell; + } + internal::maybe_compute_q_points<dim,spacedim> (DataSetDescriptor::cell (), + data, + output_data.quadrature_points); + internal::maybe_update_Jacobians<dim,spacedim> (cell_similarity, + DataSetDescriptor::cell (), + data); const UpdateFlags update_flags(data.current_update_flags()); const std::vector<double> &weights=quadrature.get_weights(); @@ -915,16 +933,25 @@ namespace internal { namespace { + /** + * Depending on what information is called for in the update flags of the + * @p data object, compute the various pieces of information that is required + * by the fill_fe_face_values() and fill_fe_subface_values() functions. + * This function simply unifies the work that would be done by + * those two functions. + * + * The resulting data is put into the @p output_data argument. + */ template <int dim, int spacedim> void - compute_fill_face (const dealii::MappingQ1<dim,spacedim> &mapping, - const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const unsigned int n_q_points, - const std::vector<double> &weights, - const typename dealii::MappingQ1<dim,spacedim>::InternalData &data, - FEValuesData<dim,spacedim> &output_data) + maybe_compute_face_data (const dealii::MappingQ1<dim,spacedim> &mapping, + const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const unsigned int n_q_points, + const std::vector<double> &weights, + const typename dealii::MappingQ1<dim,spacedim>::InternalData &data, + FEValuesData<dim,spacedim> &output_data) { const UpdateFlags update_flags(data.current_update_flags()); @@ -1060,21 +1087,38 @@ template<int dim, int spacedim> void MappingQ1<dim,spacedim>:: compute_fill_face (const typename Triangulation<dim,spacedim>::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const unsigned int n_q_points, - const DataSetDescriptor data_set, - const std::vector<double> &weights, + const unsigned int face_no, + const unsigned int subface_no, + const DataSetDescriptor data_set, + const Quadrature<dim-1> &quadrature, const InternalData &internal_data, FEValuesData<dim,spacedim> &output_data) const { - compute_fill (cell, n_q_points, data_set, CellSimilarity::none, - internal_data, - output_data.quadrature_points); - internal::compute_fill_face (*this, - cell, face_no, subface_no, n_q_points, - weights, internal_data, - output_data); + // if necessary, recompute the support points of the transformation of this cell + // (note that we need to first check the triangulation pointer, since otherwise + // the second test might trigger an exception if the triangulations are not the + // same) + if ((internal_data.mapping_support_points.size() == 0) + || + (&cell->get_triangulation() != + &internal_data.cell_of_current_support_points->get_triangulation()) + || + (cell != internal_data.cell_of_current_support_points)) + { + compute_mapping_support_points(cell, internal_data.mapping_support_points); + internal_data.cell_of_current_support_points = cell; + } + + internal::maybe_compute_q_points<dim,spacedim> (data_set, + internal_data, + output_data.quadrature_points); + internal::maybe_update_Jacobians<dim,spacedim> (CellSimilarity::none, + data_set, + internal_data); + internal::maybe_compute_face_data (*this, + cell, face_no, subface_no, quadrature.size(), + quadrature.get_weights(), internal_data, + output_data); } @@ -1096,13 +1140,12 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator & const unsigned int n_q_points = quadrature.size(); compute_fill_face (cell, face_no, numbers::invalid_unsigned_int, - n_q_points, DataSetDescriptor::face (face_no, cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), n_q_points), - quadrature.get_weights(), + quadrature, data, output_data); } @@ -1128,14 +1171,13 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato const unsigned int n_q_points = quadrature.size(); compute_fill_face (cell, face_no, subface_no, - n_q_points, DataSetDescriptor::subface (face_no, subface_no, cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), n_q_points, cell->subface_case(face_no)), - quadrature.get_weights(), + quadrature, data, output_data); } -- 2.39.5