From 7c6bc8a1df7326937b9064c06e810bef2151bc84 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 23 Aug 2015 22:58:42 -0500 Subject: [PATCH] Move initialization of MappingQ1::InternalData into this class. This allows removing the poorly named 'compute_data()' function from the interface of MappingQ1 (don't all functions somehow compute some data?). --- include/deal.II/fe/mapping_q1.h | 39 +++-- source/fe/mapping_q.cc | 59 ++++---- source/fe/mapping_q1.cc | 246 ++++++++++++++++---------------- 3 files changed, 177 insertions(+), 167 deletions(-) diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h index c735f065ae..48b9d00ffb 100644 --- a/include/deal.II/fe/mapping_q1.h +++ b/include/deal.II/fe/mapping_q1.h @@ -165,6 +165,24 @@ public: */ InternalData(const unsigned int n_shape_functions); + /** + * Initialize the object's member variables related to cell data + * based on the given arguments. + */ + void + initialize (const UpdateFlags update_flags, + const Quadrature &quadrature, + const unsigned int n_original_q_points); + + /** + * Initialize the object's member variables related to cell and + * face data based on the given arguments. + */ + void + initialize_face (const UpdateFlags update_flags, + const Quadrature &quadrature, + const unsigned int n_original_q_points); + /** * Shape function at quadrature point. Shape functions are in tensor * product order, so vertices must be reordered to obtain transformation. @@ -414,27 +432,6 @@ protected: void compute_shapes (const std::vector > &unit_points, InternalData &data) const; - /** - * Do the computations for the @p get_data functions. Here, the data vectors - * of @p InternalData are reinitialized to proper size and shape values are - * computed. - */ - void compute_data (const UpdateFlags flags, - const Quadrature &quadrature, - const unsigned int n_orig_q_points, - InternalData &data) const; - - /** - * Do the computations for the @p get_face_data functions. Here, the data - * vectors of @p InternalData are reinitialized to proper size and shape - * values and derivatives are computed. Furthermore @p unit_tangential - * vectors of the face are computed. - */ - void compute_face_data (const UpdateFlags flags, - const Quadrature &quadrature, - const unsigned int n_orig_q_points, - InternalData &data) const; - /** * Compute shape values and/or derivatives. */ diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 29fa65d00e..ad5784c1de 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -203,22 +203,26 @@ MappingQ::get_data (const UpdateFlags update_flags, // fill the data of both the Q_p and the Q_1 objects in parallel Threads::TaskGroup<> tasks; tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ::compute_data, - this, + (std_cxx11::bind(&MappingQ1::InternalData::initialize, + data, update_flags, std_cxx11::cref(quadrature), - quadrature.size(), - std_cxx11::ref(*data)))); + quadrature.size()))); if (!use_mapping_q_on_all_cells) tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ::compute_data, - this, + (std_cxx11::bind(&MappingQ1::InternalData::initialize, + &data->mapping_q1_data, update_flags, std_cxx11::cref(quadrature), - quadrature.size(), - std_cxx11::ref(data->mapping_q1_data)))); + quadrature.size()))); tasks.join_all (); + // TODO: parallelize this as well + compute_shapes (quadrature.get_points(), *data); + if (!use_mapping_q_on_all_cells) + compute_shapes (quadrature.get_points(), data->mapping_q1_data); + + return data; } @@ -235,22 +239,25 @@ MappingQ::get_face_data (const UpdateFlags update_flags, // fill the data of both the Q_p and the Q_1 objects in parallel Threads::TaskGroup<> tasks; tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ::compute_face_data, - this, + (std_cxx11::bind(&MappingQ1::InternalData::initialize_face, + data, update_flags, std_cxx11::cref(q), - quadrature.size(), - std_cxx11::ref(*data)))); + quadrature.size()))); if (!use_mapping_q_on_all_cells) tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ::compute_face_data, - this, + (std_cxx11::bind(&MappingQ1::InternalData::initialize_face, + &data->mapping_q1_data, update_flags, std_cxx11::cref(q), - quadrature.size(), - std_cxx11::ref(data->mapping_q1_data)))); + quadrature.size()))); tasks.join_all (); + // TODO: parallelize this as well + compute_shapes (q.get_points(), *data); + if (!use_mapping_q_on_all_cells) + compute_shapes (q.get_points(), data->mapping_q1_data); + return data; } @@ -267,22 +274,26 @@ MappingQ::get_subface_data (const UpdateFlags update_flags, // fill the data of both the Q_p and the Q_1 objects in parallel Threads::TaskGroup<> tasks; tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ::compute_face_data, - this, + (std_cxx11::bind(&MappingQ1::InternalData::initialize_face, + data, update_flags, std_cxx11::cref(q), - quadrature.size(), - std_cxx11::ref(*data)))); + quadrature.size()))); if (!use_mapping_q_on_all_cells) tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ::compute_face_data, - this, + (std_cxx11::bind(&MappingQ1::InternalData::initialize_face, + &data->mapping_q1_data, update_flags, std_cxx11::cref(q), - quadrature.size(), - std_cxx11::ref(data->mapping_q1_data)))); + quadrature.size()))); tasks.join_all (); + // TODO: parallelize this as well + compute_shapes (q.get_points(), *data); + if (!use_mapping_q_on_all_cells) + compute_shapes (q.get_points(), data->mapping_q1_data); + + return data; } diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index b1e4ec2a6e..4d0d1a1eb7 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -68,6 +68,116 @@ MappingQ1::InternalData::memory_consumption () const } +template +void +MappingQ1::InternalData:: +initialize (const UpdateFlags update_flags, + const Quadrature &q, + const unsigned int n_original_q_points) +{ + // store the flags in the internal data object so we can access them + // in fill_fe_*_values() + this->update_each = update_flags; + + const unsigned int n_q_points = q.size(); + + // see if we need the (transformation) shape function values + // and/or gradients and resize the necessary arrays + if (this->update_each & update_quadrature_points) + shape_values.resize(n_shape_functions * n_q_points); + + if (this->update_each & (update_covariant_transformation + | update_contravariant_transformation + | update_JxW_values + | update_boundary_forms + | update_normal_vectors + | update_jacobians + | update_jacobian_grads + | update_inverse_jacobians)) + shape_derivatives.resize(n_shape_functions * n_q_points); + + if (this->update_each & update_covariant_transformation) + covariant.resize(n_original_q_points); + + if (this->update_each & update_contravariant_transformation) + contravariant.resize(n_original_q_points); + + if (this->update_each & update_volume_elements) + volume_elements.resize(n_original_q_points); + + if (this->update_each & update_jacobian_grads) + shape_second_derivatives.resize(n_shape_functions * n_q_points); +} + + +template +void +MappingQ1::InternalData:: +initialize_face (const UpdateFlags update_flags, + const Quadrature &q, + const unsigned int n_original_q_points) +{ + initialize (update_flags, q, n_original_q_points); + + if (dim > 1) + { + if (this->update_each & update_boundary_forms) + { + aux.resize (dim-1, std::vector > (n_original_q_points)); + + // Compute tangentials to the + // unit cell. + const unsigned int nfaces = GeometryInfo::faces_per_cell; + unit_tangentials.resize (nfaces*(dim-1), + std::vector > (n_original_q_points)); + if (dim==2) + { + // ensure a counterclockwise + // orientation of tangentials + static const int tangential_orientation[4]= {-1,1,1,-1}; + for (unsigned int i=0; i tang; + tang[1-i/2]=tangential_orientation[i]; + std::fill (unit_tangentials[i].begin(), + unit_tangentials[i].end(), tang); + } + } + else if (dim==3) + { + for (unsigned int i=0; i tang1, tang2; + + const unsigned int nd= + GeometryInfo::unit_normal_direction[i]; + + // first tangential + // vector in direction + // of the (nd+1)%3 axis + // and inverted in case + // of unit inward normal + tang1[(nd+1)%dim]=GeometryInfo::unit_normal_orientation[i]; + // second tangential + // vector in direction + // of the (nd+2)%3 axis + tang2[(nd+2)%dim]=1.; + + // same unit tangents + // for all quadrature + // points on this face + std::fill (unit_tangentials[i].begin(), + unit_tangentials[i].end(), tang1); + std::fill (unit_tangentials[nfaces+i].begin(), + unit_tangentials[nfaces+i].end(), tang2); + } + } + } + } +} + + + template MappingQ1::MappingQ1 () @@ -531,51 +641,6 @@ MappingQ1::requires_update_flags (const UpdateFlags in) const } -template -void -MappingQ1::compute_data (const UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points, - InternalData &data) const -{ - // store the flags in the internal data object so we can access them - // in fill_fe_*_values(). use the transitive hull of the required - // flags - data.update_each = requires_update_flags(update_flags); - - const unsigned int n_q_points = q.size(); - - // see if we need the (transformation) shape function values - // and/or gradients and resize the necessary arrays - if (data.update_each & update_quadrature_points) - data.shape_values.resize(data.n_shape_functions * n_q_points); - - if (data.update_each & (update_covariant_transformation - | update_contravariant_transformation - | update_JxW_values - | update_boundary_forms - | update_normal_vectors - | update_jacobians - | update_jacobian_grads - | update_inverse_jacobians)) - data.shape_derivatives.resize(data.n_shape_functions * n_q_points); - - if (data.update_each & update_covariant_transformation) - data.covariant.resize(n_original_q_points); - - if (data.update_each & update_contravariant_transformation) - data.contravariant.resize(n_original_q_points); - - if (data.update_each & update_volume_elements) - data.volume_elements.resize(n_original_q_points); - - if (data.update_each & update_jacobian_grads) - data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points); - - compute_shapes (q.get_points(), data); -} - - template typename MappingQ1::InternalData * @@ -583,76 +648,10 @@ MappingQ1::get_data (const UpdateFlags update_flags, const Quadrature &q) const { InternalData *data = new InternalData(n_shape_functions); - compute_data (update_flags, q, q.size(), *data); - return data; -} - - + data->initialize (requires_update_flags(update_flags), q, q.size()); + compute_shapes (q.get_points(), *data); -template -void -MappingQ1::compute_face_data (const UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points, - InternalData &data) const -{ - compute_data (update_flags, q, n_original_q_points, data); - - if (dim > 1) - { - if (data.update_each & update_boundary_forms) - { - data.aux.resize (dim-1, std::vector > (n_original_q_points)); - - // Compute tangentials to the - // unit cell. - const unsigned int nfaces = GeometryInfo::faces_per_cell; - data.unit_tangentials.resize (nfaces*(dim-1), - std::vector > (n_original_q_points)); - if (dim==2) - { - // ensure a counterclockwise - // orientation of tangentials - static const int tangential_orientation[4]= {-1,1,1,-1}; - for (unsigned int i=0; i tang; - tang[1-i/2]=tangential_orientation[i]; - std::fill (data.unit_tangentials[i].begin(), - data.unit_tangentials[i].end(), tang); - } - } - else if (dim==3) - { - for (unsigned int i=0; i tang1, tang2; - - const unsigned int nd= - GeometryInfo::unit_normal_direction[i]; - - // first tangential - // vector in direction - // of the (nd+1)%3 axis - // and inverted in case - // of unit inward normal - tang1[(nd+1)%dim]=GeometryInfo::unit_normal_orientation[i]; - // second tangential - // vector in direction - // of the (nd+2)%3 axis - tang2[(nd+2)%dim]=1.; - - // same unit tangents - // for all quadrature - // points on this face - std::fill (data.unit_tangentials[i].begin(), - data.unit_tangentials[i].end(), tang1); - std::fill (data.unit_tangentials[nfaces+i].begin(), - data.unit_tangentials[nfaces+i].end(), tang2); - } - } - } - } + return data; } @@ -663,10 +662,11 @@ MappingQ1::get_face_data (const UpdateFlags update_flags, const Quadrature &quadrature) const { InternalData *data = new InternalData(n_shape_functions); - compute_face_data (update_flags, - QProjector::project_to_all_faces(quadrature), - quadrature.size(), - *data); + 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); return data; } @@ -679,10 +679,12 @@ MappingQ1::get_subface_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { InternalData *data = new InternalData(n_shape_functions); - compute_face_data (update_flags, - QProjector::project_to_all_subfaces(quadrature), - quadrature.size(), - *data); + 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); + return data; } -- 2.39.5