From 0faad4e08019d5f9909e70ebee390198ca978d32 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth <bangerth@math.tamu.edu> Date: Thu, 3 Sep 2015 08:15:26 -0500 Subject: [PATCH] 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. --- include/deal.II/fe/mapping_q.h | 20 +--- include/deal.II/fe/mapping_q1.h | 22 +++-- source/fe/mapping_q.cc | 131 ++---------------------- source/fe/mapping_q1.cc | 170 +++++++++++++++++++++++++------- 4 files changed, 155 insertions(+), 188 deletions(-) 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<dim,spacedim> &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<Point<dim> > &unit_points, - typename MappingQ1<dim,spacedim>::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<unsigned int> 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<dim> &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<Point<dim> > &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 <tt>data.polynomial_degree</tt> - * is one or is greater than that. - */ - void compute_shapes (const std::vector<Point<dim> > &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 <deal.II/grid/tria_iterator.h> #include <deal.II/grid/tria_boundary.h> #include <deal.II/dofs/dof_accessor.h> -#include <deal.II/fe/fe_tools.h> #include <deal.II/fe/mapping_q.h> #include <deal.II/fe/fe_q.h> #include <deal.II/fe/fe_values.h> @@ -55,22 +54,6 @@ MappingQ<dim,spacedim>::InternalData::memory_consumption () const -namespace -{ - template <int dim> - std::vector<unsigned int> - get_dpo_vector (const unsigned int degree) - { - std::vector<unsigned int> dpo(dim+1, 1U); - for (unsigned int i=1; i<dpo.size(); ++i) - dpo[i]=dpo[i-1]*(degree-1); - return dpo; - } -} - - - - template<int dim, int spacedim> MappingQ<dim,spacedim>::MappingQ (const unsigned int p, const bool use_mapping_q_on_all_cells) @@ -81,23 +64,12 @@ MappingQ<dim,spacedim>::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<dim>(degree+1)), - renumber(FETools:: - lexicographic_to_hierarchic_numbering ( - FiniteElementData<dim> (get_dpo_vector<dim>(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<dim> - (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<dim,spacedim>::MappingQ (const MappingQ<dim,spacedim> &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<dim> (*mapping.tensor_pols); laplace_on_quad_vector=mapping.laplace_on_quad_vector; laplace_on_hex_vector=mapping.laplace_on_hex_vector; } -template<int dim, int spacedim> -MappingQ<dim,spacedim>::~MappingQ () -{ - delete tensor_pols; -} - - -template<int dim, int spacedim> -void -MappingQ<dim,spacedim>::compute_shapes (const std::vector<Point<dim> > &unit_points, - typename MappingQ1<dim,spacedim>::InternalData &data) const -{ - const unsigned int n_points=unit_points.size(); - std::vector<double> values; - std::vector<Tensor<1,dim> > 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<Tensor<2,dim> > 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<Tensor<3,dim> > 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<Tensor<4,dim> > 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<n_points; ++point) - { - tensor_pols->compute(unit_points[point], values, grads, grad2, grad3, grad4); - - if (data.shape_values.size()!=0) - for (unsigned int i=0; i<n_shape_functions; ++i) - data.shape(point,renumber[i]) = values[i]; - - if (data.shape_derivatives.size()!=0) - for (unsigned int i=0; i<n_shape_functions; ++i) - data.derivative(point,renumber[i]) = grads[i]; - - if (data.shape_second_derivatives.size()!=0) - for (unsigned int i=0; i<n_shape_functions; ++i) - data.second_derivative(point,renumber[i]) = grad2[i]; - - if (data.shape_third_derivatives.size()!=0) - for (unsigned int i=0; i<n_shape_functions; ++i) - data.third_derivative(point,renumber[i]) = grad3[i]; - - if (data.shape_fourth_derivatives.size()!=0) - for (unsigned int i=0; i<n_shape_functions; ++i) - data.fourth_derivative(point,renumber[i]) = grad4[i]; - } -} - - template<int dim, int spacedim> typename MappingQ<dim,spacedim>::InternalData * @@ -243,9 +126,9 @@ MappingQ<dim,spacedim>::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<dim,spacedim>::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<dim,spacedim>::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<dim,spacedim>::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<dim,spacedim>::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<dim,spacedim>::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<dim,spacedim>::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<long double> 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 <deal.II/grid/tria.h> #include <deal.II/grid/tria_iterator.h> #include <deal.II/dofs/dof_accessor.h> +#include <deal.II/fe/fe_tools.h> +#include <deal.II/fe/fe.h> #include <deal.II/fe/fe_values.h> #include <deal.II/fe/mapping_q1.h> #include <deal.II/fe/mapping_q.h> @@ -326,9 +328,9 @@ namespace internal template <int spacedim> void - compute_shapes (const unsigned int n_shape_functions, - const std::vector<Point<1> > &unit_points, - typename dealii::MappingQ1<1,spacedim>::InternalData &data) + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector<Point<1> > &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 <int spacedim> void - compute_shapes (const unsigned int n_shape_functions, - const std::vector<Point<2> > &unit_points, - typename dealii::MappingQ1<2,spacedim>::InternalData &data) + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector<Point<2> > &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 <int spacedim> void - compute_shapes (const unsigned int n_shape_functions, - const std::vector<Point<3> > &unit_points, - typename dealii::MappingQ1<3,spacedim>::InternalData &data) + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector<Point<3> > &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 <int dim> + std::vector<unsigned int> + get_dpo_vector (const unsigned int degree) + { + std::vector<unsigned int> dpo(dim+1, 1U); + for (unsigned int i=1; i<dpo.size(); ++i) + dpo[i]=dpo[i-1]*(degree-1); + return dpo; + } +} + + + template<int dim, int spacedim> void -MappingQ1<dim,spacedim>::compute_shapes (const std::vector<Point<dim> > &unit_points, - InternalData &data) const +MappingQ1<dim,spacedim>::InternalData:: +compute_shape_function_values (const std::vector<Point<dim> > &unit_points) { - internal::MappingQ1:: - compute_shapes<spacedim> (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<spacedim> (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<dim> + 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<unsigned int> + renumber (FETools:: + lexicographic_to_hierarchic_numbering ( + FiniteElementData<dim> (get_dpo_vector<dim>(polynomial_degree), 1, + polynomial_degree))); + + std::vector<double> values; + std::vector<Tensor<1,dim> > 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<Tensor<2,dim> > 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<Tensor<3,dim> > 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<Tensor<4,dim> > 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_points; ++point) + { + tensor_pols.compute(unit_points[point], values, grads, grad2, grad3, grad4); + + if (shape_values.size()!=0) + for (unsigned int i=0; i<n_shape_functions; ++i) + shape(point,renumber[i]) = values[i]; + + if (shape_derivatives.size()!=0) + for (unsigned int i=0; i<n_shape_functions; ++i) + derivative(point,renumber[i]) = grads[i]; + + if (shape_second_derivatives.size()!=0) + for (unsigned int i=0; i<n_shape_functions; ++i) + second_derivative(point,renumber[i]) = grad2[i]; + + if (shape_third_derivatives.size()!=0) + for (unsigned int i=0; i<n_shape_functions; ++i) + third_derivative(point,renumber[i]) = grad3[i]; + + if (shape_fourth_derivatives.size()!=0) + for (unsigned int i=0; i<n_shape_functions; ++i) + fourth_derivative(point,renumber[i]) = grad4[i]; + } + } } @@ -732,7 +843,7 @@ MappingQ1<dim,spacedim>::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<dim,spacedim>::get_face_data (const UpdateFlags update_flags, data->initialize_face (requires_update_flags(update_flags), QProjector<dim>::project_to_all_faces(quadrature), quadrature.size()); - compute_shapes (QProjector<dim>::project_to_all_faces(quadrature).get_points(), - *data); + data->compute_shape_function_values (QProjector<dim>::project_to_all_faces(quadrature).get_points()); return data; } @@ -765,8 +875,7 @@ MappingQ1<dim,spacedim>::get_subface_data (const UpdateFlags update_flags, data->initialize_face (requires_update_flags(update_flags), QProjector<dim>::project_to_all_subfaces(quadrature), quadrature.size()); - compute_shapes (QProjector<dim>::project_to_all_subfaces(quadrature).get_points(), - *data); + data->compute_shape_function_values (QProjector<dim>::project_to_all_subfaces(quadrature).get_points()); return data; @@ -2513,10 +2622,7 @@ transform_real_to_unit_cell_internal Point<dim> p_unit = initial_p_unit; - if (dynamic_cast<typename MappingQ<dim,spacedim>::InternalData *>(&mdata) == 0) - this->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata); - else - dynamic_cast<const MappingQ<dim,spacedim>*>(this)->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata); + mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit)); Point<spacedim> 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<typename MappingQ<dim,spacedim>::InternalData *>(&mdata) == 0) - this->compute_shapes(std::vector<Point<dim> > (1, p_unit_trial), mdata); - else - dynamic_cast<const MappingQ<dim,spacedim>*>(this)->compute_shapes(std::vector<Point<dim> > (1, p_unit_trial), mdata); - + mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit_trial)); // f(x) Point<spacedim> 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<typename MappingQ<dim,spacedim>::InternalData *>(&mdata) == 0) - this->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata); - else - dynamic_cast<const MappingQ<dim,spacedim>*>(this)->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata); + mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit)); + for (unsigned int k=0; k<mdata.n_shape_functions; ++k) { const Tensor<1,dim1> &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<typename MappingQ<dim,spacedim>::InternalData *>(&mdata) == 0) - this->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata); - else - dynamic_cast<const MappingQ<dim,spacedim>*>(this)->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata); + mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit)); + for (unsigned int k=0; k<mdata.n_shape_functions; ++k) { const Tensor<1,dim1> &grad_phi_k = mdata.derivative(0,k); -- 2.39.5