From: Wolfgang Bangerth Date: Thu, 3 Sep 2015 13:15:26 +0000 (-0500) Subject: Move the code that initializes the Mapping{Q,Q1}::InternalData into that structure. X-Git-Tag: v8.4.0-rc2~474^2~9 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0faad4e08019d5f9909e70ebee390198ca978d32;p=dealii.git Move the code that initializes the Mapping{Q,Q1}::InternalData into that structure. This allows us to get rid of the two implementations of compute_shapes(). It also significantly simplifies the logic which of the two functions actually need to be called, since the InternalData object actually knows that itself. --- diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index f6b588a1ec..2d4914eaaf 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -72,17 +72,10 @@ public: const bool use_mapping_q_on_all_cells = false); /** - * Copy constructor. Performs a deep copy, i.e. duplicates what #tensor_pols - * points to instead of simply copying the #tensor_pols pointer as done by a - * default copy constructor. + * Copy constructor. */ MappingQ (const MappingQ &mapping); - /** - * Destructor. - */ - virtual ~MappingQ (); - /** * Transforms the point @p p on the unit cell to the point @p p_real on the * real cell @p cell and returns @p p_real. @@ -333,10 +326,7 @@ protected: /** * @} */ -public: - void - compute_shapes (const std::vector > &unit_points, - typename MappingQ1::InternalData &data) const; + protected: /** @@ -462,12 +452,6 @@ protected: */ const unsigned int n_shape_functions; - /** - * Mapping from lexicographic to to the Qp shape function numbering. Its - * size is @p dofs_per_cell. - */ - const std::vector renumber; - /** * If this flag is set @p true then @p MappingQ is used on all cells, not * only on boundary cells. diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h index 6a32f40ceb..06d4beff9b 100644 --- a/include/deal.II/fe/mapping_q1.h +++ b/include/deal.II/fe/mapping_q1.h @@ -184,6 +184,18 @@ public: const Quadrature &quadrature, const unsigned int n_original_q_points); + /** + * Compute the values and/or derivatives of the shape functions + * used for the mapping. + * + * Which values, derivatives, or higher order derivatives are + * computed is determined by which of the member arrays have + * nonzero sizes. They is typically set to their appropriate sizes + * by the initialize() and initialize_face() functions. + */ + void compute_shape_function_values (const std::vector > &unit_points); + + /** * Shape function at quadrature point. Shape functions are in tensor * product order, so vertices must be reordered to obtain transformation. @@ -465,16 +477,6 @@ protected: DataSetDescriptor; - /** - * Compute shape values and/or derivatives. - * - * Calls either the @p compute_shapes_virtual of this class or that of the - * derived class, depending on whether data.polynomial_degree - * is one or is greater than that. - */ - void compute_shapes (const std::vector > &unit_points, - InternalData &data) const; - /** * Transforms a point @p p on the unit cell to the point @p p_real on the * real cell @p cell and returns @p p_real. diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index e2d6ea4737..158437471b 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -24,7 +24,6 @@ #include #include #include -#include #include #include #include @@ -55,22 +54,6 @@ MappingQ::InternalData::memory_consumption () const -namespace -{ - template - std::vector - get_dpo_vector (const unsigned int degree) - { - std::vector dpo(dim+1, 1U); - for (unsigned int i=1; i MappingQ::MappingQ (const unsigned int p, const bool use_mapping_q_on_all_cells) @@ -81,23 +64,12 @@ MappingQ::MappingQ (const unsigned int p, ((dim==2) ? 4+4*(degree-1) : 8+12*(degree-1)+6*(degree-1)*(degree-1))), - tensor_pols(0), n_shape_functions(Utilities::fixed_power(degree+1)), - renumber(FETools:: - lexicographic_to_hierarchic_numbering ( - FiniteElementData (get_dpo_vector(degree), 1, - degree))), use_mapping_q_on_all_cells (use_mapping_q_on_all_cells || (dim != spacedim)), feq(degree), line_support_points(degree+1) { - // Construct the tensor product polynomials used as shape functions for the - // Qp mapping of cells at the boundary. - tensor_pols = new TensorProductPolynomials - (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points())); - Assert (n_shape_functions==tensor_pols->n(), - ExcInternalError()); Assert(n_inner+n_outer==n_shape_functions, ExcInternalError()); // build laplace_on_quad_vector @@ -118,105 +90,16 @@ MappingQ::MappingQ (const MappingQ &mapping) degree(mapping.degree), n_inner(mapping.n_inner), n_outer(mapping.n_outer), - tensor_pols(0), n_shape_functions(mapping.n_shape_functions), - renumber(mapping.renumber), use_mapping_q_on_all_cells (mapping.use_mapping_q_on_all_cells), feq(degree), line_support_points(degree+1) { - tensor_pols=new TensorProductPolynomials (*mapping.tensor_pols); laplace_on_quad_vector=mapping.laplace_on_quad_vector; laplace_on_hex_vector=mapping.laplace_on_hex_vector; } -template -MappingQ::~MappingQ () -{ - delete tensor_pols; -} - - -template -void -MappingQ::compute_shapes (const std::vector > &unit_points, - typename MappingQ1::InternalData &data) const -{ - const unsigned int n_points=unit_points.size(); - 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); - } - -// // dummy variable of size 0 - 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; pointcompute(unit_points[point], values, grads, grad2, grad3, grad4); - - if (data.shape_values.size()!=0) - for (unsigned int i=0; i typename MappingQ::InternalData * @@ -243,9 +126,9 @@ MappingQ::get_data (const UpdateFlags update_flags, tasks.join_all (); // TODO: parallelize this as well - this->compute_shapes (quadrature.get_points(), *data); + data->compute_shape_function_values (quadrature.get_points()); if (!use_mapping_q_on_all_cells) - this->MappingQ1::compute_shapes (quadrature.get_points(), data->mapping_q1_data); + data->mapping_q1_data.compute_shape_function_values (quadrature.get_points()); return data; @@ -279,9 +162,9 @@ MappingQ::get_face_data (const UpdateFlags update_flags, tasks.join_all (); // TODO: parallelize this as well - this->compute_shapes (q.get_points(), *data); + data->compute_shape_function_values (q.get_points()); if (!use_mapping_q_on_all_cells) - this->MappingQ1::compute_shapes (q.get_points(), data->mapping_q1_data); + data->mapping_q1_data.compute_shape_function_values (q.get_points()); return data; } @@ -314,9 +197,9 @@ MappingQ::get_subface_data (const UpdateFlags update_flags, tasks.join_all (); // TODO: parallelize this as well - this->compute_shapes (q.get_points(), *data); + data->compute_shape_function_values (q.get_points()); if (!use_mapping_q_on_all_cells) - this->MappingQ1::compute_shapes (q.get_points(), data->mapping_q1_data); + data->mapping_q1_data.compute_shape_function_values (q.get_points()); return data; @@ -616,7 +499,7 @@ MappingQ::compute_laplace_vector(Table<2,double> &lvs) const InternalData quadrature_data(degree); quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points); - this->compute_shapes(quadrature.get_points(), quadrature_data); + quadrature_data.compute_shape_function_values(quadrature.get_points()); // Compute the stiffness matrix of the inner dofs FullMatrix S(n_inner); diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index 69ce4b5d8f..ca19ad1b8e 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -24,6 +24,8 @@ #include #include #include +#include +#include #include #include #include @@ -326,9 +328,9 @@ namespace internal template void - compute_shapes (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingQ1<1,spacedim>::InternalData &data) + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector > &unit_points, + typename dealii::MappingQ1<1,spacedim>::InternalData &data) { (void)n_shape_functions; const unsigned int n_points=unit_points.size(); @@ -391,9 +393,9 @@ namespace internal template void - compute_shapes (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingQ1<2,spacedim>::InternalData &data) + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector > &unit_points, + typename dealii::MappingQ1<2,spacedim>::InternalData &data) { (void)n_shape_functions; const unsigned int n_points=unit_points.size(); @@ -469,9 +471,9 @@ namespace internal template void - compute_shapes (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingQ1<3,spacedim>::InternalData &data) + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector > &unit_points, + typename dealii::MappingQ1<3,spacedim>::InternalData &data) { (void)n_shape_functions; const unsigned int n_points=unit_points.size(); @@ -653,14 +655,123 @@ namespace internal } +namespace +{ + template + std::vector + get_dpo_vector (const unsigned int degree) + { + std::vector dpo(dim+1, 1U); + for (unsigned int i=1; i void -MappingQ1::compute_shapes (const std::vector > &unit_points, - InternalData &data) const +MappingQ1::InternalData:: +compute_shape_function_values (const std::vector > &unit_points) { - internal::MappingQ1:: - compute_shapes (n_shape_functions, - unit_points, data); + // 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, + 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 QGaussLobatto<1> line_support_points (polynomial_degree + 1); + 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::get_data (const UpdateFlags update_flags, { InternalData *data = new InternalData(1); data->initialize (requires_update_flags(update_flags), q, q.size()); - compute_shapes (q.get_points(), *data); + data->compute_shape_function_values (q.get_points()); return data; } @@ -748,8 +859,7 @@ MappingQ1::get_face_data (const UpdateFlags update_flags, data->initialize_face (requires_update_flags(update_flags), QProjector::project_to_all_faces(quadrature), quadrature.size()); - compute_shapes (QProjector::project_to_all_faces(quadrature).get_points(), - *data); + data->compute_shape_function_values (QProjector::project_to_all_faces(quadrature).get_points()); return data; } @@ -765,8 +875,7 @@ MappingQ1::get_subface_data (const UpdateFlags update_flags, data->initialize_face (requires_update_flags(update_flags), QProjector::project_to_all_subfaces(quadrature), quadrature.size()); - compute_shapes (QProjector::project_to_all_subfaces(quadrature).get_points(), - *data); + data->compute_shape_function_values (QProjector::project_to_all_subfaces(quadrature).get_points()); return data; @@ -2513,10 +2622,7 @@ transform_real_to_unit_cell_internal Point p_unit = initial_p_unit; - if (dynamic_cast::InternalData *>(&mdata) == 0) - this->compute_shapes(std::vector > (1, p_unit), mdata); - else - dynamic_cast*>(this)->compute_shapes(std::vector > (1, p_unit), mdata); + mdata.compute_shape_function_values(std::vector > (1, p_unit)); Point p_real = transform_unit_to_real_cell_internal(mdata); Tensor<1,spacedim> f = p_real-p; @@ -2604,11 +2710,7 @@ transform_real_to_unit_cell_internal // shape values and derivatives // at new p_unit point - if (dynamic_cast::InternalData *>(&mdata) == 0) - this->compute_shapes(std::vector > (1, p_unit_trial), mdata); - else - dynamic_cast*>(this)->compute_shapes(std::vector > (1, p_unit_trial), mdata); - + mdata.compute_shape_function_values(std::vector > (1, p_unit_trial)); // f(x) Point p_real_trial = transform_unit_to_real_cell_internal(mdata); @@ -2750,10 +2852,8 @@ transform_real_to_unit_cell_internal_codim1 Tensor<2,dim1> df; //Evaluate first and second derivatives - if (dynamic_cast::InternalData *>(&mdata) == 0) - this->compute_shapes(std::vector > (1, p_unit), mdata); - else - dynamic_cast*>(this)->compute_shapes(std::vector > (1, p_unit), mdata); + mdata.compute_shape_function_values(std::vector > (1, p_unit)); + for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); @@ -2805,10 +2905,8 @@ transform_real_to_unit_cell_internal_codim1 D2F[j][l].clear(); } - if (dynamic_cast::InternalData *>(&mdata) == 0) - this->compute_shapes(std::vector > (1, p_unit), mdata); - else - dynamic_cast*>(this)->compute_shapes(std::vector > (1, p_unit), mdata); + mdata.compute_shape_function_values(std::vector > (1, p_unit)); + for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k);