From f013eadc1477b049f426bb090672efdc09040ede Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Thu, 29 Aug 2019 13:35:56 +0200 Subject: [PATCH] Split MappingCartesian::compute_fill into several functions. Presently MappingCartesian::compute_fill is called from the various fill_fe*_values functions to compute quadrature points and normals. compute_fill is called with dummy-normals in fill_fe_values and with invalid face indices whenever face_index and subface_index has no meaning. This makes the code hard to understand and the asserts that check the arguments complicated. Split the compute_fill function into several smaller functions to avoid passing around dummy or invalid data. --- include/deal.II/fe/mapping_cartesian.h | 75 +++++++- source/fe/mapping_cartesian.cc | 242 ++++++++++++++----------- 2 files changed, 199 insertions(+), 118 deletions(-) diff --git a/include/deal.II/fe/mapping_cartesian.h b/include/deal.II/fe/mapping_cartesian.h index 6a51c77cc9..7b18a9e518 100644 --- a/include/deal.II/fe/mapping_cartesian.h +++ b/include/deal.II/fe/mapping_cartesian.h @@ -19,6 +19,8 @@ #include +#include + #include #include @@ -255,24 +257,77 @@ private: * @} */ + /** + * Update the cell_extents field of the incoming InternalData object with the + * size of the incoming cell. + */ + void + update_cell_extents( + const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const InternalData & data) const; + + /** + * Compute the quadrature points if the UpdateFlags of the incoming + * InternalData object say that they should be updated. + * + * Called from fill_fe_values. + */ + void + maybe_update_cell_quadrature_points( + const typename Triangulation::cell_iterator &cell, + const InternalData & data, + std::vector> &quadrature_points) const; + + /** + * Compute the quadrature points if the UpdateFlags of the incoming + * InternalData object say that they should be updated. + * + * Called from fill_fe_face_values. + */ + void + maybe_update_face_quadrature_points( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const InternalData & data, + std::vector> &quadrature_points) const; + /** + * Compute the quadrature points if the UpdateFlags of the incoming + * InternalData object say that they should be updated. + * + * Called from fill_fe_subface_values. + */ + void + maybe_update_subface_quadrature_points( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const InternalData & data, + std::vector> &quadrature_points) const; /** - * Do the computation for the fill_* functions. + * Transform quadrature points in InternalData to real space by scaling unit + * coordinates with cell_extends in each direction. + * + * Called from the various maybe_update_*_quadrature_points functions. */ void - compute_fill(const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - const CellSimilarity::Similarity cell_similarity, - const InternalData & data, - std::vector> & quadrature_points, - std::vector> & normal_vectors) const; + transform_quadrature_points( + const typename Triangulation::cell_iterator &cell, + const InternalData & data, + const typename QProjector::DataSetDescriptor & offset, + std::vector> &quadrature_points) const; /** - * Value to indicate that a given face or subface number is invalid. + * Compute the normal vectors if the UpdateFlags of the incoming InternalData + * object say that they should be updated. */ - static const unsigned int invalid_face_number = numbers::invalid_unsigned_int; + void + maybe_update_normal_vectors( + const unsigned int face_no, + const InternalData & data, + std::vector> &normal_vectors) const; }; /*@}*/ diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc index 2b084f528d..96aa3d3574 100644 --- a/source/fe/mapping_cartesian.cc +++ b/source/fe/mapping_cartesian.cc @@ -38,10 +38,6 @@ DEAL_II_NAMESPACE_OPEN -template -const unsigned int MappingCartesian::invalid_face_number; - - template MappingCartesian::InternalData::InternalData( @@ -164,64 +160,18 @@ MappingCartesian::get_subface_data( template void -MappingCartesian::compute_fill( +MappingCartesian::update_cell_extents( const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, const CellSimilarity::Similarity cell_similarity, - const InternalData & data, - std::vector> & quadrature_points, - std::vector> &normal_vectors) const + const InternalData & data) const { - const UpdateFlags update_flags = data.update_each; - - // some more sanity checks - if (face_no != invalid_face_number) - { - // Add 1 on both sides of - // assertion to avoid compiler - // warning about testing - // unsigned int < 0 in 1d. - Assert(face_no + 1 < GeometryInfo::faces_per_cell + 1, - ExcIndexRange(face_no, 0, GeometryInfo::faces_per_cell)); - - // We would like to check for - // sub_no < cell->face(face_no)->n_children(), - // but unfortunately the current - // function is also called for - // faces without children (see - // tests/fe/mapping.cc). Therefore, - // we must use following workaround - // of two separate assertions - Assert((sub_no == invalid_face_number) || - cell->face(face_no)->has_children() || - (sub_no + 1 < GeometryInfo::max_children_per_face + 1), - ExcIndexRange(sub_no, - 0, - GeometryInfo::max_children_per_face)); - Assert((sub_no == invalid_face_number) || - !cell->face(face_no)->has_children() || - (sub_no < cell->face(face_no)->n_children()), - ExcIndexRange(sub_no, 0, cell->face(face_no)->n_children())); - } - else - // invalid face number, so - // subface should be invalid as - // well - Assert(sub_no == invalid_face_number, ExcInternalError()); - - // let @p{start} be the origin of a - // local coordinate system. it is - // chosen as the (lower) left - // vertex - const Point start = cell->vertex(0); - // Compute start point and sizes // along axes. Strange vertex // numbering makes this complicated // again. if (cell_similarity != CellSimilarity::translation) { + const Point &start = cell->vertex(0); switch (dim) { case 1: @@ -240,47 +190,121 @@ MappingCartesian::compute_fill( Assert(false, ExcNotImplemented()); } } +} + - // transform quadrature point. this - // is obtained simply by scaling - // unit coordinates with lengths in - // each direction - if (update_flags & update_quadrature_points) +template +void +MappingCartesian::maybe_update_cell_quadrature_points( + const typename Triangulation::cell_iterator &cell, + const InternalData & data, + std::vector> &quadrature_points) const +{ + if (data.update_each & update_quadrature_points) { 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), - cell->face_flip(face_no), - cell->face_rotation(face_no), - quadrature_points.size()) : - // called from FESubfaceValues - QProjector::DataSetDescriptor::subface( - face_no, - sub_no, - cell->face_orientation(face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - quadrature_points.size(), - cell->subface_case(face_no)))); - - for (unsigned int i = 0; i < quadrature_points.size(); ++i) - { - quadrature_points[i] = start; - for (unsigned int d = 0; d < dim; ++d) - quadrature_points[i](d) += - data.cell_extents[d] * data.quadrature_points[i + offset](d); - } + QProjector::DataSetDescriptor::cell(); + + transform_quadrature_points(cell, data, offset, quadrature_points); + } +} + + + +template +void +MappingCartesian::maybe_update_face_quadrature_points( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const InternalData & data, + std::vector> &quadrature_points) const +{ + AssertIndexRange(face_no, GeometryInfo::faces_per_cell); + + if (data.update_each & update_quadrature_points) + { + const typename QProjector::DataSetDescriptor offset = + QProjector::DataSetDescriptor::face(face_no, + cell->face_orientation( + face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature_points.size()); + + + transform_quadrature_points(cell, data, offset, quadrature_points); + } +} + + + +template +void +MappingCartesian::maybe_update_subface_quadrature_points( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const InternalData & data, + std::vector> &quadrature_points) const +{ + AssertIndexRange(face_no, GeometryInfo::faces_per_cell); + AssertIndexRange(sub_no, GeometryInfo::max_children_per_face); + if (cell->face(face_no)->has_children()) + { + AssertIndexRange(sub_no, cell->face(face_no)->n_children()); + } + + if (data.update_each & update_quadrature_points) + { + const typename QProjector::DataSetDescriptor offset = + QProjector::DataSetDescriptor::subface( + face_no, + sub_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature_points.size(), + cell->subface_case(face_no)); + + transform_quadrature_points(cell, data, offset, quadrature_points); + } +} + + + +template +void +MappingCartesian::transform_quadrature_points( + const typename Triangulation::cell_iterator &cell, + const InternalData & data, + const typename QProjector::DataSetDescriptor & offset, + std::vector> &quadrature_points) const +{ + // let @p{start} be the origin of a local coordinate system. it is chosen + // as the (lower) left vertex + const Point &start = cell->vertex(0); + + for (unsigned int i = 0; i < quadrature_points.size(); ++i) + { + quadrature_points[i] = start; + for (unsigned int d = 0; d < dim; ++d) + quadrature_points[i](d) += + data.cell_extents[d] * data.quadrature_points[i + offset](d); } +} + +template +void +MappingCartesian::maybe_update_normal_vectors( + const unsigned int face_no, + const InternalData & data, + std::vector> &normal_vectors) const +{ // compute normal vectors. All normals on a face have the same value. - if (update_flags & update_normal_vectors) + if (data.update_each & update_normal_vectors) { Assert(face_no < GeometryInfo::faces_per_cell, ExcInternalError()); std::fill(normal_vectors.begin(), @@ -307,15 +331,12 @@ MappingCartesian::fill_fe_values( ExcInternalError()); const InternalData &data = static_cast(internal_data); - std::vector> dummy; - compute_fill(cell, - invalid_face_number, - invalid_face_number, - cell_similarity, - data, - output_data.quadrature_points, - dummy); + update_cell_extents(cell, cell_similarity, data); + + maybe_update_cell_quadrature_points(cell, + data, + output_data.quadrature_points); // compute Jacobian determinant. all values are equal and are the // product of the local lengths in each coordinate direction @@ -420,13 +441,14 @@ MappingCartesian::fill_fe_face_values( ExcInternalError()); const InternalData &data = static_cast(internal_data); - compute_fill(cell, - face_no, - invalid_face_number, - CellSimilarity::none, - data, - output_data.quadrature_points, - output_data.normal_vectors); + update_cell_extents(cell, CellSimilarity::none, data); + + maybe_update_face_quadrature_points(cell, + face_no, + data, + output_data.quadrature_points); + + maybe_update_normal_vectors(face_no, data, output_data.normal_vectors); // first compute Jacobian determinant, which is simply the product // of the local lengths since the jacobian is diagonal @@ -523,13 +545,12 @@ MappingCartesian::fill_fe_subface_values( ExcInternalError()); const InternalData &data = static_cast(internal_data); - compute_fill(cell, - face_no, - subface_no, - CellSimilarity::none, - data, - output_data.quadrature_points, - output_data.normal_vectors); + update_cell_extents(cell, CellSimilarity::none, data); + + maybe_update_subface_quadrature_points( + cell, face_no, subface_no, data, output_data.quadrature_points); + + maybe_update_normal_vectors(face_no, data, output_data.normal_vectors); // first compute Jacobian determinant, which is simply the product // of the local lengths since the jacobian is diagonal @@ -898,6 +919,7 @@ MappingCartesian::transform( } + template void MappingCartesian::transform( @@ -935,6 +957,8 @@ MappingCartesian::transform( } } + + template void MappingCartesian::transform( @@ -1022,6 +1046,7 @@ MappingCartesian::transform( } + template Point MappingCartesian::transform_unit_to_real_cell( @@ -1090,6 +1115,7 @@ MappingCartesian::transform_real_to_unit_cell( } + template std::unique_ptr> MappingCartesian::clone() const -- 2.39.5