From 18a76d3c55e2da46531e9a4906f285dfaf899bb8 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 21 Aug 2013 05:35:12 +0000 Subject: [PATCH] Cleanup some code so it compiles much quicker git-svn-id: https://svn.dealii.org/trunk@30364 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/hp/fe_values.h | 49 ++ .../deal.II/numerics/vector_tools.templates.h | 608 +++++++++--------- .../vector_tools_integrate_difference.inst.in | 72 +-- 3 files changed, 383 insertions(+), 346 deletions(-) diff --git a/deal.II/include/deal.II/hp/fe_values.h b/deal.II/include/deal.II/hp/fe_values.h index d077e4d838..cad63fdf4f 100644 --- a/deal.II/include/deal.II/hp/fe_values.h +++ b/deal.II/include/deal.II/hp/fe_values.h @@ -87,6 +87,25 @@ namespace internal const dealii::hp::FECollection & get_fe_collection () const; + /** + * Get a reference to the collection of mapping objects used + * here. + */ + const dealii::hp::MappingCollection & + get_mapping_collection () const; + + /** + * Get a reference to the collection of quadrature objects used + * here. + */ + const dealii::hp::QCollection & + get_quadrature_collection () const; + + /** + * Get the underlying update flags. + */ + UpdateFlags get_update_flags() const; + /** * Return a reference to the @p FEValues object selected by the last * call to select_fe_values(). select_fe_values() in turn is called when @@ -592,6 +611,36 @@ namespace internal { return *fe_collection; } + + + + template + inline + const dealii::hp::MappingCollection & + FEValuesBase::get_mapping_collection () const + { + return *mapping_collection; + } + + + + template + inline + const dealii::hp::QCollection & + FEValuesBase::get_quadrature_collection () const + { + return q_collection; + } + + + + template + inline + dealii::UpdateFlags + FEValuesBase::get_update_flags () const + { + return update_flags; + } } } diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index 6f9e2d5f3d..0c0799fc5d 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -4713,10 +4713,296 @@ namespace VectorTools namespace internal { + template + struct IDScratchData + { + IDScratchData (const dealii::hp::MappingCollection &mapping, + const dealii::hp::FECollection &fe, + const dealii::hp::QCollection &q, + const UpdateFlags update_flags); + + IDScratchData (const IDScratchData &data); + + void resize_vectors (const unsigned int n_q_points, + const unsigned int n_components); + + std::vector > function_values; + std::vector > > function_grads; + std::vector weight_values; + std::vector > weight_vectors; + + std::vector > psi_values; + std::vector > > psi_grads; + std::vector psi_scalar; + + std::vector tmp_values; + std::vector > tmp_gradients; + + dealii::hp::FEValues x_fe_values; + }; + + + template + IDScratchData + ::IDScratchData(const dealii::hp::MappingCollection &mapping, + const dealii::hp::FECollection &fe, + const dealii::hp::QCollection &q, + const UpdateFlags update_flags) + : + x_fe_values(mapping, fe, q, update_flags) + {} + + template + IDScratchData::IDScratchData (const IDScratchData &data) + : + x_fe_values(data.x_fe_values.get_mapping_collection(), + data.x_fe_values.get_fe_collection(), + data.x_fe_values.get_quadrature_collection(), + data.x_fe_values.get_update_flags()) + {} + + template + void + IDScratchData::resize_vectors (const unsigned int n_q_points, + const unsigned int n_components) + { + function_values.resize (n_q_points, + dealii::Vector(n_components)); + function_grads.resize (n_q_points, + std::vector >(n_components)); + + weight_values.resize (n_q_points); + weight_vectors.resize (n_q_points, + dealii::Vector(n_components)); + + psi_values.resize (n_q_points, + dealii::Vector(n_components)); + psi_grads.resize (n_q_points, + std::vector >(n_components)); + psi_scalar.resize (n_q_points); + + tmp_values.resize (n_q_points); + tmp_gradients.resize (n_q_points); + } + + + // avoid compiling inner function for many vector types when we always + // really do the same thing by putting the main work into this helper + // function + template + double + integrate_difference_inner (const Function &exact_solution, + const NormType &norm, + const Function *weight, + const UpdateFlags update_flags, + const double exponent, + const unsigned int n_components, + IDScratchData &data) + { + const bool fe_is_system = (n_components != 1); + const dealii::FEValues &fe_values = data.x_fe_values.get_present_fe_values (); + const unsigned int n_q_points = fe_values.n_quadrature_points; + + if (weight!=0) + { + if (weight->n_components>1) + weight->vector_value_list (fe_values.get_quadrature_points(), + data.weight_vectors); + else + { + weight->value_list (fe_values.get_quadrature_points(), + data.weight_values); + for (unsigned int k=0; k static void - do_integrate_difference (const dealii::hp::MappingCollection &mapping, + do_integrate_difference (const dealii::hp::MappingCollection &mapping, const DH &dof, const InVector &fe_function, const Function &exact_solution, @@ -4726,13 +5012,9 @@ namespace VectorTools const Function *weight, const double exponent_1) { - // we mark the "exponent" parameter - // to this function "const" since - // it is strictly incoming, but we - // need to set it to something - // different later on, if - // necessary, so have a read-write - // version of it: + // we mark the "exponent" parameter to this function "const" since it is + // strictly incoming, but we need to set it to something different later + // on, if necessary, so have a read-write version of it: double exponent = exponent_1; const unsigned int n_components = dof.get_fe().n_components(); @@ -4783,32 +5065,7 @@ namespace VectorTools } dealii::hp::FECollection fe_collection (dof.get_fe()); - dealii::hp::FEValues x_fe_values(mapping, fe_collection, q, update_flags); - - const unsigned int max_n_q_points = q.max_n_quadrature_points (); - - std::vector< dealii::Vector > - function_values (max_n_q_points, dealii::Vector(n_components)); - std::vector > > - function_grads (max_n_q_points, std::vector >(n_components)); - - std::vector - weight_values (max_n_q_points); - std::vector > - weight_vectors (max_n_q_points, dealii::Vector(n_components)); - - std::vector > - psi_values (max_n_q_points, dealii::Vector(n_components)); - std::vector > > - psi_grads (max_n_q_points, std::vector >(n_components)); - std::vector - psi_scalar (max_n_q_points); - - // tmp vector when we use the - // Function functions for - // scalar functions - std::vector tmp_values (max_n_q_points); - std::vector > tmp_gradients (max_n_q_points); + IDScratchData data(mapping, fe_collection, q, update_flags); // loop over all cells typename DH::active_cell_iterator cell = dof.begin_active(), @@ -4816,285 +5073,26 @@ namespace VectorTools for (unsigned int index=0; cell != endc; ++cell, ++index) if (cell->is_locally_owned()) { - double diff=0; // initialize for this cell - x_fe_values.reinit (cell); + data.x_fe_values.reinit (cell); - const dealii::FEValues &fe_values = x_fe_values.get_present_fe_values (); + const dealii::FEValues &fe_values = data.x_fe_values.get_present_fe_values (); const unsigned int n_q_points = fe_values.n_quadrature_points; - - // resize all out scratch - // arrays to the number of - // quadrature points we use - // for the present cell - function_values.resize (n_q_points, - dealii::Vector(n_components)); - function_grads.resize (n_q_points, - std::vector >(n_components)); - - weight_values.resize (n_q_points); - weight_vectors.resize (n_q_points, - dealii::Vector(n_components)); - - psi_values.resize (n_q_points, - dealii::Vector(n_components)); - psi_grads.resize (n_q_points, - std::vector >(n_components)); - psi_scalar.resize (n_q_points); - - tmp_values.resize (n_q_points); - tmp_gradients.resize (n_q_points); - - if (weight!=0) - { - if (weight->n_components>1) - weight->vector_value_list (fe_values.get_quadrature_points(), - weight_vectors); - else - { - weight->value_list (fe_values.get_quadrature_points(), - weight_values); - for (unsigned int k=0; k::mapping_collection, + ::do_integrate_difference(hp::StaticMappingQ1::mapping_collection, dof, fe_function, exact_solution, difference, q, norm, weight, exponent); diff --git a/deal.II/source/numerics/vector_tools_integrate_difference.inst.in b/deal.II/source/numerics/vector_tools_integrate_difference.inst.in index 08cf8b9d22..4ad4fe0ec4 100644 --- a/deal.II/source/numerics/vector_tools_integrate_difference.inst.in +++ b/deal.II/source/numerics/vector_tools_integrate_difference.inst.in @@ -22,7 +22,8 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens template void integrate_difference, deal_II_space_dimension> - (const DoFHandler&, + (const Mapping&, + const DoFHandler&, const VEC&, const Function&, Vector&, @@ -32,23 +33,23 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const double); template - void integrate_difference, deal_II_space_dimension > + void integrate_difference, deal_II_space_dimension> (const DoFHandler&, const VEC&, const Function&, - Vector&, + Vector&, const Quadrature&, const NormType&, const Function*, const double); template - void integrate_difference, deal_II_space_dimension> + void integrate_difference, deal_II_space_dimension > (const Mapping&, const DoFHandler&, const VEC&, const Function&, - Vector&, + Vector&, const Quadrature&, const NormType&, const Function*, @@ -56,8 +57,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens template void integrate_difference, deal_II_space_dimension > - (const Mapping&, - const DoFHandler&, + (const DoFHandler&, const VEC&, const Function&, Vector&, @@ -66,62 +66,52 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const Function*, const double); - \} -#endif - } - - - -//TODO[SP]: replace by -// where applicable and move to codimension cases above also when applicable -for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) - { -#if deal_II_dimension == deal_II_space_dimension - - namespace VectorTools \{ - template - void integrate_difference - (const hp::DoFHandler&, + void integrate_difference, deal_II_space_dimension> + (const hp::MappingCollection&, + const hp::DoFHandler&, const VEC&, - const Function&, + const Function&, Vector&, const hp::QCollection&, const NormType&, - const Function*, + const Function*, const double); + template - void integrate_difference - (const hp::DoFHandler&, + void integrate_difference, deal_II_space_dimension> + (const hp::DoFHandler&, const VEC&, - const Function&, - Vector&, + const Function&, + Vector&, const hp::QCollection&, const NormType&, - const Function*, + const Function*, const double); + template - void integrate_difference - (const hp::MappingCollection&, - const hp::DoFHandler&, + void integrate_difference, deal_II_space_dimension> + (const hp::MappingCollection&, + const hp::DoFHandler&, const VEC&, - const Function&, - Vector&, + const Function&, + Vector&, const hp::QCollection&, const NormType&, - const Function*, + const Function*, const double); + template - void integrate_difference - (const hp::MappingCollection&, - const hp::DoFHandler&, + void integrate_difference, deal_II_space_dimension> + (const hp::DoFHandler&, const VEC&, - const Function&, + const Function&, Vector&, const hp::QCollection&, const NormType&, - const Function*, + const Function*, const double); + \} #endif } -- 2.39.5