From 6310b2dee3e7526561b785c66a32567bdb93988c Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 8 Jul 2017 14:48:48 -0400 Subject: [PATCH] Reorganize the internal mapping functions. If we want to ultimately use a unity build then we will need to avoid duplicate function names. This commit adds a new policy for avoiding these (a namespacing scheme) and applies the fix to the mapping classes. --- doc/doxygen/headers/coding_conventions.h | 20 + source/fe/mapping_fe_field.cc | 1535 +++++----- source/fe/mapping_manifold.cc | 1182 ++++---- source/fe/mapping_q_generic.cc | 3439 +++++++++++----------- 4 files changed, 3119 insertions(+), 3057 deletions(-) diff --git a/doc/doxygen/headers/coding_conventions.h b/doc/doxygen/headers/coding_conventions.h index 0ae51f4d93..3b2ce05ffe 100644 --- a/doc/doxygen/headers/coding_conventions.h +++ b/doc/doxygen/headers/coding_conventions.h @@ -130,6 +130,26 @@ source files. handle run time errors (like I/O failures) which must be on in any case, not only in debug mode. +
  • Sometimes it makes sense to implement a class by using several + non-member functions that are not part of the public interface and are only + meant to be called in the current source file. Such free functions should be + put in an anonymous namespace structured in the following way: + + 
    +  namespace internal
