From 5c6e0f35ff9f2afb5fe0cb6ce76e70ef687cac76 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 2 Sep 2015 17:44:59 -0500 Subject: [PATCH] Get rid of the awkward dispatch mechanism in MappingQ*::compute_shapes(). In particular, remove the manual dispatch to either the Q1 or Qp functions. --- include/deal.II/fe/mapping_q.h | 15 +++---- include/deal.II/fe/mapping_q1.h | 17 +++---- source/fe/mapping_q.cc | 28 +++++------- source/fe/mapping_q1.cc | 80 +++++++++++++++------------------ 4 files changed, 62 insertions(+), 78 deletions(-) diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 317eb7953a..f6b588a1ec 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -261,8 +261,7 @@ public: /** * Constructor. */ - InternalData (const unsigned int polynomial_degree, - const unsigned int n_shape_functions); + InternalData (const unsigned int polynomial_degree); /** @@ -334,13 +333,11 @@ protected: /** * @} */ - - /** - * Compute shape values and/or derivatives. - */ - virtual void - compute_shapes_virtual (const std::vector > &unit_points, - typename MappingQ1::InternalData &data) const; +public: + void + compute_shapes (const std::vector > &unit_points, + typename MappingQ1::InternalData &data) const; +protected: /** * This function is needed by the constructor of diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h index 082f5e7466..6a32f40ceb 100644 --- a/include/deal.II/fe/mapping_q1.h +++ b/include/deal.II/fe/mapping_q1.h @@ -161,10 +161,10 @@ public: { public: /** - * Constructor. Pass the number of shape functions. + * Constructor. The argument denotes the polynomial degree of + * the mapping to which this object will correspond. */ - InternalData(const unsigned int polynomial_degree, - const unsigned int n_shape_functions); + InternalData(const unsigned int polynomial_degree); /** * Initialize the object's member variables related to cell data @@ -315,8 +315,11 @@ public: * the number of vertices per cell. However, since also derived classes * use this class (e.g. the Mapping_Q() class), the number of shape * functions may also be different. + * + * In general, it is $(p+1)^\text{dim}$, where $p$ is the + * polynomial degree of the mapping. */ - unsigned int n_shape_functions; + const unsigned int n_shape_functions; /** * Tensors of covariant transformation at each of the quadrature points. @@ -472,12 +475,6 @@ protected: void compute_shapes (const std::vector > &unit_points, InternalData &data) const; - /** - * Compute shape values and/or derivatives. - */ - virtual void compute_shapes_virtual (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 08068cc01c..e2d6ea4737 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -35,14 +35,11 @@ DEAL_II_NAMESPACE_OPEN template -MappingQ::InternalData::InternalData (const unsigned int polynomial_degree, - const unsigned int n_shape_functions) +MappingQ::InternalData::InternalData (const unsigned int polynomial_degree) : - MappingQ1::InternalData(polynomial_degree, - n_shape_functions), + MappingQ1::InternalData(polynomial_degree), use_mapping_q1_on_current_cell(false), - mapping_q1_data(1, - GeometryInfo::vertices_per_cell) + mapping_q1_data(1) {} @@ -141,11 +138,10 @@ MappingQ::~MappingQ () } - template void -MappingQ::compute_shapes_virtual (const std::vector > &unit_points, - typename MappingQ1::InternalData &data) const +MappingQ::compute_shapes (const std::vector > &unit_points, + typename MappingQ1::InternalData &data) const { const unsigned int n_points=unit_points.size(); std::vector values; @@ -227,7 +223,7 @@ typename MappingQ::InternalData * MappingQ::get_data (const UpdateFlags update_flags, const Quadrature &quadrature) const { - InternalData *data = new InternalData(degree, n_shape_functions); + InternalData *data = new InternalData(degree); // fill the data of both the Q_p and the Q_1 objects in parallel Threads::TaskGroup<> tasks; @@ -249,7 +245,7 @@ MappingQ::get_data (const UpdateFlags update_flags, // TODO: parallelize this as well this->compute_shapes (quadrature.get_points(), *data); if (!use_mapping_q_on_all_cells) - this->compute_shapes (quadrature.get_points(), data->mapping_q1_data); + this->MappingQ1::compute_shapes (quadrature.get_points(), data->mapping_q1_data); return data; @@ -262,7 +258,7 @@ typename Mapping::InternalDataBase * MappingQ::get_face_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { - InternalData *data = new InternalData(degree, n_shape_functions); + InternalData *data = new InternalData(degree); const Quadrature q (QProjector::project_to_all_faces(quadrature)); // fill the data of both the Q_p and the Q_1 objects in parallel @@ -285,7 +281,7 @@ MappingQ::get_face_data (const UpdateFlags update_flags, // TODO: parallelize this as well this->compute_shapes (q.get_points(), *data); if (!use_mapping_q_on_all_cells) - this->compute_shapes (q.get_points(), data->mapping_q1_data); + this->MappingQ1::compute_shapes (q.get_points(), data->mapping_q1_data); return data; } @@ -297,7 +293,7 @@ typename Mapping::InternalDataBase * MappingQ::get_subface_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { - InternalData *data = new InternalData(degree, n_shape_functions); + InternalData *data = new InternalData(degree); const Quadrature q (QProjector::project_to_all_subfaces(quadrature)); // fill the data of both the Q_p and the Q_1 objects in parallel @@ -320,7 +316,7 @@ MappingQ::get_subface_data (const UpdateFlags update_flags, // TODO: parallelize this as well this->compute_shapes (q.get_points(), *data); if (!use_mapping_q_on_all_cells) - this->compute_shapes (q.get_points(), data->mapping_q1_data); + this->MappingQ1::compute_shapes (q.get_points(), data->mapping_q1_data); return data; @@ -618,7 +614,7 @@ MappingQ::compute_laplace_vector(Table<2,double> &lvs) const const QGauss quadrature(degree+1); const unsigned int n_q_points=quadrature.size(); - InternalData quadrature_data(degree, n_shape_functions); + InternalData quadrature_data(degree); quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points); this->compute_shapes(quadrature.get_points(), quadrature_data); diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index e6a176f48a..69ce4b5d8f 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -42,16 +42,11 @@ const unsigned int MappingQ1::n_shape_functions; template -MappingQ1::InternalData::InternalData (const unsigned int polynomial_degree, - const unsigned int n_shape_functions) +MappingQ1::InternalData::InternalData (const unsigned int polynomial_degree) : polynomial_degree (polynomial_degree), - n_shape_functions (n_shape_functions) -{ - Assert (n_shape_functions == - Utilities::fixed_power(polynomial_degree+1), - ExcInternalError()); -} + n_shape_functions (Utilities::fixed_power(polynomial_degree+1)) +{} @@ -205,20 +200,6 @@ MappingQ1::MappingQ1 () -template -void -MappingQ1::compute_shapes (const std::vector > &unit_points, - InternalData &data) const -{ - // choose either the function implemented in this class, or whatever - // a virtual function call resolves to - if (dynamic_cast::InternalData *>(&data) == 0) - MappingQ1::compute_shapes_virtual(unit_points, data); - else - compute_shapes_virtual(unit_points, data); -} - - namespace internal { namespace MappingQ1 @@ -345,9 +326,9 @@ namespace internal template void - compute_shapes_virtual (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingQ1<1,spacedim>::InternalData &data) + compute_shapes (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(); @@ -410,9 +391,9 @@ namespace internal template void - compute_shapes_virtual (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingQ1<2,spacedim>::InternalData &data) + compute_shapes (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(); @@ -488,9 +469,9 @@ namespace internal template void - compute_shapes_virtual (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingQ1<3,spacedim>::InternalData &data) + compute_shapes (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(); @@ -674,13 +655,12 @@ namespace internal template void -MappingQ1:: -compute_shapes_virtual (const std::vector > &unit_points, - InternalData &data) const +MappingQ1::compute_shapes (const std::vector > &unit_points, + InternalData &data) const { internal::MappingQ1:: - compute_shapes_virtual (n_shape_functions, - unit_points, data); + compute_shapes (n_shape_functions, + unit_points, data); } @@ -750,7 +730,7 @@ typename MappingQ1::InternalData * MappingQ1::get_data (const UpdateFlags update_flags, const Quadrature &q) const { - InternalData *data = new InternalData(1, n_shape_functions); + InternalData *data = new InternalData(1); data->initialize (requires_update_flags(update_flags), q, q.size()); compute_shapes (q.get_points(), *data); @@ -764,7 +744,7 @@ typename Mapping::InternalDataBase * MappingQ1::get_face_data (const UpdateFlags update_flags, const Quadrature &quadrature) const { - InternalData *data = new InternalData(1, n_shape_functions); + InternalData *data = new InternalData(1); data->initialize_face (requires_update_flags(update_flags), QProjector::project_to_all_faces(quadrature), quadrature.size()); @@ -781,7 +761,7 @@ typename Mapping::InternalDataBase * MappingQ1::get_subface_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { - InternalData *data = new InternalData(1, n_shape_functions); + InternalData *data = new InternalData(1); data->initialize_face (requires_update_flags(update_flags), QProjector::project_to_all_subfaces(quadrature), quadrature.size()); @@ -2533,7 +2513,11 @@ transform_real_to_unit_cell_internal Point p_unit = initial_p_unit; - compute_shapes(std::vector > (1, p_unit), mdata); + 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); + Point p_real = transform_unit_to_real_cell_internal(mdata); Tensor<1,spacedim> f = p_real-p; @@ -2620,7 +2604,11 @@ transform_real_to_unit_cell_internal // shape values and derivatives // at new p_unit point - compute_shapes(std::vector > (1, p_unit_trial), mdata); + 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); + // f(x) Point p_real_trial = transform_unit_to_real_cell_internal(mdata); @@ -2762,7 +2750,10 @@ transform_real_to_unit_cell_internal_codim1 Tensor<2,dim1> df; //Evaluate first and second derivatives - compute_shapes(std::vector > (1, p_unit), mdata); + 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); for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); @@ -2814,7 +2805,10 @@ transform_real_to_unit_cell_internal_codim1 D2F[j][l].clear(); } - compute_shapes(std::vector > (1, p_unit), mdata); + 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); for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); -- 2.39.5