From a332d1db2803d48c3720c45b63b0bd067138b5f5 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 2 Sep 2015 14:46:45 -0500 Subject: [PATCH] Attempt to disentangle MappingQ1 and MappingQ. Store the polynomial degree in MappingQ1::InternalData, rather than which kind of mapping created it. --- include/deal.II/fe/mapping_q.h | 7 +++- include/deal.II/fe/mapping_q1.h | 16 ++++---- source/fe/mapping_q.cc | 71 ++++++++++++++------------------- source/fe/mapping_q1.cc | 27 ++++++++----- 4 files changed, 60 insertions(+), 61 deletions(-) diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index ec29075344..317eb7953a 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -238,7 +238,7 @@ private: * @{ */ -protected: +public: /** * Storage for internal data of this mapping. See Mapping::InternalDataBase @@ -261,7 +261,8 @@ protected: /** * Constructor. */ - InternalData (const unsigned int n_shape_functions); + InternalData (const unsigned int polynomial_degree, + const unsigned int n_shape_functions); /** @@ -284,6 +285,8 @@ protected: typename MappingQ1::InternalData mapping_q1_data; }; +protected: + // documentation can be found in Mapping::get_data() virtual InternalData * diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h index 1b4c05bc77..082f5e7466 100644 --- a/include/deal.II/fe/mapping_q1.h +++ b/include/deal.II/fe/mapping_q1.h @@ -163,7 +163,8 @@ public: /** * Constructor. Pass the number of shape functions. */ - InternalData(const unsigned int n_shape_functions); + InternalData(const unsigned int polynomial_degree, + const unsigned int n_shape_functions); /** * Initialize the object's member variables related to cell data @@ -303,12 +304,11 @@ public: std::vector > > unit_tangentials; /** - * Default value of this flag is @p true. If *this is an object - * of a derived class, this flag is set to @p false. (This is, for example, - * the case for MappingQ, which derives MappingQ::InternalData from the - * current MappingQ1::InternalData.) + * The polynomial degree of the mapping. Since the objects here + * are also used (with minor adjustments) by MappingQ, we need to + * store this. */ - bool is_mapping_q1_data; + unsigned int polynomial_degree; /** * Number of shape functions. If this is a Q1 mapping, then it is simply @@ -466,8 +466,8 @@ protected: * 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.is_mapping_q1_data - * equals @p true or @p false. + * 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; diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 4f811729ed..08068cc01c 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -35,14 +35,15 @@ DEAL_II_NAMESPACE_OPEN template -MappingQ::InternalData::InternalData (const unsigned int n_shape_functions) +MappingQ::InternalData::InternalData (const unsigned int polynomial_degree, + const unsigned int n_shape_functions) : - MappingQ1::InternalData(n_shape_functions), + MappingQ1::InternalData(polynomial_degree, + n_shape_functions), use_mapping_q1_on_current_cell(false), - mapping_q1_data(1 << dim) -{ - this->is_mapping_q1_data=false; -} + mapping_q1_data(1, + GeometryInfo::vertices_per_cell) +{} @@ -226,7 +227,7 @@ typename MappingQ::InternalData * MappingQ::get_data (const UpdateFlags update_flags, const Quadrature &quadrature) const { - InternalData *data = new InternalData(n_shape_functions); + InternalData *data = new InternalData(degree, n_shape_functions); // fill the data of both the Q_p and the Q_1 objects in parallel Threads::TaskGroup<> tasks; @@ -261,7 +262,7 @@ typename Mapping::InternalDataBase * MappingQ::get_face_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { - InternalData *data = new InternalData(n_shape_functions); + InternalData *data = new InternalData(degree, n_shape_functions); const Quadrature q (QProjector::project_to_all_faces(quadrature)); // fill the data of both the Q_p and the Q_1 objects in parallel @@ -296,7 +297,7 @@ typename Mapping::InternalDataBase * MappingQ::get_subface_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { - InternalData *data = new InternalData(n_shape_functions); + InternalData *data = new InternalData(degree, n_shape_functions); const Quadrature q (QProjector::project_to_all_subfaces(quadrature)); // fill the data of both the Q_p and the Q_1 objects in parallel @@ -617,7 +618,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(n_shape_functions); + InternalData quadrature_data(degree, n_shape_functions); quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points); this->compute_shapes(quadrature.get_points(), quadrature_data); @@ -1009,15 +1010,13 @@ transform (const VectorSlice > > input, Assert(q1_data!=0, ExcInternalError()); // If it is a genuine MappingQ::InternalData, we have to test further - if (!q1_data->is_mapping_q1_data) + if (const InternalData *data = dynamic_cast(&mapping_data)) { - Assert (dynamic_cast(&mapping_data) != 0, - ExcInternalError()); - const InternalData &data = static_cast(mapping_data); // If we only use the Q1-portion, we have to extract that data object - if (data.use_mapping_q1_on_current_cell) - q1_data = &data.mapping_q1_data; + if (data->use_mapping_q1_on_current_cell) + q1_data = &data->mapping_q1_data; } + // Now, q1_data should have the right tensors in it and we call the base // classes transform function MappingQ1::transform(input, mapping_type, *q1_data, output); @@ -1041,15 +1040,13 @@ transform (const VectorSlice Assert(q1_data!=0, ExcInternalError()); // If it is a genuine MappingQ::InternalData, we have to test further - if (!q1_data->is_mapping_q1_data) + if (const InternalData *data = dynamic_cast(&mapping_data)) { - Assert (dynamic_cast(&mapping_data) != 0, - ExcInternalError()); - const InternalData &data = static_cast(mapping_data); // If we only use the Q1-portion, we have to extract that data object - if (data.use_mapping_q1_on_current_cell) - q1_data = &data.mapping_q1_data; + if (data->use_mapping_q1_on_current_cell) + q1_data = &data->mapping_q1_data; } + // Now, q1_data should have the right tensors in it and we call the base // classes transform function MappingQ1::transform(input, mapping_type, *q1_data, output); @@ -1071,15 +1068,13 @@ transform (const VectorSlice > > input, Assert(q1_data!=0, ExcInternalError()); // If it is a genuine MappingQ::InternalData, we have to test further - if (!q1_data->is_mapping_q1_data) + if (const InternalData *data = dynamic_cast(&mapping_data)) { - Assert (dynamic_cast(&mapping_data) != 0, - ExcInternalError()); - const InternalData &data = static_cast(mapping_data); // If we only use the Q1-portion, we have to extract that data object - if (data.use_mapping_q1_on_current_cell) - q1_data = &data.mapping_q1_data; + if (data->use_mapping_q1_on_current_cell) + q1_data = &data->mapping_q1_data; } + // Now, q1_data should have the right tensors in it and we call the base // classes transform function MappingQ1::transform(input, mapping_type, *q1_data, output); @@ -1103,15 +1098,13 @@ transform (const VectorSlice Assert(q1_data!=0, ExcInternalError()); // If it is a genuine MappingQ::InternalData, we have to test further - if (!q1_data->is_mapping_q1_data) + if (const InternalData *data = dynamic_cast(&mapping_data)) { - Assert (dynamic_cast(&mapping_data) != 0, - ExcInternalError()); - const InternalData &data = static_cast(mapping_data); // If we only use the Q1-portion, we have to extract that data object - if (data.use_mapping_q1_on_current_cell) - q1_data = &data.mapping_q1_data; + if (data->use_mapping_q1_on_current_cell) + q1_data = &data->mapping_q1_data; } + // Now, q1_data should have the right tensors in it and we call the base // classes transform function MappingQ1::transform(input, mapping_type, *q1_data, output); @@ -1133,15 +1126,13 @@ transform (const VectorSlice > > input, Assert(q1_data!=0, ExcInternalError()); // If it is a genuine MappingQ::InternalData, we have to test further - if (!q1_data->is_mapping_q1_data) + if (const InternalData *data = dynamic_cast(&mapping_data)) { - Assert (dynamic_cast(&mapping_data) != 0, - ExcInternalError()); - const InternalData &data = static_cast(mapping_data); // If we only use the Q1-portion, we have to extract that data object - if (data.use_mapping_q1_on_current_cell) - q1_data = &data.mapping_q1_data; + if (data->use_mapping_q1_on_current_cell) + q1_data = &data->mapping_q1_data; } + // Now, q1_data should have the right tensors in it and we call the base // classes transform function MappingQ1::transform(input, mapping_type, *q1_data, output); diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index 56e83e29fa..e6a176f48a 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -26,6 +26,7 @@ #include #include #include +#include #include #include @@ -41,11 +42,16 @@ const unsigned int MappingQ1::n_shape_functions; template -MappingQ1::InternalData::InternalData (const unsigned int n_shape_functions) +MappingQ1::InternalData::InternalData (const unsigned int polynomial_degree, + const unsigned int n_shape_functions) : - is_mapping_q1_data(true), + polynomial_degree (polynomial_degree), n_shape_functions (n_shape_functions) -{} +{ + Assert (n_shape_functions == + Utilities::fixed_power(polynomial_degree+1), + ExcInternalError()); +} @@ -63,7 +69,7 @@ MappingQ1::InternalData::memory_consumption () const MemoryConsumption::memory_consumption (mapping_support_points) + MemoryConsumption::memory_consumption (cell_of_current_support_points) + MemoryConsumption::memory_consumption (volume_elements) + - MemoryConsumption::memory_consumption (is_mapping_q1_data) + + MemoryConsumption::memory_consumption (polynomial_degree) + MemoryConsumption::memory_consumption (n_shape_functions)); } @@ -204,10 +210,9 @@ 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 (data.is_mapping_q1_data) + // 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); @@ -745,7 +750,7 @@ typename MappingQ1::InternalData * MappingQ1::get_data (const UpdateFlags update_flags, const Quadrature &q) const { - InternalData *data = new InternalData(n_shape_functions); + InternalData *data = new InternalData(1, n_shape_functions); data->initialize (requires_update_flags(update_flags), q, q.size()); compute_shapes (q.get_points(), *data); @@ -759,7 +764,7 @@ typename Mapping::InternalDataBase * MappingQ1::get_face_data (const UpdateFlags update_flags, const Quadrature &quadrature) const { - InternalData *data = new InternalData(n_shape_functions); + InternalData *data = new InternalData(1, n_shape_functions); data->initialize_face (requires_update_flags(update_flags), QProjector::project_to_all_faces(quadrature), quadrature.size()); @@ -776,7 +781,7 @@ typename Mapping::InternalDataBase * MappingQ1::get_subface_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { - InternalData *data = new InternalData(n_shape_functions); + InternalData *data = new InternalData(1, n_shape_functions); data->initialize_face (requires_update_flags(update_flags), QProjector::project_to_all_subfaces(quadrature), quadrature.size()); -- 2.39.5