+  {
+    namespace ClassName
+    {
+      namespace
+      {
+        // free functions go here
+      }
+    }
+  }
+
    +     + where ClassName is the name of the calling class. +
  • Classes and types generally are named using uppercase letters to denote word beginnings (e.g. TriaIterator) — sometimes called camel diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index a91a96b861..e4cad13d82 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -520,774 +520,777 @@ MappingFEField::get_subface_data namespace internal { - namespace + namespace MappingFEField { - /** - * 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 - void - maybe_compute_q_points (const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField::InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector > &quadrature_points) + namespace { - const UpdateFlags update_flags = data.update_each; + /** + * 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 + void + maybe_compute_q_points (const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField::InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector > &quadrature_points) + { + const UpdateFlags update_flags = data.update_each; - if (update_flags & update_quadrature_points) - { - for (unsigned int point=0; point result; - const double *shape = &data.shape(point+data_set,0); + if (update_flags & update_quadrature_points) + { + for (unsigned int point=0; point result; + const double *shape = &data.shape(point+data_set,0); - for (unsigned int k=0; k - void - maybe_update_Jacobians (const CellSimilarity::Similarity cell_similarity, - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField::InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real) - { - const UpdateFlags update_flags = data.update_each; + /** + * 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 + void + maybe_update_Jacobians (const CellSimilarity::Similarity cell_similarity, + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField::InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real) + { + const UpdateFlags update_flags = data.update_each; - // then Jacobians - if (update_flags & update_contravariant_transformation) - { + // then Jacobians + 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(); + // 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(); - Assert (data.n_shape_functions > 0, ExcInternalError()); + Assert (data.n_shape_functions > 0, ExcInternalError()); - for (unsigned int point=0; point *data_derv = - &data.derivative(point+data_set, 0); + for (unsigned int point=0; point *data_derv = + &data.derivative(point+data_set, 0); - Tensor<1, dim> result[spacedim]; + Tensor<1, dim> result[spacedim]; - for (unsigned int k=0; k - void - maybe_update_jacobian_grads (const CellSimilarity::Similarity cell_similarity, - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField::InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector > &jacobian_grads) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_grads) - { - const unsigned int n_q_points = jacobian_grads.size(); + /** + * Update the Hessian of the transformation from unit to real cell, the + * Jacobian gradients. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + void + maybe_update_jacobian_grads (const CellSimilarity::Similarity cell_similarity, + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField::InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector > &jacobian_grads) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_grads) + { + const unsigned int n_q_points = jacobian_grads.size(); - if (cell_similarity != CellSimilarity::translation) - { - for (unsigned int point=0; point *second = - &data.second_derivative(point+data_set, 0); + if (cell_similarity != CellSimilarity::translation) + { + for (unsigned int point=0; point *second = + &data.second_derivative(point+data_set, 0); - DerivativeForm<2,dim,spacedim> result; + DerivativeForm<2,dim,spacedim> result; - for (unsigned int k=0; k - void - maybe_update_jacobian_pushed_forward_grads ( - const CellSimilarity::Similarity cell_similarity, - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField::InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector > &jacobian_pushed_forward_grads ) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_grads) - { - const unsigned int n_q_points = jacobian_pushed_forward_grads.size(); + /** + * Update the Hessian of the transformation from unit to real cell, the + * Jacobian gradients, pushed forward to the real cell coordinates. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + void + maybe_update_jacobian_pushed_forward_grads + (const CellSimilarity::Similarity cell_similarity, + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField::InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector > &jacobian_pushed_forward_grads ) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_grads) + { + const unsigned int n_q_points = jacobian_pushed_forward_grads.size(); - if (cell_similarity != CellSimilarity::translation) - { - double tmp[spacedim][spacedim][spacedim]; - for (unsigned int point=0; point *second = - &data.second_derivative(point+data_set, 0); + if (cell_similarity != CellSimilarity::translation) + { + double tmp[spacedim][spacedim][spacedim]; + for (unsigned int point=0; point *second = + &data.second_derivative(point+data_set, 0); - DerivativeForm<2,dim,spacedim> result; + DerivativeForm<2,dim,spacedim> result; - for (unsigned int k=0; k - void - maybe_update_jacobian_2nd_derivatives (const CellSimilarity::Similarity cell_similarity, - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField::InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector > &jacobian_2nd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_2nd_derivatives) - { - const unsigned int n_q_points = jacobian_2nd_derivatives.size(); + /** + * Update the third derivative of the transformation from unit to real + * cell, the Jacobian hessians. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + void + maybe_update_jacobian_2nd_derivatives (const CellSimilarity::Similarity cell_similarity, + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField::InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector > &jacobian_2nd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_2nd_derivatives) + { + const unsigned int n_q_points = jacobian_2nd_derivatives.size(); - if (cell_similarity != CellSimilarity::translation) - { - for (unsigned int point=0; point *third = - &data.third_derivative(point+data_set, 0); + if (cell_similarity != CellSimilarity::translation) + { + for (unsigned int point=0; point *third = + &data.third_derivative(point+data_set, 0); - DerivativeForm<3,dim,spacedim> result; + DerivativeForm<3,dim,spacedim> result; - for (unsigned int k=0; k - void - maybe_update_jacobian_pushed_forward_2nd_derivatives ( - const CellSimilarity::Similarity cell_similarity, - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField::InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector > &jacobian_pushed_forward_2nd_derivatives ) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) - { - const unsigned int n_q_points = jacobian_pushed_forward_2nd_derivatives.size(); + /** + * Update the third derivative of the transformation from unit to real cell, + * the Jacobian hessians, pushed forward to the real cell coordinates. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + void + maybe_update_jacobian_pushed_forward_2nd_derivatives + (const CellSimilarity::Similarity cell_similarity, + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField::InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector > &jacobian_pushed_forward_2nd_derivatives ) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) + { + const unsigned int n_q_points = jacobian_pushed_forward_2nd_derivatives.size(); - if (cell_similarity != CellSimilarity::translation) - { - double tmp[spacedim][spacedim][spacedim][spacedim]; - for (unsigned int point=0; point *third = - &data.third_derivative(point+data_set, 0); + if (cell_similarity != CellSimilarity::translation) + { + double tmp[spacedim][spacedim][spacedim][spacedim]; + for (unsigned int point=0; point *third = + &data.third_derivative(point+data_set, 0); - DerivativeForm<3,dim,spacedim> result; + DerivativeForm<3,dim,spacedim> result; - for (unsigned int k=0; k - void - maybe_update_jacobian_3rd_derivatives (const CellSimilarity::Similarity cell_similarity, - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField::InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector > &jacobian_3rd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_3rd_derivatives) + /** + * Update the fourth derivative of the transformation from unit to real + * cell, the Jacobian hessian gradients. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + void + maybe_update_jacobian_3rd_derivatives (const CellSimilarity::Similarity cell_similarity, + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField::InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector > &jacobian_3rd_derivatives) { - const unsigned int n_q_points = jacobian_3rd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_3rd_derivatives) { - for (unsigned int point=0; point *fourth = - &data.fourth_derivative(point+data_set, 0); + for (unsigned int point=0; point *fourth = + &data.fourth_derivative(point+data_set, 0); - DerivativeForm<4,dim,spacedim> result; + DerivativeForm<4,dim,spacedim> result; - for (unsigned int k=0; k - void - maybe_update_jacobian_pushed_forward_3rd_derivatives ( - const CellSimilarity::Similarity cell_similarity, - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField::InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector > &jacobian_pushed_forward_3rd_derivatives ) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) + /** + * Update the fourth derivative of the transformation from unit to real cell, + * the Jacobian hessian gradients, pushed forward to the real cell + * coordinates. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + void + maybe_update_jacobian_pushed_forward_3rd_derivatives + (const CellSimilarity::Similarity cell_similarity, + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField::InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector > &jacobian_pushed_forward_3rd_derivatives ) { - const unsigned int n_q_points = jacobian_pushed_forward_3rd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) { - double tmp[spacedim][spacedim][spacedim][spacedim][spacedim]; - for (unsigned int point=0; point *fourth = - &data.fourth_derivative(point+data_set, 0); + double tmp[spacedim][spacedim][spacedim][spacedim][spacedim]; + for (unsigned int point=0; point *fourth = + &data.fourth_derivative(point+data_set, 0); - DerivativeForm<4,dim,spacedim> result; + DerivativeForm<4,dim,spacedim> result; - for (unsigned int k=0; k - void - maybe_compute_face_data (const dealii::Mapping &mapping, - const typename dealii::Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const std::vector &weights, - const typename dealii::MappingFEField::InternalData &data, - internal::FEValues::MappingRelatedData &output_data) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_boundary_forms) + /** + * 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 + void + maybe_compute_face_data (const dealii::Mapping &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const std::vector &weights, + const typename dealii::MappingFEField::InternalData &data, + internal::FEValues::MappingRelatedData &output_data) { - const unsigned int n_q_points = output_data.boundary_forms.size(); - if (update_flags & update_normal_vectors) - AssertDimension (output_data.normal_vectors.size(), n_q_points); - if (update_flags & update_JxW_values) - AssertDimension (output_data.JxW_values.size(), n_q_points); - - // map the unit tangentials to the real cell. checking for d!=dim-1 - // eliminates compiler warnings regarding unsigned int expressions < - // 0. - for (unsigned int d=0; d!=dim-1; ++d) - { - Assert (face_no+GeometryInfo::faces_per_cell*d < - data.unit_tangentials.size(), - ExcInternalError()); - Assert (data.aux[d].size() <= - data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d].size(), - ExcInternalError()); - - mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d]), - mapping_contravariant, - data, - make_array_view(data.aux[d])); - } + const UpdateFlags update_flags = data.update_each; - // if dim==spacedim, we can use the unit tangentials to compute the - // boundary form by simply taking the cross product - if (dim == spacedim) + if (update_flags & update_boundary_forms) { - for (unsigned int i=0; i::faces_per_cell*d < + data.unit_tangentials.size(), + ExcInternalError()); + Assert (data.aux[d].size() <= + data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d].size(), + ExcInternalError()); + + mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d]), + mapping_contravariant, + data, + make_array_view(data.aux[d])); + } + + // if dim==spacedim, we can use the unit tangentials to compute the + // boundary form by simply taking the cross product + if (dim == spacedim) + { + for (unsigned int i=0; i DX_t = - data.contravariant[point].transpose(); + if (dim==2) + { + const DerivativeForm<1,spacedim,dim> DX_t = + data.contravariant[point].transpose(); - Tensor<1, spacedim> cell_normal = - cross_product_3d(DX_t[0], DX_t[1]); - cell_normal /= cell_normal.norm(); + Tensor<1, spacedim> cell_normal = + cross_product_3d(DX_t[0], DX_t[1]); + cell_normal /= cell_normal.norm(); - // then compute the face normal from the face tangent - // and the cell normal: - output_data.boundary_forms[point] = - cross_product_3d(data.aux[0][point], cell_normal); - } + // then compute the face normal from the face tangent + // and the cell normal: + output_data.boundary_forms[point] = + cross_product_3d(data.aux[0][point], cell_normal); + } + } } - } - if (update_flags & (update_normal_vectors | update_JxW_values)) - for (unsigned int i=0; i::subface_ratio( - cell->subface_case(face_no), subface_no); - output_data.JxW_values[i] *= area_ratio; + output_data.JxW_values[i] = output_data.boundary_forms[i].norm() * weights[i]; + + if (subface_no != numbers::invalid_unsigned_int) + { + const double area_ratio=GeometryInfo::subface_ratio( + cell->subface_case(face_no), subface_no); + output_data.JxW_values[i] *= area_ratio; + } } - } - if (update_flags & update_normal_vectors) - output_data.normal_vectors[i] = Point(output_data.boundary_forms[i] / output_data.boundary_forms[i].norm()); - } + if (update_flags & update_normal_vectors) + output_data.normal_vectors[i] = Point(output_data.boundary_forms[i] / output_data.boundary_forms[i].norm()); + } - if (update_flags & update_jacobians) - for (unsigned int point=0; point - void - do_fill_fe_face_values (const dealii::Mapping &mapping, - const typename dealii::Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const typename dealii::QProjector::DataSetDescriptor data_set, - const Quadrature &quadrature, - const typename dealii::MappingFEField::InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - internal::FEValues::MappingRelatedData &output_data) - { - maybe_compute_q_points - (data_set, - data, - fe, fe_mask, fe_to_real, - output_data.quadrature_points); - - maybe_update_Jacobians - (CellSimilarity::none, - data_set, - data, - fe, fe_mask, fe_to_real); - - maybe_update_jacobian_grads - (CellSimilarity::none, - data_set, - data, - fe, fe_mask, fe_to_real, - output_data.jacobian_grads); - - maybe_update_jacobian_pushed_forward_grads - (CellSimilarity::none, - data_set, - data, - fe, fe_mask, fe_to_real, - output_data.jacobian_pushed_forward_grads); - - maybe_update_jacobian_2nd_derivatives - (CellSimilarity::none, - data_set, - data, - fe, fe_mask, fe_to_real, - output_data.jacobian_2nd_derivatives); - - maybe_update_jacobian_pushed_forward_2nd_derivatives - (CellSimilarity::none, - data_set, - data, - fe, fe_mask, fe_to_real, - output_data.jacobian_pushed_forward_2nd_derivatives); - - maybe_update_jacobian_3rd_derivatives - (CellSimilarity::none, - data_set, - data, - fe, fe_mask, fe_to_real, - output_data.jacobian_3rd_derivatives); - - maybe_update_jacobian_pushed_forward_3rd_derivatives - (CellSimilarity::none, - data_set, - data, - fe, fe_mask, fe_to_real, - output_data.jacobian_pushed_forward_3rd_derivatives); - - maybe_compute_face_data - (mapping, - cell, face_no, subface_no, - quadrature.get_weights(), data, - output_data); + /** + * Do the work of MappingFEField::fill_fe_face_values() and + * MappingFEField::fill_fe_subface_values() in a generic way, using the + * 'data_set' to differentiate whether we will work on a face (and if so, + * which one) or subface. + */ + template + void + do_fill_fe_face_values (const dealii::Mapping &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const typename dealii::QProjector::DataSetDescriptor data_set, + const Quadrature &quadrature, + const typename dealii::MappingFEField::InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + internal::FEValues::MappingRelatedData &output_data) + { + maybe_compute_q_points + (data_set, + data, + fe, fe_mask, fe_to_real, + output_data.quadrature_points); + + maybe_update_Jacobians + (CellSimilarity::none, + data_set, + data, + fe, fe_mask, fe_to_real); + + maybe_update_jacobian_grads + (CellSimilarity::none, + data_set, + data, + fe, fe_mask, fe_to_real, + output_data.jacobian_grads); + + maybe_update_jacobian_pushed_forward_grads + (CellSimilarity::none, + data_set, + data, + fe, fe_mask, fe_to_real, + output_data.jacobian_pushed_forward_grads); + + maybe_update_jacobian_2nd_derivatives + (CellSimilarity::none, + data_set, + data, + fe, fe_mask, fe_to_real, + output_data.jacobian_2nd_derivatives); + + maybe_update_jacobian_pushed_forward_2nd_derivatives + (CellSimilarity::none, + data_set, + data, + fe, fe_mask, fe_to_real, + output_data.jacobian_pushed_forward_2nd_derivatives); + + maybe_update_jacobian_3rd_derivatives + (CellSimilarity::none, + data_set, + data, + fe, fe_mask, fe_to_real, + output_data.jacobian_3rd_derivatives); + + maybe_update_jacobian_pushed_forward_3rd_derivatives + (CellSimilarity::none, + data_set, + data, + fe, fe_mask, fe_to_real, + output_data.jacobian_pushed_forward_3rd_derivatives); + + maybe_compute_face_data + (mapping, + cell, face_no, subface_no, + quadrature.get_weights(), data, + output_data); + } + } } } @@ -1318,12 +1321,12 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, update_internal_dofs(cell, data); - internal::maybe_compute_q_points + internal::MappingFEField::maybe_compute_q_points (QProjector::DataSetDescriptor::cell (), data, *fe, fe_mask, fe_to_real, output_data.quadrature_points); - internal::maybe_update_Jacobians + internal::MappingFEField::maybe_update_Jacobians (cell_similarity, QProjector::DataSetDescriptor::cell (), data, *fe, fe_mask, fe_to_real); @@ -1427,35 +1430,35 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, } // calculate derivatives of the Jacobians - internal::maybe_update_jacobian_grads + internal::MappingFEField::maybe_update_jacobian_grads (cell_similarity, QProjector::DataSetDescriptor::cell(), data, *fe, fe_mask, fe_to_real, output_data.jacobian_grads); // calculate derivatives of the Jacobians pushed forward to real cell coordinates - internal::maybe_update_jacobian_pushed_forward_grads + internal::MappingFEField::maybe_update_jacobian_pushed_forward_grads (cell_similarity, QProjector::DataSetDescriptor::cell(), data, *fe, fe_mask, fe_to_real, output_data.jacobian_pushed_forward_grads); // calculate hessians of the Jacobians - internal::maybe_update_jacobian_2nd_derivatives + internal::MappingFEField::maybe_update_jacobian_2nd_derivatives (cell_similarity, QProjector::DataSetDescriptor::cell(), data, *fe, fe_mask, fe_to_real, output_data.jacobian_2nd_derivatives); // calculate hessians of the Jacobians pushed forward to real cell coordinates - internal::maybe_update_jacobian_pushed_forward_2nd_derivatives + internal::MappingFEField::maybe_update_jacobian_pushed_forward_2nd_derivatives (cell_similarity, QProjector::DataSetDescriptor::cell(), data, *fe, fe_mask, fe_to_real, output_data.jacobian_pushed_forward_2nd_derivatives); // calculate gradients of the hessians of the Jacobians - internal::maybe_update_jacobian_3rd_derivatives + internal::MappingFEField::maybe_update_jacobian_3rd_derivatives (cell_similarity, QProjector::DataSetDescriptor::cell(), data, *fe, fe_mask, fe_to_real, @@ -1463,7 +1466,7 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, // calculate gradients of the hessians of the Jacobians pushed forward to real // cell coordinates - internal::maybe_update_jacobian_pushed_forward_3rd_derivatives + internal::MappingFEField::maybe_update_jacobian_pushed_forward_3rd_derivatives (cell_similarity, QProjector::DataSetDescriptor::cell(), data, *fe, fe_mask, fe_to_real, @@ -1491,7 +1494,7 @@ fill_fe_face_values (const typename Triangulation::cell_iterator & update_internal_dofs(cell, data); - internal::do_fill_fe_face_values + internal::MappingFEField::do_fill_fe_face_values (*this, cell, face_no, numbers::invalid_unsigned_int, QProjector::DataSetDescriptor:: @@ -1525,7 +1528,7 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato update_internal_dofs(cell, data); - internal::do_fill_fe_face_values + internal::MappingFEField::do_fill_fe_face_values (*this, cell, face_no, numbers::invalid_unsigned_int, QProjector::DataSetDescriptor:: @@ -1542,101 +1545,107 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato } -namespace +namespace internal { - template - void - transform_fields(const ArrayView > &input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView > &output) + namespace MappingFEField { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), - ExcInternalError()); - const typename MappingFEField::InternalData - &data = static_cast::InternalData &>(mapping_data); - - switch (mapping_type) - { - case mapping_contravariant: + namespace + { + template + void + transform_fields(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) { - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingFEField::InternalData + &data = static_cast::InternalData &>(mapping_data); - for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); - return; - } + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); - Assert (data.update_each & update_volume_elements, - typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); - Assert (rank==1, ExcMessage("Only for rank 1")); - for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + Assert (rank==1, ExcMessage("Only for rank 1")); + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); + //We still allow this operation as in the + //reference cell Derivatives are Tensor + //rather than DerivativeForm + case mapping_covariant: + { + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - for (unsigned int i=0; i - void - transform_differential_forms - (const ArrayView > &input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView > &output) - { + template + void + transform_differential_forms + (const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) + { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), - ExcInternalError()); - const typename MappingFEField::InternalData - &data = static_cast::InternalData &>(mapping_data); + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingFEField::InternalData + &data = static_cast::InternalData &>(mapping_data); - switch (mapping_type) - { - case mapping_covariant: - { - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + switch (mapping_type) + { + case mapping_covariant: + { + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - for (unsigned int i=0; i > &input, { AssertDimension (input.size(), output.size()); - transform_fields(input, mapping_type, mapping_data, output); + internal::MappingFEField::transform_fields(input, mapping_type, mapping_data, output); } @@ -1667,7 +1676,7 @@ transform (const ArrayView > &input, { AssertDimension (input.size(), output.size()); - transform_differential_forms(input, mapping_type, mapping_data, output); + internal::MappingFEField::transform_differential_forms(input, mapping_type, mapping_data, output); } diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 2ac37d2cb9..8949115094 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -349,121 +349,126 @@ MappingManifold::get_subface_data (const UpdateFlags update_flags, namespace internal { - namespace + namespace MappingManifold { - /** - * 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 - void - maybe_compute_q_points (const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingManifold::InternalData &data, - std::vector > &quadrature_points) + namespace { - const UpdateFlags update_flags = data.update_each; + /** + * 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 + void + maybe_compute_q_points + (const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingManifold::InternalData &data, + std::vector > &quadrature_points) + { + const UpdateFlags update_flags = data.update_each; - AssertDimension(data.vertices.size(), GeometryInfo::vertices_per_cell); + AssertDimension(data.vertices.size(), GeometryInfo::vertices_per_cell); - if (update_flags & update_quadrature_points) - { - for (unsigned int point=0; point - get_new_point(data.vertices, - data.cell_manifold_quadrature_weights[point+data_set]); - } - } - } + if (update_flags & update_quadrature_points) + { + for (unsigned int point=0; point + get_new_point(data.vertices, + data.cell_manifold_quadrature_weights[point+data_set]); + } + } + } - /** - * 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. - */ - template - void - maybe_update_Jacobians (const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingManifold::InternalData &data) - { - const UpdateFlags update_flags = data.update_each; + /** + * 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. + */ + template + void + maybe_update_Jacobians + (const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingManifold::InternalData &data) + { + const UpdateFlags update_flags = data.update_each; - if (update_flags & update_contravariant_transformation) - { - const unsigned int n_q_points = data.contravariant.size(); + if (update_flags & update_contravariant_transformation) + { + const unsigned int n_q_points = data.contravariant.size(); - std::fill(data.contravariant.begin(), data.contravariant.end(), - DerivativeForm<1,dim,spacedim>()); + std::fill(data.contravariant.begin(), data.contravariant.end(), + DerivativeForm<1,dim,spacedim>()); - AssertDimension(GeometryInfo::vertices_per_cell, - data.vertices.size()); - for (unsigned int point=0; point &p = data.quad.point(point+data_set); - - // And get its image on the manifold: - const Point P = data.manifold-> - get_new_point(data.vertices, - data.cell_manifold_quadrature_weights[point+data_set]); - - // To compute the Jacobian, we choose dim points aligned - // with the dim reference axes, which are still in the - // given cell, and ask for the tangent vector in these - // directions. Choosing the points is somewhat arbitrary, - // so we try to be smart and we pick points which are - // on the opposite quadrant w.r.t. the evaluation - // point. - for (unsigned int i=0; i ei = Point::unit_vector(i); - const double pi = p[i]; - Assert(pi >=0 && pi <= 1.0, - ExcInternalError("Was expecting a quadrature point " - "inside the unit reference element.")); + AssertDimension(GeometryInfo::vertices_per_cell, + data.vertices.size()); + for (unsigned int point=0; point &p = data.quad.point(point+data_set); + + // And get its image on the manifold: + const Point P = data.manifold-> + get_new_point(data.vertices, + data.cell_manifold_quadrature_weights[point+data_set]); + + // To compute the Jacobian, we choose dim points aligned + // with the dim reference axes, which are still in the + // given cell, and ask for the tangent vector in these + // directions. Choosing the points is somewhat arbitrary, + // so we try to be smart and we pick points which are + // on the opposite quadrant w.r.t. the evaluation + // point. + for (unsigned int i=0; i ei = Point::unit_vector(i); + const double pi = p[i]; + Assert(pi >=0 && pi <= 1.0, + ExcInternalError("Was expecting a quadrature point " + "inside the unit reference element.")); - // In the length L, we store also the direction sign, - // which is positive, if the coordinate is < .5, - const double L = pi > .5 ? -pi: 1-pi; + // In the length L, we store also the direction sign, + // which is positive, if the coordinate is < .5, + const double L = pi > .5 ? -pi: 1-pi; - const Point np(p + L*ei); + const Point np(p + L*ei); - // Get the weights to compute the np point in real space - for (unsigned int j=0; j::vertices_per_cell; ++j) - data.vertex_weights[j] = GeometryInfo::d_linear_shape_function(np, j); + // Get the weights to compute the np point in real space + for (unsigned int j=0; j::vertices_per_cell; ++j) + data.vertex_weights[j] = GeometryInfo::d_linear_shape_function(np, j); - const Point NP= - data.manifold->get_new_point(data.vertices, - data.vertex_weights); + const Point NP= + data.manifold->get_new_point(data.vertices, + data.vertex_weights); - const Tensor<1,spacedim> T = data.manifold->get_tangent_vector(P, NP); + const Tensor<1,spacedim> T = data.manifold->get_tangent_vector(P, NP); - for (unsigned int d=0; d::cell_iterator &cell, data.store_vertices(cell); data.manifold = &(cell->get_manifold()); - internal::maybe_compute_q_points (QProjector::DataSetDescriptor::cell (), - data, - output_data.quadrature_points); + internal::MappingManifold::maybe_compute_q_points + (QProjector::DataSetDescriptor::cell (), + data, + output_data.quadrature_points); - internal::maybe_update_Jacobians (QProjector::DataSetDescriptor::cell (), - data); + internal::MappingManifold::maybe_update_Jacobians + (QProjector::DataSetDescriptor::cell (), + data); const UpdateFlags update_flags = data.update_each; const std::vector &weights=quadrature.get_weights(); @@ -616,580 +623,582 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, namespace internal { - namespace + namespace MappingManifold { - /** - * 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 - void - maybe_compute_face_data (const dealii::MappingManifold &mapping, - const typename dealii::Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const unsigned int n_q_points, - const std::vector &weights, - const typename dealii::MappingManifold::InternalData &data, - internal::FEValues::MappingRelatedData &output_data) + namespace { - const UpdateFlags update_flags = data.update_each; - - if (update_flags & update_boundary_forms) - { - AssertDimension (output_data.boundary_forms.size(), n_q_points); - if (update_flags & update_normal_vectors) - AssertDimension (output_data.normal_vectors.size(), n_q_points); - if (update_flags & update_JxW_values) - AssertDimension (output_data.JxW_values.size(), n_q_points); - - // map the unit tangentials to the real cell. checking for d!=dim-1 - // eliminates compiler warnings regarding unsigned int expressions < - // 0. - for (unsigned int d=0; d!=dim-1; ++d) - { - Assert (face_no+GeometryInfo::faces_per_cell*d < - data.unit_tangentials.size(), - ExcInternalError()); - Assert (data.aux[d].size() <= - data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d].size(), - ExcInternalError()); - - mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d]), - mapping_contravariant, - data, - make_array_view(data.aux[d])); - } + /** + * 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 + void + maybe_compute_face_data + (const dealii::MappingManifold &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const unsigned int n_q_points, + const std::vector &weights, + const typename dealii::MappingManifold::InternalData &data, + internal::FEValues::MappingRelatedData &output_data) + { + const UpdateFlags update_flags = data.update_each; - // if dim==spacedim, we can use the unit tangentials to compute the - // boundary form by simply taking the cross product - if (dim == spacedim) - { - for (unsigned int i=0; i::faces_per_cell*d < + data.unit_tangentials.size(), + ExcInternalError()); + Assert (data.aux[d].size() <= + data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d].size(), + ExcInternalError()); + + mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d]), + mapping_contravariant, + data, + make_array_view(data.aux[d])); + } - for (unsigned int point=0; point DX_t = - data.contravariant[point].transpose(); - - Tensor<1, spacedim> cell_normal = - cross_product_3d(DX_t[0], DX_t[1]); - cell_normal /= cell_normal.norm(); - - // then compute the face normal from the face tangent - // and the cell normal: - output_data.boundary_forms[point] = - cross_product_3d(data.aux[0][point], cell_normal); - + output_data.boundary_forms[i] = + cross_product_2d(data.aux[0][i]); + break; + case 3: + output_data.boundary_forms[i] = + cross_product_3d(data.aux[0][i], data.aux[1][i]); break; - } - default: - Assert (false, ExcNotImplemented()); + Assert(false, ExcNotImplemented()); } - } - } - - if (update_flags & (update_normal_vectors - | update_JxW_values)) - for (unsigned int i=0; i::subface_ratio(cell->subface_case(face_no), - subface_no); - output_data.JxW_values[i] *= area_ratio; - } - } - - if (update_flags & update_normal_vectors) - output_data.normal_vectors[i] = Point(output_data.boundary_forms[i] / - output_data.boundary_forms[i].norm()); - } - - if (update_flags & update_jacobians) - for (unsigned int point=0; point - void - do_fill_fe_face_values (const dealii::MappingManifold &mapping, - const typename dealii::Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const typename QProjector::DataSetDescriptor data_set, - const Quadrature &quadrature, - const typename dealii::MappingManifold::InternalData &data, - internal::FEValues::MappingRelatedData &output_data) - { - data.store_vertices(cell); + break; + } - data.manifold = &cell->face(face_no)->get_manifold(); + case 2: + { + const DerivativeForm<1,spacedim,dim> DX_t = + data.contravariant[point].transpose(); - maybe_compute_q_points (data_set, - data, - output_data.quadrature_points); - maybe_update_Jacobians (data_set, - data); + Tensor<1, spacedim> cell_normal = + cross_product_3d(DX_t[0], DX_t[1]); + cell_normal /= cell_normal.norm(); - maybe_compute_face_data (mapping, - cell, face_no, subface_no, quadrature.size(), - quadrature.get_weights(), data, - output_data); - } - } -} + // then compute the face normal from the face tangent + // and the cell normal: + output_data.boundary_forms[point] = + cross_product_3d(data.aux[0][point], cell_normal); + break; + } + default: + Assert (false, ExcNotImplemented()); + } + } + } -template -void -MappingManifold:: -fill_fe_face_values (const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const Quadrature &quadrature, - const typename Mapping::InternalDataBase &internal_data, - internal::FEValues::MappingRelatedData &output_data) const -{ - // ensure that the following cast is really correct: - Assert ((dynamic_cast(&internal_data) != nullptr), - ExcInternalError()); - const InternalData &data - = static_cast(internal_data); + if (update_flags & (update_normal_vectors + | update_JxW_values)) + for (unsigned int i=0; i::subface_ratio(cell->subface_case(face_no), + subface_no); + output_data.JxW_values[i] *= area_ratio; + } + } - internal::do_fill_fe_face_values (*this, - cell, face_no, numbers::invalid_unsigned_int, - QProjector::DataSetDescriptor::face (face_no, - cell->face_orientation(face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - quadrature.size()), - quadrature, - data, - output_data); -} + if (update_flags & update_normal_vectors) + output_data.normal_vectors[i] = Point(output_data.boundary_forms[i] / + output_data.boundary_forms[i].norm()); + } + if (update_flags & update_jacobians) + for (unsigned int point=0; point -void -MappingManifold:: -fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const Quadrature &quadrature, - const typename Mapping::InternalDataBase &internal_data, - internal::FEValues::MappingRelatedData &output_data) const -{ - // ensure that the following cast is really correct: - Assert ((dynamic_cast(&internal_data) != nullptr), - ExcInternalError()); - const InternalData &data - = static_cast(internal_data); - internal::do_fill_fe_face_values (*this, - cell, face_no, subface_no, - QProjector::DataSetDescriptor::subface (face_no, subface_no, - cell->face_orientation(face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - quadrature.size(), - cell->subface_case(face_no)), - quadrature, - data, - output_data); -} + /** + * Do the work of MappingManifold::fill_fe_face_values() and + * MappingManifold::fill_fe_subface_values() in a generic way, + * using the 'data_set' to differentiate whether we will + * work on a face (and if so, which one) or subface. + */ + template + void + do_fill_fe_face_values + (const dealii::MappingManifold &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const typename QProjector::DataSetDescriptor data_set, + const Quadrature &quadrature, + const typename dealii::MappingManifold::InternalData &data, + internal::FEValues::MappingRelatedData &output_data) + { + data.store_vertices(cell); + data.manifold = &cell->face(face_no)->get_manifold(); + maybe_compute_q_points (data_set, + data, + output_data.quadrature_points); + maybe_update_Jacobians (data_set, + data); -namespace -{ - template - void - transform_fields(const ArrayView > &input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView > &output) - { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), - ExcInternalError()); - const typename MappingManifold::InternalData - &data = static_cast::InternalData &>(mapping_data); + maybe_compute_face_data (mapping, + cell, face_no, subface_no, quadrature.size(), + quadrature.get_weights(), data, + output_data); + } - switch (mapping_type) + template + void + transform_fields(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) { - case mapping_contravariant: - { - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingManifold::InternalData + &data = static_cast::InternalData &>(mapping_data); - for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); - return; - } + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); - Assert (data.update_each & update_volume_elements, - typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); - Assert (rank==1, ExcMessage("Only for rank 1")); - if (rank!=1) - return; - - for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + Assert (rank==1, ExcMessage("Only for rank 1")); + if (rank!=1) + return; + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_covariant_transformation")); + //We still allow this operation as in the + //reference cell Derivatives are Tensor + //rather than DerivativeForm + case mapping_covariant: + { + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - for (unsigned int i=0; i - void - transform_gradients(const ArrayView > &input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView > &output) - { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), - ExcInternalError()); - const typename MappingManifold::InternalData - &data = static_cast::InternalData &>(mapping_data); - switch (mapping_type) + template + void + transform_gradients(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) { - case mapping_contravariant_gradient: - { - Assert (data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - Assert (rank==2, ExcMessage("Only for rank 2")); + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingManifold::InternalData + &data = static_cast::InternalData &>(mapping_data); - for (unsigned int i=0; i A = - apply_transformation(data.contravariant[i], transpose(input[i]) ); - output[i] = apply_transformation(data.covariant[i], A.transpose() ); - } + case mapping_contravariant_gradient: + { + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (rank==2, ExcMessage("Only for rank 2")); - return; - } + for (unsigned int i=0; i A = + apply_transformation(data.contravariant[i], transpose(input[i]) ); + output[i] = apply_transformation(data.covariant[i], A.transpose() ); + } - case mapping_covariant_gradient: - { - Assert (data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - Assert (rank==2, ExcMessage("Only for rank 2")); + return; + } - for (unsigned int i=0; i A = - apply_transformation(data.covariant[i], transpose(input[i]) ); - output[i] = apply_transformation(data.covariant[i], A.transpose() ); - } + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (rank==2, ExcMessage("Only for rank 2")); - return; - } + for (unsigned int i=0; i A = + apply_transformation(data.covariant[i], transpose(input[i]) ); + output[i] = apply_transformation(data.covariant[i], A.transpose() ); + } - case mapping_piola_gradient: - { - Assert (data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - Assert (data.update_each & update_volume_elements, - typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); - Assert (rank==2, ExcMessage("Only for rank 2")); - - for (unsigned int i=0; i A = - apply_transformation(data.covariant[i], input[i] ); - Tensor<2,spacedim> T = - apply_transformation(data.contravariant[i], A.transpose() ); + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + Assert (rank==2, ExcMessage("Only for rank 2")); + + for (unsigned int i=0; i A = + apply_transformation(data.covariant[i], input[i] ); + Tensor<2,spacedim> T = + apply_transformation(data.contravariant[i], A.transpose() ); - output[i] = transpose(T); - output[i] /= data.volume_elements[i]; - } + output[i] = transpose(T); + output[i] /= data.volume_elements[i]; + } - return; - } + return; + } - default: - Assert(false, ExcNotImplemented()); + default: + Assert(false, ExcNotImplemented()); + } } - } - template - void - transform_hessians(const ArrayView > &input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView > &output) - { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), - ExcInternalError()); - const typename MappingManifold::InternalData - &data = static_cast::InternalData &>(mapping_data); - - switch (mapping_type) - { - case mapping_contravariant_hessian: + template + void + transform_hessians(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) { - Assert (data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingManifold::InternalData + &data = static_cast::InternalData &>(mapping_data); - for (unsigned int q=0; q::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + + for (unsigned int q=0; q::ExcAccessToUninitializedField("update_covariant_transformation")); + case mapping_covariant_hessian: + { + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - for (unsigned int q=0; q::ExcAccessToUninitializedField("update_covariant_transformation")); - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - Assert (data.update_each & update_volume_elements, - typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); - - for (unsigned int q=0; q::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + + for (unsigned int q=0; q - void - transform_differential_forms(const ArrayView > &input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView > &output) - { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), - ExcInternalError()); - const typename MappingManifold::InternalData - &data = static_cast::InternalData &>(mapping_data); - switch (mapping_type) + template + void + transform_differential_forms(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) { - case mapping_covariant: - { - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingManifold::InternalData + &data = static_cast::InternalData &>(mapping_data); - for (unsigned int i=0; i::ExcAccessToUninitializedField("update_covariant_transformation")); - return; - } - default: - Assert(false, ExcNotImplemented()); + for (unsigned int i=0; i +void +MappingManifold:: +fill_fe_face_values (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValues::MappingRelatedData &output_data) const +{ + // ensure that the following cast is really correct: + Assert ((dynamic_cast(&internal_data) != nullptr), + ExcInternalError()); + const InternalData &data + = static_cast(internal_data); + + internal::MappingManifold::do_fill_fe_face_values + (*this, + cell, face_no, numbers::invalid_unsigned_int, + QProjector::DataSetDescriptor::face (face_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature.size()), + quadrature, + data, + output_data); +} + + + +template +void +MappingManifold:: +fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValues::MappingRelatedData &output_data) const +{ + // ensure that the following cast is really correct: + Assert ((dynamic_cast(&internal_data) != nullptr), + ExcInternalError()); + const InternalData &data + = static_cast(internal_data); + + internal::MappingManifold::do_fill_fe_face_values + (*this, + cell, face_no, subface_no, + QProjector::DataSetDescriptor::subface (face_no, subface_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature.size(), + cell->subface_case(face_no)), + quadrature, + data, + output_data); +} + + + template void MappingManifold:: @@ -1198,7 +1207,7 @@ transform (const ArrayView > &input, const typename Mapping::InternalDataBase &mapping_data, const ArrayView > &output) const { - transform_fields(input, mapping_type, mapping_data, output); + internal::MappingManifold::transform_fields(input, mapping_type, mapping_data, output); } @@ -1211,7 +1220,7 @@ transform (const ArrayView > &input, const typename Mapping::InternalDataBase &mapping_data, const ArrayView > &output) const { - transform_differential_forms(input, mapping_type, mapping_data, output); + internal::MappingManifold::transform_differential_forms(input, mapping_type, mapping_data, output); } @@ -1227,13 +1236,13 @@ transform (const ArrayView > &input, switch (mapping_type) { case mapping_contravariant: - transform_fields(input, mapping_type, mapping_data, output); + internal::MappingManifold::transform_fields(input, mapping_type, mapping_data, output); return; case mapping_piola_gradient: case mapping_contravariant_gradient: case mapping_covariant_gradient: - transform_gradients(input, mapping_type, mapping_data, output); + internal::MappingManifold::transform_gradients(input, mapping_type, mapping_data, output); return; default: Assert(false, ExcNotImplemented()); @@ -1304,7 +1313,7 @@ transform (const ArrayView > &input, case mapping_piola_hessian: case mapping_contravariant_hessian: case mapping_covariant_hessian: - transform_hessians(input, mapping_type, mapping_data, output); + internal::MappingManifold::transform_hessians(input, mapping_type, mapping_data, output); return; default: Assert(false, ExcNotImplemented()); @@ -1316,4 +1325,3 @@ transform (const ArrayView > &input, DEAL_II_NAMESPACE_CLOSE - diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 97f3d31512..292af5d8b9 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -42,26 +43,28 @@ DEAL_II_NAMESPACE_OPEN -namespace +namespace internal { - template - std::vector - get_dpo_vector (const unsigned int degree) + namespace MappingQGeneric { - std::vector dpo(dim+1, 1U); - for (unsigned int i=1; i + std::vector + get_dpo_vector (const unsigned int degree) + { + std::vector dpo(dim+1, 1U); + for (unsigned int i=1; i renumber (FETools:: - lexicographic_to_hierarchic_numbering ( - FiniteElementData (get_dpo_vector(data.polynomial_degree), 1, - data.polynomial_degree))); + lexicographic_to_hierarchic_numbering + (FiniteElementData (internal::MappingQGeneric::get_dpo_vector + (data.polynomial_degree), 1, data.polynomial_degree))); std::vector values; std::vector > grads; @@ -830,633 +833,501 @@ compute_shape_function_values (const std::vector > &unit_points) } -namespace +namespace internal { - /** - * Compute the support_point_weights_on_quad(hex) arrays. - * - * Called by the compute_support_point_weights_on_quad(hex) functions if the - * data is not yet hardcoded. - * - * For the definition of the support_point_weights_on_quad(hex) please - * refer to equation (8) of the `mapping' report. - */ - template - Table<2,double> - compute_laplace_vector(const unsigned int polynomial_degree) + namespace MappingQGeneric { - Table<2,double> lvs; - - Assert(lvs.n_rows()==0, ExcInternalError()); - Assert(dim==2 || dim==3, ExcNotImplemented()); - - // for degree==1, we shouldn't have to compute any support points, since all - // of them are on the vertices - Assert(polynomial_degree>1, ExcInternalError()); - - const unsigned int n_inner = Utilities::fixed_power(polynomial_degree-1); - const unsigned int n_outer = (dim==1) ? 2 : - ((dim==2) ? - 4+4*(polynomial_degree-1) : - 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1)); - - - // compute the shape gradients at the quadrature points on the unit cell - const QGauss quadrature(polynomial_degree+1); - const unsigned int n_q_points=quadrature.size(); - - typename MappingQGeneric::InternalData quadrature_data(polynomial_degree); - quadrature_data.shape_derivatives.resize(quadrature_data.n_shape_functions * - n_q_points); - quadrature_data.compute_shape_function_values(quadrature.get_points()); + namespace + { + /** + * Compute the support_point_weights_on_quad(hex) arrays. + * + * Called by the compute_support_point_weights_on_quad(hex) functions if the + * data is not yet hardcoded. + * + * For the definition of the support_point_weights_on_quad(hex) please + * refer to equation (8) of the `mapping' report. + */ + template + dealii::Table<2,double> + compute_laplace_vector(const unsigned int polynomial_degree) + { + dealii::Table<2,double> lvs; - // Compute the stiffness matrix of the inner dofs - FullMatrix S(n_inner); - for (unsigned int point=0; point1, ExcInternalError()); - // Compute the components of T to be the product of gradients of inner and - // outer shape functions. - FullMatrix T(n_inner, n_outer); - for (unsigned int point=0; point(polynomial_degree-1); + const unsigned int n_outer = (dim==1) ? 2 : + ((dim==2) ? + 4+4*(polynomial_degree-1) : + 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1)); - T(i,k) += res *(long double)quadrature.weight(point); - } - FullMatrix S_1(n_inner); - S_1.invert(S); + // compute the shape gradients at the quadrature points on the unit cell + const QGauss quadrature(polynomial_degree+1); + const unsigned int n_q_points=quadrature.size(); - FullMatrix S_1_T(n_inner, n_outer); + typename dealii::MappingQGeneric::InternalData quadrature_data(polynomial_degree); + quadrature_data.shape_derivatives.resize(quadrature_data.n_shape_functions * + n_q_points); + quadrature_data.compute_shape_function_values(quadrature.get_points()); - // S:=S_1*T - S_1.mmult(S_1_T,T); + // Compute the stiffness matrix of the inner dofs + FullMatrix S(n_inner); + for (unsigned int point=0; point T(n_inner, n_outer); + for (unsigned int point=0; pointMappingQ for dim= 2 and 3. - * - * For degree<4 this function sets the @p support_point_weights_on_quad to - * the hardcoded data. For degree>=4 and MappingQ<2> this vector is - * computed. - * - * For the definition of the @p support_point_weights_on_quad please refer to - * equation (8) of the `mapping' report. - */ - Table<2,double> - compute_support_point_weights_on_quad(const unsigned int polynomial_degree) - { - Table<2,double> loqvs; + FullMatrix S_1(n_inner); + S_1.invert(S); - // we are asked to compute weights for interior support points, but - // there are no interior points if degree==1 - if (polynomial_degree == 1) - return loqvs; + FullMatrix S_1_T(n_inner, n_outer); - const unsigned int n_inner_2d=(polynomial_degree-1)*(polynomial_degree-1); - const unsigned int n_outer_2d=4+4*(polynomial_degree-1); + // S:=S_1*T + S_1.mmult(S_1_T,T); - // first check whether we have precomputed the values for some polynomial - // degree; the sizes of arrays is n_inner_2d*n_outer_2d - if (polynomial_degree == 2) - { - // (checked these values against the output of compute_laplace_vector - // again, and found they're indeed right -- just in case someone wonders - // where they come from -- WB) - static const double loqv2[1*8] - = {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.}; - Assert (sizeof(loqv2)/sizeof(loqv2[0]) == - n_inner_2d * n_outer_2d, - ExcInternalError()); + // Resize and initialize the lvs + lvs.reinit (n_inner, n_outer); + for (unsigned int i=0; i(polynomial_degree); + return lvs; } - // the sum of weights of the points at the outer rim should be one. check - // this - for (unsigned int unit_point=0; unit_pointMappingQ<3>. - * - * For degree==2 this function sets the @p support_point_weights_on_hex to - * the hardcoded data. For degree>2 this vector is computed. - * - * For the definition of the @p support_point_weights_on_hex please refer to - * equation (8) of the `mapping' report. - */ - Table<2,double> - compute_support_point_weights_on_hex(const unsigned int polynomial_degree) - { - Table<2,double> lohvs; - - // we are asked to compute weights for interior support points, but - // there are no interior points if degree==1 - if (polynomial_degree == 1) - return lohvs; - const unsigned int n_inner = Utilities::fixed_power<3>(polynomial_degree-1); - const unsigned int n_outer = 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1); - // first check whether we have precomputed the values for some polynomial - // degree; the sizes of arrays is n_inner_2d*n_outer_2d - if (polynomial_degree == 2) + /** + * This function is needed by the constructor of + * MappingQ for dim= 2 and 3. + * + * For degree<4 this function sets the @p support_point_weights_on_quad to + * the hardcoded data. For degree>=4 and MappingQ<2> this vector is + * computed. + * + * For the definition of the @p support_point_weights_on_quad please refer to + * equation (8) of the `mapping' report. + */ + dealii::Table<2,double> + compute_support_point_weights_on_quad(const unsigned int polynomial_degree) { - static const double lohv2[26] - = {1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., - 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., - 7/192., 7/192., 7/192., 7/192., - 1/12., 1/12., 1/12., 1/12., 1/12., 1/12. - }; - - // copy and return - lohvs.reinit(n_inner, n_outer); - for (unsigned int unit_point=0; unit_point(polynomial_degree); - } - - // the sum of weights of the points at the outer rim should be one. check - // this - for (unsigned int unit_point=0; unit_point loqvs; - /** - * This function collects the output of - * compute_support_point_weights_on_{quad,hex} in a single data structure. - */ - std::vector > - compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree, - const unsigned int dim) - { - Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); - std::vector > output(dim); - if (polynomial_degree <= 1) - return output; - - // fill the 1D interior weights - QGaussLobatto<1> quadrature(polynomial_degree+1); - output[0].reinit(polynomial_degree-1, GeometryInfo<1>::vertices_per_cell); - for (unsigned int q=0; q::vertices_per_cell; ++i) - output[0](q,i) = GeometryInfo<1>::d_linear_shape_function(quadrature.point(q+1), - i); - - if (dim > 1) - output[1] = compute_support_point_weights_on_quad(polynomial_degree); - - if (dim > 2) - output[2] = compute_support_point_weights_on_hex(polynomial_degree); - - return output; - } - - /** - * Collects all interior points for the various dimensions. - */ - template - Table<2,double> - compute_support_point_weights_cell(const unsigned int polynomial_degree) - { - Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); - if (polynomial_degree <= 1) - return Table<2,double>(); + // we are asked to compute weights for interior support points, but + // there are no interior points if degree==1 + if (polynomial_degree == 1) + return loqvs; - QGaussLobatto quadrature(polynomial_degree+1); - std::vector h2l(quadrature.size()); - FETools::hierarchic_to_lexicographic_numbering(polynomial_degree, h2l); + const unsigned int n_inner_2d=(polynomial_degree-1)*(polynomial_degree-1); + const unsigned int n_outer_2d=4+4*(polynomial_degree-1); - Table<2,double> output(quadrature.size() - GeometryInfo::vertices_per_cell, - GeometryInfo::vertices_per_cell); - for (unsigned int q=0; q::vertices_per_cell; ++i) - output(q,i) = GeometryInfo::d_linear_shape_function(quadrature.point(h2l[q+GeometryInfo::vertices_per_cell]), - i); + // first check whether we have precomputed the values for some polynomial + // degree; the sizes of arrays is n_inner_2d*n_outer_2d + if (polynomial_degree == 2) + { + // (checked these values against the output of compute_laplace_vector + // again, and found they're indeed right -- just in case someone wonders + // where they come from -- WB) + static const double loqv2[1*8] + = {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.}; + Assert (sizeof(loqv2)/sizeof(loqv2[0]) == + n_inner_2d * n_outer_2d, + ExcInternalError()); + + // copy and return + loqvs.reinit(n_inner_2d, n_outer_2d); + for (unsigned int unit_point=0; unit_point(polynomial_degree); + } - return output; - } -} + // the sum of weights of the points at the outer rim should be one. check + // this + for (unsigned int unit_point=0; unit_point -MappingQGeneric::MappingQGeneric (const unsigned int p) - : - polynomial_degree(p), - line_support_points(this->polynomial_degree+1), - fe_q(dim == 3 ? new FE_Q(this->polynomial_degree) : nullptr), - support_point_weights_perimeter_to_interior (compute_support_point_weights_perimeter_to_interior(this->polynomial_degree, dim)), - support_point_weights_cell (compute_support_point_weights_cell(this->polynomial_degree)) -{ - Assert (p >= 1, ExcMessage ("It only makes sense to create polynomial mappings " - "with a polynomial degree greater or equal to one.")); -} + /** + * This function is needed by the constructor of MappingQ<3>. + * + * For degree==2 this function sets the @p support_point_weights_on_hex to + * the hardcoded data. For degree>2 this vector is computed. + * + * For the definition of the @p support_point_weights_on_hex please refer to + * equation (8) of the `mapping' report. + */ + dealii::Table<2,double> + compute_support_point_weights_on_hex(const unsigned int polynomial_degree) + { + dealii::Table<2,double> lohvs; + // we are asked to compute weights for interior support points, but + // there are no interior points if degree==1 + if (polynomial_degree == 1) + return lohvs; + const unsigned int n_inner = Utilities::fixed_power<3>(polynomial_degree-1); + const unsigned int n_outer = 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1); -template -MappingQGeneric::MappingQGeneric (const MappingQGeneric &mapping) - : - polynomial_degree(mapping.polynomial_degree), - line_support_points(mapping.line_support_points), - fe_q(dim == 3 ? new FE_Q(*mapping.fe_q) : nullptr), - support_point_weights_perimeter_to_interior (mapping.support_point_weights_perimeter_to_interior), - support_point_weights_cell (mapping.support_point_weights_cell) -{} + // first check whether we have precomputed the values for some polynomial + // degree; the sizes of arrays is n_inner_2d*n_outer_2d + if (polynomial_degree == 2) + { + static const double lohv2[26] + = {1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., + 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., + 7/192., 7/192., 7/192., 7/192., + 1/12., 1/12., 1/12., 1/12., 1/12., 1/12. + }; + + // copy and return + lohvs.reinit(n_inner, n_outer); + for (unsigned int unit_point=0; unit_point(polynomial_degree); + } + // the sum of weights of the points at the outer rim should be one. check + // this + for (unsigned int unit_point=0; unit_point -Mapping * -MappingQGeneric::clone () const -{ - return new MappingQGeneric(*this); -} + /** + * This function collects the output of + * compute_support_point_weights_on_{quad,hex} in a single data structure. + */ + std::vector > + compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree, + const unsigned int dim) + { + Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); + std::vector > output(dim); + if (polynomial_degree <= 1) + return output; + + // fill the 1D interior weights + QGaussLobatto<1> quadrature(polynomial_degree+1); + output[0].reinit(polynomial_degree-1, GeometryInfo<1>::vertices_per_cell); + for (unsigned int q=0; q::vertices_per_cell; ++i) + output[0](q,i) = GeometryInfo<1>::d_linear_shape_function(quadrature.point(q+1), + i); + + if (dim > 1) + output[1] = compute_support_point_weights_on_quad(polynomial_degree); + + if (dim > 2) + output[2] = compute_support_point_weights_on_hex(polynomial_degree); + + return output; + } + /** + * Collects all interior points for the various dimensions. + */ + template + dealii::Table<2,double> + compute_support_point_weights_cell(const unsigned int polynomial_degree) + { + Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); + if (polynomial_degree <= 1) + return dealii::Table<2,double>(); + + QGaussLobatto quadrature(polynomial_degree+1); + std::vector h2l(quadrature.size()); + FETools::hierarchic_to_lexicographic_numbering(polynomial_degree, h2l); + + dealii::Table<2,double> output(quadrature.size() - GeometryInfo::vertices_per_cell, + GeometryInfo::vertices_per_cell); + for (unsigned int q=0; q::vertices_per_cell; ++i) + output(q,i) = GeometryInfo::d_linear_shape_function(quadrature.point(h2l[q+GeometryInfo::vertices_per_cell]), + i); + + return output; + } -template -unsigned int -MappingQGeneric::get_degree() const -{ - return polynomial_degree; -} + /** + * Using the relative weights of the shape functions evaluated at + * one point on the reference cell (and stored in data.shape_values + * and accessed via data.shape(0,i)) and the locations of mapping + * support points (stored in data.mapping_support_points), compute + * the mapped location of that point in real space. + */ + template + Point + compute_mapped_location_of_point + (const typename dealii::MappingQGeneric::InternalData &data) + { + AssertDimension (data.shape_values.size(), + data.mapping_support_points.size()); + // use now the InternalData to compute the point in real space. + Point p_real; + for (unsigned int i=0; i -Point -MappingQGeneric:: -transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, - const Point &p) const -{ - // set up the polynomial space - const TensorProductPolynomials - tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points())); - Assert (tensor_pols.n() == Utilities::fixed_power(polynomial_degree+1), - ExcInternalError()); + return p_real; + } - // then also construct the mapping from lexicographic to the Qp shape function numbering - const std::vector - renumber (FETools:: - lexicographic_to_hierarchic_numbering ( - FiniteElementData (get_dpo_vector(polynomial_degree), 1, - polynomial_degree))); - const std::vector > support_points - = this->compute_mapping_support_points(cell); - Point mapped_point; - for (unsigned int i=0; i + Point + do_transform_real_to_unit_cell_internal + (const typename dealii::Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + typename dealii::MappingQGeneric::InternalData &mdata) + { + const unsigned int spacedim = dim; - return mapped_point; -} + const unsigned int n_shapes=mdata.shape_values.size(); + (void)n_shapes; + Assert(n_shapes!=0, ExcInternalError()); + AssertDimension (mdata.shape_derivatives.size(), n_shapes); + std::vector > &points=mdata.mapping_support_points; + AssertDimension (points.size(), n_shapes); -// In the code below, GCC tries to instantiate MappingQGeneric<3,4> when -// seeing which of the overloaded versions of -// do_transform_real_to_unit_cell_internal() to call. This leads to bad -// error messages and, generally, nothing very good. Avoid this by ensuring -// that this class exists, but does not have an inner InternalData -// type, thereby ruling out the codim-1 version of the function -// below when doing overload resolution. -template <> -class MappingQGeneric<3,4> -{}; -namespace -{ - /** - * Using the relative weights of the shape functions evaluated at - * one point on the reference cell (and stored in data.shape_values - * and accessed via data.shape(0,i)) and the locations of mapping - * support points (stored in data.mapping_support_points), compute - * the mapped location of that point in real space. - */ - template - Point - compute_mapped_location_of_point (const typename MappingQGeneric::InternalData &data) - { - AssertDimension (data.shape_values.size(), - data.mapping_support_points.size()); + // Newton iteration to solve + // f(x)=p(x)-p=0 + // where we are looking for 'x' and p(x) is the forward transformation + // from unit to real cell. We solve this using a Newton iteration + // x_{n+1}=x_n-[f'(x)]^{-1}f(x) + // The start value is set to be the linear approximation to the cell - // use now the InternalData to compute the point in real space. - Point p_real; - for (unsigned int i=0; i p_unit = initial_p_unit; + mdata.compute_shape_function_values(std::vector > (1, p_unit)); - /** - * Implementation of transform_real_to_unit_cell for dim==spacedim - */ - template - Point - do_transform_real_to_unit_cell_internal - (const typename Triangulation::cell_iterator &cell, - const Point &p, - const Point &initial_p_unit, - typename MappingQGeneric::InternalData &mdata) - { - const unsigned int spacedim = dim; - - const unsigned int n_shapes=mdata.shape_values.size(); - (void)n_shapes; - Assert(n_shapes!=0, ExcInternalError()); - AssertDimension (mdata.shape_derivatives.size(), n_shapes); - - std::vector > &points=mdata.mapping_support_points; - AssertDimension (points.size(), n_shapes); - - - // Newton iteration to solve - // f(x)=p(x)-p=0 - // where we are looking for 'x' and p(x) is the forward transformation - // from unit to real cell. We solve this using a Newton iteration - // x_{n+1}=x_n-[f'(x)]^{-1}f(x) - // The start value is set to be the linear approximation to the cell - - // The shape values and derivatives of the mapping at this point are - // previously computed. - - Point p_unit = initial_p_unit; - - mdata.compute_shape_function_values(std::vector > (1, p_unit)); - - Point p_real = compute_mapped_location_of_point(mdata); - Tensor<1,spacedim> f = p_real-p; - - // early out if we already have our point - if (f.norm_square() < 1e-24 * cell->diameter() * cell->diameter()) - return p_unit; - - // we need to compare the position of the computed p(x) against the given - // point 'p'. We will terminate the iteration and return 'x' if they are - // less than eps apart. The question is how to choose eps -- or, put maybe - // more generally: in which norm we want these 'p' and 'p(x)' to be eps - // apart. - // - // the question is difficult since we may have to deal with very elongated - // cells where we may achieve 1e-12*h for the distance of these two points - // in the 'long' direction, but achieving this tolerance in the 'short' - // direction of the cell may not be possible - // - // what we do instead is then to terminate iterations if - // \| p(x) - p \|_A < eps - // where the A-norm is somehow induced by the transformation of the cell. - // in particular, we want to measure distances relative to the sizes of - // the cell in its principal directions. - // - // to define what exactly A should be, note that to first order we have - // the following (assuming that x* is the solution of the problem, i.e., - // p(x*)=p): - // p(x) - p = p(x) - p(x*) - // = -grad p(x) * (x*-x) + higher order terms - // This suggest to measure with a norm that corresponds to - // A = {[grad p(x]^T [grad p(x)]}^{-1} - // because then - // \| p(x) - p \|_A \approx \| x - x* \| - // Consequently, we will try to enforce that - // \| p(x) - p \|_A = \| f \| <= eps - // - // Note that using this norm is a bit dangerous since the norm changes - // in every iteration (A isn't fixed by depends on xk). However, if the - // cell is not too deformed (it may be stretched, but not twisted) then - // the mapping is almost linear and A is indeed constant or nearly so. - const double eps = 1.e-11; - const unsigned int newton_iteration_limit = 20; - - unsigned int newton_iteration = 0; - double last_f_weighted_norm; - do - { + Point p_real = compute_mapped_location_of_point(mdata); + Tensor<1,spacedim> f = p_real-p; + + // early out if we already have our point + if (f.norm_square() < 1e-24 * cell->diameter() * cell->diameter()) + return p_unit; + + // we need to compare the position of the computed p(x) against the given + // point 'p'. We will terminate the iteration and return 'x' if they are + // less than eps apart. The question is how to choose eps -- or, put maybe + // more generally: in which norm we want these 'p' and 'p(x)' to be eps + // apart. + // + // the question is difficult since we may have to deal with very elongated + // cells where we may achieve 1e-12*h for the distance of these two points + // in the 'long' direction, but achieving this tolerance in the 'short' + // direction of the cell may not be possible + // + // what we do instead is then to terminate iterations if + // \| p(x) - p \|_A < eps + // where the A-norm is somehow induced by the transformation of the cell. + // in particular, we want to measure distances relative to the sizes of + // the cell in its principal directions. + // + // to define what exactly A should be, note that to first order we have + // the following (assuming that x* is the solution of the problem, i.e., + // p(x*)=p): + // p(x) - p = p(x) - p(x*) + // = -grad p(x) * (x*-x) + higher order terms + // This suggest to measure with a norm that corresponds to + // A = {[grad p(x]^T [grad p(x)]}^{-1} + // because then + // \| p(x) - p \|_A \approx \| x - x* \| + // Consequently, we will try to enforce that + // \| p(x) - p \|_A = \| f \| <= eps + // + // Note that using this norm is a bit dangerous since the norm changes + // in every iteration (A isn't fixed by depends on xk). However, if the + // cell is not too deformed (it may be stretched, but not twisted) then + // the mapping is almost linear and A is indeed constant or nearly so. + const double eps = 1.e-11; + const unsigned int newton_iteration_limit = 20; + + unsigned int newton_iteration = 0; + double last_f_weighted_norm; + do + { #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << "Newton iteration " << newton_iteration << std::endl; + std::cout << "Newton iteration " << newton_iteration << std::endl; #endif - // f'(x) - Tensor<2,spacedim> df; - for (unsigned int k=0; k &grad_transform=mdata.derivative(0,k); - const Point &point=points[k]; + // f'(x) + Tensor<2,spacedim> df; + for (unsigned int k=0; k &grad_transform=mdata.derivative(0,k); + const Point &point=points[k]; - for (unsigned int i=0; i 0, - (typename Mapping::ExcTransformationFailed())); - Tensor<2,spacedim> df_inverse = invert(df); - const Tensor<1,spacedim> delta = df_inverse * static_cast&>(f); + // Solve [f'(x)]d=f(x) + AssertThrow(determinant(df) > 0, + (typename Mapping::ExcTransformationFailed())); + Tensor<2,spacedim> df_inverse = invert(df); + const Tensor<1,spacedim> delta = df_inverse * static_cast&>(f); #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << " delta=" << delta << std::endl; + std::cout << " delta=" << delta << std::endl; #endif - // do a line search - double step_length = 1; - do - { - // update of p_unit. The spacedim-th component of transformed point - // is simply ignored in codimension one case. When this component is - // not zero, then we are projecting the point to the surface or - // curve identified by the cell. - Point p_unit_trial = p_unit; - for (unsigned int i=0; i > (1, p_unit_trial)); - - // f(x) - Point p_real_trial = compute_mapped_location_of_point(mdata); - const Tensor<1,spacedim> f_trial = p_real_trial-p; + // do a line search + double step_length = 1; + do + { + // update of p_unit. The spacedim-th component of transformed point + // is simply ignored in codimension one case. When this component is + // not zero, then we are projecting the point to the surface or + // curve identified by the cell. + Point p_unit_trial = p_unit; + for (unsigned int i=0; i > (1, p_unit_trial)); + + // f(x) + Point p_real_trial = internal::MappingQGeneric::compute_mapped_location_of_point(mdata); + const Tensor<1,spacedim> f_trial = p_real_trial-p; #ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << " step_length=" << step_length << std::endl - << " ||f || =" << f.norm() << std::endl - << " ||f*|| =" << f_trial.norm() << std::endl - << " ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl; + std::cout << " step_length=" << step_length << std::endl + << " ||f || =" << f.norm() << std::endl + << " ||f*|| =" << f_trial.norm() << std::endl + << " ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl; #endif - // see if we are making progress with the current step length - // and if not, reduce it by a factor of two and try again - // - // strictly speaking, we should probably use the same norm as we use - // for the outer algorithm. in practice, line search is just a - // crutch to find a "reasonable" step length, and so using the l2 - // norm is probably just fine - if (f_trial.norm() < f.norm()) - { - p_real = p_real_trial; - p_unit = p_unit_trial; - f = f_trial; - break; + // see if we are making progress with the current step length + // and if not, reduce it by a factor of two and try again + // + // strictly speaking, we should probably use the same norm as we use + // for the outer algorithm. in practice, line search is just a + // crutch to find a "reasonable" step length, and so using the l2 + // norm is probably just fine + if (f_trial.norm() < f.norm()) + { + p_real = p_real_trial; + p_unit = p_unit_trial; + f = f_trial; + break; + } + else if (step_length > 0.05) + step_length /= 2; + else + AssertThrow (false, + (typename Mapping::ExcTransformationFailed())); } - else if (step_length > 0.05) - step_length /= 2; - else + while (true); + + ++newton_iteration; + if (newton_iteration > newton_iteration_limit) AssertThrow (false, (typename Mapping::ExcTransformationFailed())); + last_f_weighted_norm = (df_inverse * f).norm(); } - while (true); - - ++newton_iteration; - if (newton_iteration > newton_iteration_limit) - AssertThrow (false, - (typename Mapping::ExcTransformationFailed())); - last_f_weighted_norm = (df_inverse * f).norm(); - } - while (last_f_weighted_norm > eps); - - return p_unit; - } - - - - /** - * Implementation of transform_real_to_unit_cell for dim==spacedim-1 - */ - template - Point - do_transform_real_to_unit_cell_internal_codim1 - (const typename Triangulation::cell_iterator &cell, - const Point &p, - const Point &initial_p_unit, - typename MappingQGeneric::InternalData &mdata) - { - const unsigned int spacedim = dim+1; - - const unsigned int n_shapes=mdata.shape_values.size(); - (void)n_shapes; - Assert(n_shapes!=0, ExcInternalError()); - Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError()); - Assert(mdata.shape_second_derivatives.size()==n_shapes, ExcInternalError()); + while (last_f_weighted_norm > eps); - std::vector > &points=mdata.mapping_support_points; - Assert(points.size()==n_shapes, ExcInternalError()); - - Point p_minus_F; - - Tensor<1,spacedim> DF[dim]; - Tensor<1,spacedim> D2F[dim][dim]; - - Point p_unit = initial_p_unit; - Point f; - Tensor<2,dim> df; - - // Evaluate first and second derivatives - mdata.compute_shape_function_values(std::vector > (1, p_unit)); - - for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); - const Tensor<2,dim> &hessian_k = mdata.second_derivative(0,k); - const Point &point_k = points[k]; - - for (unsigned int j=0; j(mdata); - - for (unsigned int j=0; j + Point + do_transform_real_to_unit_cell_internal_codim1 + (const typename dealii::Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + typename dealii::MappingQGeneric::InternalData &mdata) { - f[j] = DF[j] * p_minus_F; - for (unsigned int l=0; ldiameter(); - const unsigned int loop_limit = 10; + std::vector > &points=mdata.mapping_support_points; + Assert(points.size()==n_shapes, ExcInternalError()); - unsigned int loop=0; + Point p_minus_F; - while (f.norm()>eps && loop++ d = invert(df) * static_cast&>(f); - p_unit -= d; + Tensor<1,spacedim> DF[dim]; + Tensor<1,spacedim> D2F[dim][dim]; - for (unsigned int j=0; j p_unit = initial_p_unit; + Point f; + Tensor<2,dim> df; + // Evaluate first and second derivatives mdata.compute_shape_function_values(std::vector > (1, p_unit)); for (unsigned int k=0; k(mdata); + + for (unsigned int j=0; jdiameter(); + const unsigned int loop_limit = 10; - // Here we check that in the last execution of while the first - // condition was already wrong, meaning the residual was below - // eps. Only if the first condition failed, loop will have been - // increased and tested, and thus have reached the limit. - AssertThrow (loop::ExcTransformationFailed())); + unsigned int loop=0; - return p_unit; - } + while (f.norm()>eps && loop++ d = invert(df) * static_cast&>(f); + p_unit -= d; + for (unsigned int j=0; j > (1, p_unit)); + for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); + const Tensor<2,dim> &hessian_k = mdata.second_derivative(0,k); + const Point &point_k = points[k]; + for (unsigned int j=0; j(mdata); -template -Point -MappingQGeneric:: -transform_real_to_unit_cell_internal -(const typename Triangulation::cell_iterator &, - const Point &, - const Point &) const -{ - // default implementation (should never be called) - Assert(false, ExcInternalError()); - return Point(); -} + for (unsigned int j=0; j -Point<1> -MappingQGeneric<1,1>:: -transform_real_to_unit_cell_internal -(const Triangulation<1,1>::cell_iterator &cell, - const Point<1> &p, - const Point<1> &initial_p_unit) const -{ - const int dim = 1; - const int spacedim = 1; + } - const Quadrature point_quadrature(initial_p_unit); - UpdateFlags update_flags = update_quadrature_points | update_jacobians; - if (spacedim>dim) - update_flags |= update_jacobian_grads; - std::unique_ptr mdata (get_data(update_flags, - point_quadrature)); + // Here we check that in the last execution of while the first + // condition was already wrong, meaning the residual was below + // eps. Only if the first condition failed, loop will have been + // increased and tested, and thus have reached the limit. + AssertThrow (loop::ExcTransformationFailed())); - mdata->mapping_support_points = this->compute_mapping_support_points (cell); + return p_unit; + } + + /** + * 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 + void + maybe_compute_q_points + (const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector > &quadrature_points) + { + const UpdateFlags update_flags = data.update_each; + + if (update_flags & update_quadrature_points) + { + for (unsigned int point=0; point result = (shape[0] * + data.mapping_support_points[0]); + for (unsigned int k=1; k + void + maybe_update_Jacobians + (const CellSimilarity::Similarity cell_similarity, + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data) + { + const UpdateFlags update_flags = data.update_each; + + 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 *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 + void + maybe_update_jacobian_grads + (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector > &jacobian_grads) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_grads) + { + const unsigned int n_q_points = jacobian_grads.size(); + + if (cell_similarity != CellSimilarity::translation) + { + for (unsigned int point=0; point *second = + &data.second_derivative(point+data_set, 0); + double result [spacedim][dim][dim]; + for (unsigned int i=0; i + void + maybe_update_jacobian_pushed_forward_grads + (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector > &jacobian_pushed_forward_grads) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_grads) + { + const unsigned int n_q_points = jacobian_pushed_forward_grads.size(); + + if (cell_similarity != CellSimilarity::translation) + { + double tmp[spacedim][spacedim][spacedim]; + for (unsigned int point=0; point *second = + &data.second_derivative(point+data_set, 0); + double result [spacedim][dim][dim]; + for (unsigned int i=0; i + void + maybe_update_jacobian_2nd_derivatives + (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector > &jacobian_2nd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_2nd_derivatives) + { + const unsigned int n_q_points = jacobian_2nd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + for (unsigned int point=0; point *third = + &data.third_derivative(point+data_set, 0); + double result [spacedim][dim][dim][dim]; + for (unsigned int i=0; i + void + maybe_update_jacobian_pushed_forward_2nd_derivatives + (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector > &jacobian_pushed_forward_2nd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) + { + const unsigned int n_q_points = jacobian_pushed_forward_2nd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + double tmp[spacedim][spacedim][spacedim][spacedim]; + for (unsigned int point=0; point *third = + &data.third_derivative(point+data_set, 0); + double result [spacedim][dim][dim][dim]; + for (unsigned int i=0; i + void + maybe_update_jacobian_3rd_derivatives + (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector > &jacobian_3rd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_3rd_derivatives) + { + const unsigned int n_q_points = jacobian_3rd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + for (unsigned int point=0; point *fourth = + &data.fourth_derivative(point+data_set, 0); + double result [spacedim][dim][dim][dim][dim]; + for (unsigned int i=0; i + void + maybe_update_jacobian_pushed_forward_3rd_derivatives + (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector > &jacobian_pushed_forward_3rd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) + { + const unsigned int n_q_points = jacobian_pushed_forward_3rd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + double tmp[spacedim][spacedim][spacedim][spacedim][spacedim]; + for (unsigned int point=0; point *fourth = + &data.fourth_derivative(point+data_set, 0); + double result [spacedim][dim][dim][dim][dim]; + for (unsigned int i=0; i +MappingQGeneric::MappingQGeneric (const unsigned int p) + : + polynomial_degree(p), + line_support_points(this->polynomial_degree+1), + fe_q(dim == 3 ? new FE_Q(this->polynomial_degree) : nullptr), + support_point_weights_perimeter_to_interior (internal::MappingQGeneric::compute_support_point_weights_perimeter_to_interior(this->polynomial_degree, dim)), + support_point_weights_cell (internal::MappingQGeneric::compute_support_point_weights_cell(this->polynomial_degree)) +{ + Assert (p >= 1, ExcMessage ("It only makes sense to create polynomial mappings " + "with a polynomial degree greater or equal to one.")); +} + + + +template +MappingQGeneric::MappingQGeneric (const MappingQGeneric &mapping) + : + polynomial_degree(mapping.polynomial_degree), + line_support_points(mapping.line_support_points), + fe_q(dim == 3 ? new FE_Q(*mapping.fe_q) : nullptr), + support_point_weights_perimeter_to_interior (mapping.support_point_weights_perimeter_to_interior), + support_point_weights_cell (mapping.support_point_weights_cell) +{} + + + + +template +Mapping * +MappingQGeneric::clone () const +{ + return new MappingQGeneric(*this); +} + + + + +template +unsigned int +MappingQGeneric::get_degree() const +{ + return polynomial_degree; +} + + + +template +Point +MappingQGeneric:: +transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, + const Point &p) const +{ + // set up the polynomial space + const TensorProductPolynomials + tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points())); + Assert (tensor_pols.n() == Utilities::fixed_power(polynomial_degree+1), + ExcInternalError()); + + // then also construct the mapping from lexicographic to the Qp shape function numbering + const std::vector + renumber (FETools:: + lexicographic_to_hierarchic_numbering + (FiniteElementData (internal::MappingQGeneric::get_dpo_vector + (polynomial_degree), 1, polynomial_degree))); + + const std::vector > support_points + = this->compute_mapping_support_points(cell); + + Point mapped_point; + for (unsigned int i=0; i when +// seeing which of the overloaded versions of +// do_transform_real_to_unit_cell_internal() to call. This leads to bad +// error messages and, generally, nothing very good. Avoid this by ensuring +// that this class exists, but does not have an inner InternalData +// type, thereby ruling out the codim-1 version of the function +// below when doing overload resolution. +template <> +class MappingQGeneric<3,4> +{}; + + + +// visual studio freaks out when trying to determine if +// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good +// candidate. So instead of letting the compiler pick the correct overload, we +// use template specialization to make sure we pick up the right function to +// call: + +template +Point +MappingQGeneric:: +transform_real_to_unit_cell_internal +(const typename Triangulation::cell_iterator &, + const Point &, + const Point &) const +{ + // default implementation (should never be called) + Assert(false, ExcInternalError()); + return Point(); +} + +template <> +Point<1> +MappingQGeneric<1,1>:: +transform_real_to_unit_cell_internal +(const Triangulation<1,1>::cell_iterator &cell, + const Point<1> &p, + const Point<1> &initial_p_unit) const +{ + const int dim = 1; + const int spacedim = 1; + + const Quadrature point_quadrature(initial_p_unit); + + UpdateFlags update_flags = update_quadrature_points | update_jacobians; + if (spacedim>dim) + update_flags |= update_jacobian_grads; + std::unique_ptr mdata (get_data(update_flags, + point_quadrature)); + + mdata->mapping_support_points = this->compute_mapping_support_points (cell); // dispatch to the various specializations for spacedim=dim, // spacedim=dim+1, etc - return do_transform_real_to_unit_cell_internal<1>(cell, p, initial_p_unit, *mdata); + return internal::MappingQGeneric::do_transform_real_to_unit_cell_internal<1>(cell, p, initial_p_unit, *mdata); } template <> @@ -1570,7 +2124,7 @@ transform_real_to_unit_cell_internal // dispatch to the various specializations for spacedim=dim, // spacedim=dim+1, etc - return do_transform_real_to_unit_cell_internal<2>(cell, p, initial_p_unit, *mdata); + return internal::MappingQGeneric::do_transform_real_to_unit_cell_internal<2>(cell, p, initial_p_unit, *mdata); } template <> @@ -1596,7 +2150,7 @@ transform_real_to_unit_cell_internal // dispatch to the various specializations for spacedim=dim, // spacedim=dim+1, etc - return do_transform_real_to_unit_cell_internal<3>(cell, p, initial_p_unit, *mdata); + return internal::MappingQGeneric::do_transform_real_to_unit_cell_internal<3>(cell, p, initial_p_unit, *mdata); } template <> @@ -1622,7 +2176,7 @@ transform_real_to_unit_cell_internal // dispatch to the various specializations for spacedim=dim, // spacedim=dim+1, etc - return do_transform_real_to_unit_cell_internal_codim1<1>(cell, p, initial_p_unit, *mdata); + return internal::MappingQGeneric::do_transform_real_to_unit_cell_internal_codim1<1>(cell, p, initial_p_unit, *mdata); } template <> @@ -1648,7 +2202,7 @@ transform_real_to_unit_cell_internal // dispatch to the various specializations for spacedim=dim, // spacedim=dim+1, etc - return do_transform_real_to_unit_cell_internal_codim1<2>(cell, p, initial_p_unit, *mdata); + return internal::MappingQGeneric::do_transform_real_to_unit_cell_internal_codim1<2>(cell, p, initial_p_unit, *mdata); } template <> @@ -1866,584 +2420,39 @@ MappingQGeneric::get_data (const UpdateFlags update_flags, data->initialize (this->requires_update_flags(update_flags), q, q.size()); return data; -} - - - -template -typename MappingQGeneric::InternalData * -MappingQGeneric::get_face_data (const UpdateFlags update_flags, - const Quadrature &quadrature) const -{ - InternalData *data = new InternalData(polynomial_degree); - data->initialize_face (this->requires_update_flags(update_flags), - QProjector::project_to_all_faces(quadrature), - quadrature.size()); - - return data; -} - - - -template -typename MappingQGeneric::InternalData * -MappingQGeneric::get_subface_data (const UpdateFlags update_flags, - const Quadrature& quadrature) const -{ - InternalData *data = new InternalData(polynomial_degree); - data->initialize_face (this->requires_update_flags(update_flags), - QProjector::project_to_all_subfaces(quadrature), - quadrature.size()); - - return data; -} - - - -namespace internal -{ - 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 - void - maybe_compute_q_points (const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData &data, - std::vector > &quadrature_points) - { - const UpdateFlags update_flags = data.update_each; - - if (update_flags & update_quadrature_points) - { - for (unsigned int point=0; point result = (shape[0] * - data.mapping_support_points[0]); - for (unsigned int k=1; k - void - maybe_update_Jacobians (const CellSimilarity::Similarity cell_similarity, - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData &data) - { - const UpdateFlags update_flags = data.update_each; - - 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 *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 - void - maybe_update_jacobian_grads (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData &data, - std::vector > &jacobian_grads) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_grads) - { - const unsigned int n_q_points = jacobian_grads.size(); - - if (cell_similarity != CellSimilarity::translation) - { - for (unsigned int point=0; point *second = - &data.second_derivative(point+data_set, 0); - double result [spacedim][dim][dim]; - for (unsigned int i=0; i - void - maybe_update_jacobian_pushed_forward_grads (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData &data, - std::vector > &jacobian_pushed_forward_grads) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_grads) - { - const unsigned int n_q_points = jacobian_pushed_forward_grads.size(); - - if (cell_similarity != CellSimilarity::translation) - { - double tmp[spacedim][spacedim][spacedim]; - for (unsigned int point=0; point *second = - &data.second_derivative(point+data_set, 0); - double result [spacedim][dim][dim]; - for (unsigned int i=0; i - void - maybe_update_jacobian_2nd_derivatives (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData &data, - std::vector > &jacobian_2nd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_2nd_derivatives) - { - const unsigned int n_q_points = jacobian_2nd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - for (unsigned int point=0; point *third = - &data.third_derivative(point+data_set, 0); - double result [spacedim][dim][dim][dim]; - for (unsigned int i=0; i - void - maybe_update_jacobian_pushed_forward_2nd_derivatives (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData &data, - std::vector > &jacobian_pushed_forward_2nd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) - { - const unsigned int n_q_points = jacobian_pushed_forward_2nd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - double tmp[spacedim][spacedim][spacedim][spacedim]; - for (unsigned int point=0; point *third = - &data.third_derivative(point+data_set, 0); - double result [spacedim][dim][dim][dim]; - for (unsigned int i=0; i - void - maybe_update_jacobian_3rd_derivatives (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData &data, - std::vector > &jacobian_3rd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_3rd_derivatives) - { - const unsigned int n_q_points = jacobian_3rd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - for (unsigned int point=0; point *fourth = - &data.fourth_derivative(point+data_set, 0); - double result [spacedim][dim][dim][dim][dim]; - for (unsigned int i=0; i - void - maybe_update_jacobian_pushed_forward_3rd_derivatives (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingQGeneric::InternalData &data, - std::vector > &jacobian_pushed_forward_3rd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) - { - const unsigned int n_q_points = jacobian_pushed_forward_3rd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - double tmp[spacedim][spacedim][spacedim][spacedim][spacedim]; - for (unsigned int point=0; point *fourth = - &data.fourth_derivative(point+data_set, 0); - double result [spacedim][dim][dim][dim][dim]; - for (unsigned int i=0; i +typename MappingQGeneric::InternalData * +MappingQGeneric::get_face_data (const UpdateFlags update_flags, + const Quadrature &quadrature) const +{ + InternalData *data = new InternalData(polynomial_degree); + data->initialize_face (this->requires_update_flags(update_flags), + QProjector::project_to_all_faces(quadrature), + quadrature.size()); - // push-forward the n-coordinate - for (unsigned int i=0; i +typename MappingQGeneric::InternalData * +MappingQGeneric::get_subface_data (const UpdateFlags update_flags, + const Quadrature& quadrature) const +{ + InternalData *data = new InternalData(polynomial_degree); + data->initialize_face (this->requires_update_flags(update_flags), + QProjector::project_to_all_subfaces(quadrature), + quadrature.size()); + + return data; +} + + template CellSimilarity::Similarity @@ -2481,12 +2490,14 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, const CellSimilarity::Similarity computed_cell_similarity = (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none); - internal::maybe_compute_q_points (QProjector::DataSetDescriptor::cell (), - data, - output_data.quadrature_points); - internal::maybe_update_Jacobians (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data); + internal::MappingQGeneric::maybe_compute_q_points + (QProjector::DataSetDescriptor::cell (), + data, + output_data.quadrature_points); + internal::MappingQGeneric::maybe_update_Jacobians + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data); const UpdateFlags update_flags = data.update_each; const std::vector &weights=quadrature.get_weights(); @@ -2599,35 +2610,41 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, output_data.inverse_jacobians[point] = data.covariant[point].transpose(); } - internal::maybe_update_jacobian_grads (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_grads); - - internal::maybe_update_jacobian_pushed_forward_grads (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_pushed_forward_grads); - - internal::maybe_update_jacobian_2nd_derivatives (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_2nd_derivatives); - - internal::maybe_update_jacobian_pushed_forward_2nd_derivatives (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_pushed_forward_2nd_derivatives); - - internal::maybe_update_jacobian_3rd_derivatives (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_3rd_derivatives); - - internal::maybe_update_jacobian_pushed_forward_3rd_derivatives (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_pushed_forward_3rd_derivatives); + internal::MappingQGeneric::maybe_update_jacobian_grads + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_grads); + + internal::MappingQGeneric::maybe_update_jacobian_pushed_forward_grads + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_pushed_forward_grads); + + internal::MappingQGeneric::maybe_update_jacobian_2nd_derivatives + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_2nd_derivatives); + + internal::MappingQGeneric::maybe_update_jacobian_pushed_forward_2nd_derivatives + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_pushed_forward_2nd_derivatives); + + internal::MappingQGeneric::maybe_update_jacobian_3rd_derivatives + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_3rd_derivatives); + + internal::MappingQGeneric::maybe_update_jacobian_pushed_forward_3rd_derivatives + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_pushed_forward_3rd_derivatives); return computed_cell_similarity; } @@ -2639,214 +2656,217 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, namespace internal { - namespace + namespace MappingQGeneric { - /** - * 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 - void - maybe_compute_face_data (const dealii::MappingQGeneric &mapping, - const typename dealii::Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const unsigned int n_q_points, - const std::vector &weights, - const typename dealii::MappingQGeneric::InternalData &data, - internal::FEValues::MappingRelatedData &output_data) + namespace { - const UpdateFlags update_flags = data.update_each; + /** + * 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 + void + maybe_compute_face_data (const dealii::MappingQGeneric &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const unsigned int n_q_points, + const std::vector &weights, + const typename dealii::MappingQGeneric::InternalData &data, + internal::FEValues::MappingRelatedData &output_data) + { + const UpdateFlags update_flags = data.update_each; - if (update_flags & (update_boundary_forms | - update_normal_vectors | - update_jacobians | - update_JxW_values | - update_inverse_jacobians)) - { - if (update_flags & update_boundary_forms) - AssertDimension (output_data.boundary_forms.size(), n_q_points); - if (update_flags & update_normal_vectors) - AssertDimension (output_data.normal_vectors.size(), n_q_points); - if (update_flags & update_JxW_values) - AssertDimension (output_data.JxW_values.size(), n_q_points); - - Assert (data.aux.size()+1 >= dim, ExcInternalError()); - - // first compute some common data that is used for evaluating - // all of the flags below - - // map the unit tangentials to the real cell. checking for d!=dim-1 - // eliminates compiler warnings regarding unsigned int expressions < - // 0. - for (unsigned int d=0; d!=dim-1; ++d) - { - Assert (face_no+GeometryInfo::faces_per_cell*d < - data.unit_tangentials.size(), - ExcInternalError()); - Assert (data.aux[d].size() <= - data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d].size(), - ExcInternalError()); - - mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d]), - mapping_contravariant, - data, - make_array_view(data.aux[d])); - } + if (update_flags & (update_boundary_forms | + update_normal_vectors | + update_jacobians | + update_JxW_values | + update_inverse_jacobians)) + { + if (update_flags & update_boundary_forms) + AssertDimension (output_data.boundary_forms.size(), n_q_points); + if (update_flags & update_normal_vectors) + AssertDimension (output_data.normal_vectors.size(), n_q_points); + if (update_flags & update_JxW_values) + AssertDimension (output_data.JxW_values.size(), n_q_points); + + Assert (data.aux.size()+1 >= dim, ExcInternalError()); + + // first compute some common data that is used for evaluating + // all of the flags below + + // map the unit tangentials to the real cell. checking for d!=dim-1 + // eliminates compiler warnings regarding unsigned int expressions < + // 0. + for (unsigned int d=0; d!=dim-1; ++d) + { + Assert (face_no+GeometryInfo::faces_per_cell*d < + data.unit_tangentials.size(), + ExcInternalError()); + Assert (data.aux[d].size() <= + data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d].size(), + ExcInternalError()); + + mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d]), + mapping_contravariant, + data, + make_array_view(data.aux[d])); + } - if (update_flags & update_boundary_forms) - { - // if dim==spacedim, we can use the unit tangentials to compute the - // boundary form by simply taking the cross product - if (dim == spacedim) - { - for (unsigned int i=0; i DX_t = - data.contravariant[point].transpose(); - - Tensor<1, spacedim> cell_normal = - cross_product_3d(DX_t[0], DX_t[1]); - cell_normal /= cell_normal.norm(); - - // then compute the face normal from the face tangent - // and the cell normal: - output_data.boundary_forms[point] = - cross_product_3d(data.aux[0][point], cell_normal); - } - } - } - } + if (dim==2) + { + const DerivativeForm<1,spacedim,dim> DX_t = + data.contravariant[point].transpose(); - if (update_flags & update_JxW_values) - for (unsigned int i=0; i cell_normal = + cross_product_3d(DX_t[0], DX_t[1]); + cell_normal /= cell_normal.norm(); - if (subface_no != numbers::invalid_unsigned_int) - { - const double area_ratio = GeometryInfo::subface_ratio(cell->subface_case(face_no), - subface_no); - output_data.JxW_values[i] *= area_ratio; + // then compute the face normal from the face tangent + // and the cell normal: + output_data.boundary_forms[point] = + cross_product_3d(data.aux[0][point], cell_normal); + } + } } } - if (update_flags & update_normal_vectors) - for (unsigned int i=0; i(output_data.boundary_forms[i] / - output_data.boundary_forms[i].norm()); + if (update_flags & update_JxW_values) + for (unsigned int i=0; i::subface_ratio(cell->subface_case(face_no), + subface_no); + output_data.JxW_values[i] *= area_ratio; + } + } + + if (update_flags & update_normal_vectors) + for (unsigned int i=0; i(output_data.boundary_forms[i] / + output_data.boundary_forms[i].norm()); - if (update_flags & update_jacobians) - for (unsigned int point=0; point - void - do_fill_fe_face_values (const dealii::MappingQGeneric &mapping, - const typename dealii::Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const typename QProjector::DataSetDescriptor data_set, - const Quadrature &quadrature, - const typename dealii::MappingQGeneric::InternalData &data, - internal::FEValues::MappingRelatedData &output_data) - { - maybe_compute_q_points (data_set, - data, - output_data.quadrature_points); - maybe_update_Jacobians (CellSimilarity::none, - data_set, - data); - maybe_update_jacobian_grads (CellSimilarity::none, - data_set, - data, - output_data.jacobian_grads); - maybe_update_jacobian_pushed_forward_grads (CellSimilarity::none, - data_set, - data, - output_data.jacobian_pushed_forward_grads); - maybe_update_jacobian_2nd_derivatives (CellSimilarity::none, - data_set, - data, - output_data.jacobian_2nd_derivatives); - maybe_update_jacobian_pushed_forward_2nd_derivatives (CellSimilarity::none, - data_set, - data, - output_data.jacobian_pushed_forward_2nd_derivatives); - maybe_update_jacobian_3rd_derivatives (CellSimilarity::none, - data_set, - data, - output_data.jacobian_3rd_derivatives); - maybe_update_jacobian_pushed_forward_3rd_derivatives (CellSimilarity::none, - data_set, - data, - output_data.jacobian_pushed_forward_3rd_derivatives); - - maybe_compute_face_data (mapping, - cell, face_no, subface_no, quadrature.size(), - quadrature.get_weights(), data, - output_data); + /** + * Do the work of MappingQGeneric::fill_fe_face_values() and + * MappingQGeneric::fill_fe_subface_values() in a generic way, + * using the 'data_set' to differentiate whether we will + * work on a face (and if so, which one) or subface. + */ + template + void + do_fill_fe_face_values (const dealii::MappingQGeneric &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const typename QProjector::DataSetDescriptor data_set, + const Quadrature &quadrature, + const typename dealii::MappingQGeneric::InternalData &data, + internal::FEValues::MappingRelatedData &output_data) + { + maybe_compute_q_points (data_set, + data, + output_data.quadrature_points); + maybe_update_Jacobians (CellSimilarity::none, + data_set, + data); + maybe_update_jacobian_grads (CellSimilarity::none, + data_set, + data, + output_data.jacobian_grads); + maybe_update_jacobian_pushed_forward_grads (CellSimilarity::none, + data_set, + data, + output_data.jacobian_pushed_forward_grads); + maybe_update_jacobian_2nd_derivatives (CellSimilarity::none, + data_set, + data, + output_data.jacobian_2nd_derivatives); + maybe_update_jacobian_pushed_forward_2nd_derivatives (CellSimilarity::none, + data_set, + data, + output_data.jacobian_pushed_forward_2nd_derivatives); + maybe_update_jacobian_3rd_derivatives (CellSimilarity::none, + data_set, + data, + output_data.jacobian_3rd_derivatives); + maybe_update_jacobian_pushed_forward_3rd_derivatives (CellSimilarity::none, + data_set, + data, + output_data.jacobian_pushed_forward_3rd_derivatives); + + maybe_compute_face_data (mapping, + cell, face_no, subface_no, quadrature.size(), + quadrature.get_weights(), data, + output_data); + } } } } @@ -2883,16 +2903,17 @@ fill_fe_face_values (const typename Triangulation::cell_iterator & data.cell_of_current_support_points = cell; } - internal::do_fill_fe_face_values (*this, - cell, face_no, numbers::invalid_unsigned_int, - QProjector::DataSetDescriptor::face (face_no, - cell->face_orientation(face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - quadrature.size()), - quadrature, - data, - output_data); + internal::MappingQGeneric::do_fill_fe_face_values + (*this, + cell, face_no, numbers::invalid_unsigned_int, + QProjector::DataSetDescriptor::face (face_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature.size()), + quadrature, + data, + output_data); } @@ -2928,337 +2949,363 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato data.cell_of_current_support_points = cell; } - internal::do_fill_fe_face_values (*this, - cell, face_no, subface_no, - QProjector::DataSetDescriptor::subface (face_no, subface_no, - cell->face_orientation(face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - quadrature.size(), - cell->subface_case(face_no)), - quadrature, - data, - output_data); + internal::MappingQGeneric::do_fill_fe_face_values + (*this, + cell, face_no, subface_no, + QProjector::DataSetDescriptor::subface (face_no, subface_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature.size(), + cell->subface_case(face_no)), + quadrature, + data, + output_data); } -namespace +namespace internal { - template - void - transform_fields(const ArrayView > &input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView > &output) + namespace MappingQGeneric { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), - ExcInternalError()); - const typename MappingQGeneric::InternalData - &data = static_cast::InternalData &>(mapping_data); - - switch (mapping_type) - { - case mapping_contravariant: + namespace + { + // We cannot query a manifold from the faces of a 1D elements (i.e., + // vertices), which is why we add a specialization for the 3D case here + template + bool check_identical_manifolds_of_quads(const Iterator &) { - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - - for (unsigned int i=0; i::cell_iterator &cell) { - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - Assert (data.update_each & update_volume_elements, - typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); - Assert (rank==1, ExcMessage("Only for rank 1")); - if (rank!=1) - return; - - for (unsigned int i=0; i::faces_per_cell; ++f) + if (&cell->face(f)->get_manifold() != &cell->get_manifold()) + return false; + return true; } - //We still allow this operation as in the - //reference cell Derivatives are Tensor - //rather than DerivativeForm - case mapping_covariant: - { - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - for (unsigned int i=0; i + void + transform_fields(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) + { + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingQGeneric::InternalData + &data = static_cast::InternalData &>(mapping_data); + switch (mapping_type) + { + case mapping_contravariant: + { + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - template - void - transform_gradients(const ArrayView > &input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView > &output) - { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), - ExcInternalError()); - const typename MappingQGeneric::InternalData - &data = static_cast::InternalData &>(mapping_data); + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_covariant_transformation")); - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - Assert (rank==2, ExcMessage("Only for rank 2")); + return; + } - for (unsigned int i=0; i A = - apply_transformation(data.contravariant[i], transpose(input[i]) ); - output[i] = apply_transformation(data.covariant[i], A.transpose() ); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + Assert (rank==1, ExcMessage("Only for rank 1")); + if (rank!=1) + return; + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_covariant_transformation")); + + for (unsigned int i=0; i + void + transform_gradients(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) { - Assert (data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - Assert (rank==2, ExcMessage("Only for rank 2")); + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingQGeneric::InternalData + &data = static_cast::InternalData &>(mapping_data); - for (unsigned int i=0; i A = - apply_transformation(data.covariant[i], transpose(input[i]) ); - output[i] = apply_transformation(data.covariant[i], A.transpose() ); - } + case mapping_contravariant_gradient: + { + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (rank==2, ExcMessage("Only for rank 2")); - return; - } + for (unsigned int i=0; i A = + apply_transformation(data.contravariant[i], transpose(input[i]) ); + output[i] = apply_transformation(data.covariant[i], A.transpose() ); + } - case mapping_piola_gradient: - { - Assert (data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - Assert (data.update_each & update_volume_elements, - typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); - Assert (rank==2, ExcMessage("Only for rank 2")); - - for (unsigned int i=0; i A = - apply_transformation(data.covariant[i], input[i] ); - Tensor<2,spacedim> T = - apply_transformation(data.contravariant[i], A.transpose() ); + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (rank==2, ExcMessage("Only for rank 2")); - output[i] = transpose(T); - output[i] /= data.volume_elements[i]; + for (unsigned int i=0; i A = + apply_transformation(data.covariant[i], transpose(input[i]) ); + output[i] = apply_transformation(data.covariant[i], A.transpose() ); + } + + return; } - return; - } + case mapping_piola_gradient: + { + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + Assert (rank==2, ExcMessage("Only for rank 2")); + + for (unsigned int i=0; i A = + apply_transformation(data.covariant[i], input[i] ); + Tensor<2,spacedim> T = + apply_transformation(data.contravariant[i], A.transpose() ); - default: - Assert(false, ExcNotImplemented()); - } - } + output[i] = transpose(T); + output[i] /= data.volume_elements[i]; + } + + return; + } + default: + Assert(false, ExcNotImplemented()); + } + } - template - void - transform_hessians(const ArrayView > &input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView > &output) - { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), - ExcInternalError()); - const typename MappingQGeneric::InternalData - &data = static_cast::InternalData &>(mapping_data); - switch (mapping_type) - { - case mapping_contravariant_hessian: + template + void + transform_hessians(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) { - Assert (data.update_each & update_covariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingQGeneric::InternalData + &data = static_cast::InternalData &>(mapping_data); - for (unsigned int q=0; q::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + + for (unsigned int q=0; q::ExcAccessToUninitializedField("update_covariant_transformation")); + case mapping_covariant_hessian: + { + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - for (unsigned int q=0; q::ExcAccessToUninitializedField("update_covariant_transformation")); - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - Assert (data.update_each & update_volume_elements, - typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); - - for (unsigned int q=0; q::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + + for (unsigned int q=0; q - void - transform_differential_forms(const ArrayView > &input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView > &output) - { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), - ExcInternalError()); - const typename MappingQGeneric::InternalData - &data = static_cast::InternalData &>(mapping_data); - switch (mapping_type) - { - case mapping_covariant: + template + void + transform_differential_forms(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) { - Assert (data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingQGeneric::InternalData + &data = static_cast::InternalData &>(mapping_data); - for (unsigned int i=0; i::ExcAccessToUninitializedField("update_covariant_transformation")); - return; - } - default: - Assert(false, ExcNotImplemented()); + for (unsigned int i=0; i > &input, const typename Mapping::InternalDataBase &mapping_data, const ArrayView > &output) const { - transform_fields(input, mapping_type, mapping_data, output); + internal::MappingQGeneric::transform_fields(input, mapping_type, mapping_data, output); } @@ -3285,7 +3332,7 @@ transform (const ArrayView > &input, const typename Mapping::InternalDataBase &mapping_data, const ArrayView > &output) const { - transform_differential_forms(input, mapping_type, mapping_data, output); + internal::MappingQGeneric::transform_differential_forms(input, mapping_type, mapping_data, output); } @@ -3301,13 +3348,13 @@ transform (const ArrayView > &input, switch (mapping_type) { case mapping_contravariant: - transform_fields(input, mapping_type, mapping_data, output); + internal::MappingQGeneric::transform_fields(input, mapping_type, mapping_data, output); return; case mapping_piola_gradient: case mapping_contravariant_gradient: case mapping_covariant_gradient: - transform_gradients(input, mapping_type, mapping_data, output); + internal::MappingQGeneric::transform_gradients(input, mapping_type, mapping_data, output); return; default: Assert(false, ExcNotImplemented()); @@ -3378,7 +3425,7 @@ transform (const ArrayView > &input, case mapping_piola_hessian: case mapping_contravariant_hessian: case mapping_covariant_hessian: - transform_hessians(input, mapping_type, mapping_data, output); + internal::MappingQGeneric::transform_hessians(input, mapping_type, mapping_data, output); return; default: Assert(false, ExcNotImplemented()); @@ -3387,28 +3434,6 @@ transform (const ArrayView > &input, -namespace -{ - // We cannot query a manifold from the faces of a 1D elements (i.e., - // vertices), which is why we add a specialization for the 3D case here - template - bool check_identical_manifolds_of_quads(const Iterator &) - { - Assert(false, ExcNotImplemented()); - return true; - } - - bool check_identical_manifolds_of_quads(const Triangulation<3,3>::cell_iterator &cell) - { - for (unsigned int f=0; f::faces_per_cell; ++f) - if (&cell->face(f)->get_manifold() != &cell->get_manifold()) - return false; - return true; - } -} - - - template void MappingQGeneric:: @@ -3642,7 +3667,7 @@ compute_mapping_support_points(const typename Triangulation::cell_ if (&cell->line(l)->get_manifold() != &cell->get_manifold()) all_manifold_ids_are_equal = false; if (dim == 3) - if (check_identical_manifolds_of_quads(cell) == false) + if (internal::MappingQGeneric::check_identical_manifolds_of_quads(cell) == false) all_manifold_ids_are_equal = false; } -- 2.39.5