From 64d0b68c934c055ed945244f52101b75ab7a8c59 Mon Sep 17 00:00:00 2001 From: Karl Ljungkvist Date: Wed, 23 Nov 2016 14:43:07 +0100 Subject: [PATCH] change a runtime branch into a compile time one by templatization this avoids the special hardcoded case for when dim==spacedim to be instantiated also for other combinations of dim and spacedim. --- source/fe/mapping_q_generic.cc | 308 ++++++++++++++++++--------------- 1 file changed, 172 insertions(+), 136 deletions(-) diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 3cae9d3564..ae48043d15 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -40,6 +40,20 @@ DEAL_II_NAMESPACE_OPEN + +namespace +{ + template + std::vector + get_dpo_vector (const unsigned int degree) + { + std::vector dpo(dim+1, 1U); + for (unsigned int i=1; i + + + template + void compute_shape_function_values_general (const unsigned int n_shape_functions, + const std::vector > &unit_points, + typename dealii::MappingQGeneric::InternalData &data) + { + const unsigned int n_points=unit_points.size(); + + // Construct the tensor product polynomials used as shape functions for the + // Qp mapping of cells at the boundary. + const TensorProductPolynomials + tensor_pols (Polynomials::generate_complete_Lagrange_basis(data.line_support_points.get_points())); + Assert (n_shape_functions==tensor_pols.n(), + ExcInternalError()); + + // 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(data.polynomial_degree), 1, + data.polynomial_degree))); + + std::vector values; + std::vector > grads; + if (data.shape_values.size()!=0) + { + Assert(data.shape_values.size()==n_shape_functions*n_points, + ExcInternalError()); + values.resize(n_shape_functions); + } + if (data.shape_derivatives.size()!=0) + { + Assert(data.shape_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + grads.resize(n_shape_functions); + } + + std::vector > grad2; + if (data.shape_second_derivatives.size()!=0) + { + Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + grad2.resize(n_shape_functions); + } + + std::vector > grad3; + if (data.shape_third_derivatives.size()!=0) + { + Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + grad3.resize(n_shape_functions); + } + + std::vector > grad4; + if (data.shape_fourth_derivatives.size()!=0) + { + Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + grad4.resize(n_shape_functions); + } + + + if (data.shape_values.size()!=0 || + data.shape_derivatives.size()!=0 || + data.shape_second_derivatives.size()!=0 || + data.shape_third_derivatives.size()!=0 || + data.shape_fourth_derivatives.size()!=0 ) + for (unsigned int point=0; point > &unit_points, - typename dealii::MappingQGeneric<1,spacedim>::InternalData &data) + compute_shape_function_values_hardcode (const unsigned int n_shape_functions, + const std::vector > &unit_points, + typename dealii::MappingQGeneric<1,1>::InternalData &data) { (void)n_shape_functions; const unsigned int n_points=unit_points.size(); @@ -356,10 +464,6 @@ namespace internal } if (data.shape_second_derivatives.size()!=0) { - // the following may or may not - // work if dim != spacedim - Assert (spacedim == 1, ExcNotImplemented()); - Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, ExcInternalError()); data.second_derivative(k,0)[0][0] = 0; @@ -367,9 +471,6 @@ namespace internal } if (data.shape_third_derivatives.size()!=0) { - // if lower order derivative don't work, neither should this - Assert (spacedim == 1, ExcNotImplemented()); - Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points, ExcInternalError()); @@ -379,9 +480,6 @@ namespace internal } if (data.shape_fourth_derivatives.size()!=0) { - // if lower order derivative don't work, neither should this - Assert (spacedim == 1, ExcNotImplemented()); - Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points, ExcInternalError()); @@ -393,12 +491,12 @@ namespace internal } - template void - compute_shape_function_values (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingQGeneric<2,spacedim>::InternalData &data) + compute_shape_function_values_hardcode (const unsigned int n_shape_functions, + const std::vector > &unit_points, + typename dealii::MappingQGeneric<2,2>::InternalData &data) { + (void)n_shape_functions; const unsigned int n_points=unit_points.size(); for (unsigned int k = 0 ; k < n_points ; ++k) @@ -471,11 +569,10 @@ namespace internal - template void - compute_shape_function_values (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingQGeneric<3,spacedim>::InternalData &data) + compute_shape_function_values_hardcode (const unsigned int n_shape_functions, + const std::vector > &unit_points, + typename dealii::MappingQGeneric<3,3>::InternalData &data) { (void)n_shape_functions; const unsigned int n_points=unit_points.size(); @@ -529,10 +626,6 @@ namespace internal } if (data.shape_second_derivatives.size()!=0) { - // the following may or may not - // work if dim != spacedim - Assert (spacedim == 3, ExcNotImplemented()); - Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, ExcInternalError()); data.second_derivative(k,0)[0][0] = 0; @@ -613,9 +706,6 @@ namespace internal } if (data.shape_third_derivatives.size()!=0) { - // if lower order derivative don't work, neither should this - Assert (spacedim == 3, ExcNotImplemented()); - Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points, ExcInternalError()); @@ -642,9 +732,6 @@ namespace internal } if (data.shape_fourth_derivatives.size()!=0) { - // if lower order derivative don't work, neither should this - Assert (spacedim == 3, ExcNotImplemented()); - Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points, ExcInternalError()); Tensor<4,3> zero; @@ -829,124 +916,73 @@ initialize_face (const UpdateFlags update_flags, -namespace -{ - template - std::vector - get_dpo_vector (const unsigned int degree) - { - std::vector dpo(dim+1, 1U); - for (unsigned int i=1; i +template<> void -MappingQGeneric::InternalData:: -compute_shape_function_values (const std::vector > &unit_points) +MappingQGeneric<1,1>::InternalData:: +compute_shape_function_values (const std::vector > &unit_points) { // if the polynomial degree is one, then we can simplify code a bit // by using hard-coded shape functions. - if ((polynomial_degree == 1) - && - (dim == spacedim)) - internal::MappingQ1::compute_shape_function_values (n_shape_functions, + if (polynomial_degree == 1) + internal::MappingQ1::compute_shape_function_values_hardcode (n_shape_functions, unit_points, *this); else - // otherwise ask an object that describes the polynomial space { - const unsigned int n_points=unit_points.size(); - - // Construct the tensor product polynomials used as shape functions for the - // Qp mapping of cells at the boundary. - const TensorProductPolynomials - tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points())); - Assert (n_shape_functions==tensor_pols.n(), - ExcInternalError()); - - // 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))); - - std::vector values; - std::vector > grads; - if (shape_values.size()!=0) - { - Assert(shape_values.size()==n_shape_functions*n_points, - ExcInternalError()); - values.resize(n_shape_functions); - } - if (shape_derivatives.size()!=0) - { - Assert(shape_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - grads.resize(n_shape_functions); - } - - std::vector > grad2; - if (shape_second_derivatives.size()!=0) - { - Assert(shape_second_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - grad2.resize(n_shape_functions); - } - - std::vector > grad3; - if (shape_third_derivatives.size()!=0) - { - Assert(shape_third_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - grad3.resize(n_shape_functions); - } - - std::vector > grad4; - if (shape_fourth_derivatives.size()!=0) - { - Assert(shape_fourth_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - grad4.resize(n_shape_functions); - } - - - if (shape_values.size()!=0 || - shape_derivatives.size()!=0 || - shape_second_derivatives.size()!=0 || - shape_third_derivatives.size()!=0 || - shape_fourth_derivatives.size()!=0 ) - for (unsigned int point=0; point(n_shape_functions, + unit_points,*this); + } +} - if (shape_third_derivatives.size()!=0) - for (unsigned int i=0; i +void +MappingQGeneric<2,2>::InternalData:: +compute_shape_function_values (const std::vector > &unit_points) +{ + // if the polynomial degree is one, then we can simplify code a bit + // by using hard-coded shape functions. + if (polynomial_degree == 1) + internal::MappingQ1::compute_shape_function_values_hardcode (n_shape_functions, + unit_points, *this); + else + { + // otherwise ask an object that describes the polynomial space + internal::MappingQ1::compute_shape_function_values_general<2,2>(n_shape_functions, + unit_points,*this); + } +} - if (shape_fourth_derivatives.size()!=0) - for (unsigned int i=0; i +void +MappingQGeneric<3,3>::InternalData:: +compute_shape_function_values (const std::vector > &unit_points) +{ + // if the polynomial degree is one, then we can simplify code a bit + // by using hard-coded shape functions. + if (polynomial_degree == 1) + internal::MappingQ1::compute_shape_function_values_hardcode (n_shape_functions, + unit_points, *this); + else + { + // otherwise ask an object that describes the polynomial space + internal::MappingQ1::compute_shape_function_values_general<3,3>(n_shape_functions, + unit_points,*this); } } +template +void +MappingQGeneric::InternalData:: +compute_shape_function_values (const std::vector > &unit_points) +{ + // for non-matching combinations of dim and spacedim, just run the general + // case + internal::MappingQ1::compute_shape_function_values_general(n_shape_functions, + unit_points,*this); +} + namespace { -- 2.39.5