From b9a02ce66992892dd840079a1cd987e1ccc9fe2c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 5 Apr 2015 13:32:44 -0500 Subject: [PATCH] A few minor cleanups. 1/ In integrate_difference(), use the correct type. 2/ In point_gradient(), return an object of the correct type. --- include/deal.II/numerics/vector_tools.h | 16 +- .../deal.II/numerics/vector_tools.templates.h | 156 ++++++++++-------- .../vector_tools_point_gradient.inst.in | 16 +- 3 files changed, 99 insertions(+), 89 deletions(-) diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 2da0dec20c..5d953427d5 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -2122,7 +2122,7 @@ namespace VectorTools point_gradient (const DoFHandler &dof, const InVector &fe_function, const Point &point, - std::vector > &value); + std::vector > &value); /** * Same as above for hp. @@ -2135,7 +2135,7 @@ namespace VectorTools point_gradient (const hp::DoFHandler &dof, const InVector &fe_function, const Point &point, - std::vector > &value); + std::vector > &value); /** * Evaluate a scalar finite element function defined by the given DoFHandler @@ -2149,7 +2149,7 @@ namespace VectorTools * exception of type VectorTools::ExcPointNotAvailableHere is thrown. */ template - Tensor<1, spacedim> + Tensor<1, spacedim, typename InVector::value_type> point_gradient (const DoFHandler &dof, const InVector &fe_function, const Point &point); @@ -2161,7 +2161,7 @@ namespace VectorTools * exception of type VectorTools::ExcPointNotAvailableHere is thrown. */ template - Tensor<1, spacedim> + Tensor<1, spacedim, typename InVector::value_type> point_gradient (const hp::DoFHandler &dof, const InVector &fe_function, const Point &point); @@ -2183,7 +2183,7 @@ namespace VectorTools const DoFHandler &dof, const InVector &fe_function, const Point &point, - std::vector > &value); + std::vector > &value); /** * Same as above for hp. @@ -2197,7 +2197,7 @@ namespace VectorTools const hp::DoFHandler &dof, const InVector &fe_function, const Point &point, - std::vector > &value); + std::vector > &value); /** * Evaluate a scalar finite element function defined by the given DoFHandler @@ -2211,7 +2211,7 @@ namespace VectorTools * exception of type VectorTools::ExcPointNotAvailableHere is thrown. */ template - Tensor<1, spacedim> + Tensor<1, spacedim, typename InVector::value_type> point_gradient (const Mapping &mapping, const DoFHandler &dof, const InVector &fe_function, @@ -2224,7 +2224,7 @@ namespace VectorTools * exception of type VectorTools::ExcPointNotAvailableHere is thrown. */ template - Tensor<1, spacedim> + Tensor<1, spacedim, typename InVector::value_type> point_gradient (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const InVector &fe_function, diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 23d525bde0..b78abbb44d 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -6027,7 +6027,7 @@ namespace VectorTools namespace internal { - template + template struct IDScratchData { IDScratchData (const dealii::hp::MappingCollection &mapping, @@ -6040,24 +6040,26 @@ namespace VectorTools void resize_vectors (const unsigned int n_q_points, const unsigned int n_components); - std::vector > function_values; - std::vector > > function_grads; + 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 > psi_values; + std::vector > > psi_grads; + std::vector psi_scalar; std::vector tmp_values; + std::vector > tmp_vector_values; std::vector > tmp_gradients; + std::vector > > tmp_vector_gradients; dealii::hp::FEValues x_fe_values; }; - template - IDScratchData + template + IDScratchData ::IDScratchData(const dealii::hp::MappingCollection &mapping, const dealii::hp::FECollection &fe, const dealii::hp::QCollection &q, @@ -6066,8 +6068,8 @@ namespace VectorTools x_fe_values(mapping, fe, q, update_flags) {} - template - IDScratchData::IDScratchData (const IDScratchData &data) + template + IDScratchData::IDScratchData (const IDScratchData &data) : x_fe_values(data.x_fe_values.get_mapping_collection(), data.x_fe_values.get_fe_collection(), @@ -6075,35 +6077,39 @@ namespace VectorTools data.x_fe_values.get_update_flags()) {} - template + template void - IDScratchData::resize_vectors (const unsigned int n_q_points, - const unsigned int n_components) + IDScratchData::resize_vectors (const unsigned int n_q_points, + const unsigned int n_components) { function_values.resize (n_q_points, - dealii::Vector(n_components)); + dealii::Vector(n_components)); function_grads.resize (n_q_points, - std::vector >(n_components)); + 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)); + dealii::Vector(n_components)); psi_grads.resize (n_q_points, - std::vector >(n_components)); + std::vector >(n_components)); psi_scalar.resize (n_q_points); tmp_values.resize (n_q_points); + tmp_vector_values.resize (n_q_points, + dealii::Vector(n_components)); tmp_gradients.resize (n_q_points); + tmp_vector_gradients.resize (n_q_points, + std::vector >(n_components)); } // 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 + template double integrate_difference_inner (const Function &exact_solution, const NormType &norm, @@ -6111,7 +6117,7 @@ namespace VectorTools const UpdateFlags update_flags, const double exponent, const unsigned int n_components, - IDScratchData &data) + IDScratchData &data) { const bool fe_is_system = (n_components != 1); const dealii::FEValues &fe_values = data.x_fe_values.get_present_fe_values (); @@ -6140,12 +6146,22 @@ namespace VectorTools if (update_flags & update_values) { // first compute the exact solution (vectors) at the quadrature - // points try to do this as efficient as possible by avoiding a + // points. try to do this as efficient as possible by avoiding a // second virtual function call in case the function really has only // one component + // + // TODO: we have to work a bit here because the Function + // interface of the argument denoting the exact function only + // provides us with double/Tensor<1,dim> values, rather than + // with the correct data type. so evaluate into a temp + // object, then copy around if (fe_is_system) - exact_solution.vector_value_list (fe_values.get_quadrature_points(), - data.psi_values); + { + exact_solution.vector_value_list (fe_values.get_quadrature_points(), + data.tmp_vector_values); + for (unsigned int i=0; i(fe_values.normal_vector(q)))* - fe_values.normal_vector(q)); + { + // compute (f.n) n + const Number f_dot_n + = (data.psi_grads[q][k] * Tensor<1,spacedim,Number>(fe_values.normal_vector(q))); + const Tensor<1,spacedim,Number> f_dot_n_times_n + = f_dot_n * Tensor<1,spacedim,Number>(fe_values.normal_vector(q)); + + data.psi_grads[q][k] -= (data.function_grads[q][k] + f_dot_n_times_n); + } else for (unsigned int k=0; k fe_collection (dof.get_fe()); - IDScratchData data(mapping, fe_collection, q, update_flags); - - // FIXME - // temporary vectors of consistent with InVector type. - // Need these because IDScratchData does not have a template type for the InVector - std::vector> function_values; - std::vector > > function_grads; + IDScratchData data(mapping, fe_collection, q, update_flags); // loop over all cells for (typename DH::active_cell_iterator cell = dof.begin_active(); @@ -6425,29 +6448,15 @@ namespace VectorTools const unsigned int n_q_points = fe_values.n_quadrature_points; data.resize_vectors (n_q_points, n_components); - //FIXME: - function_values.resize (n_q_points, - dealii::Vector(n_components)); - function_grads.resize (n_q_points, - std::vector >(n_components)); - if (update_flags & update_values) - fe_values.get_function_values (fe_function, function_values); + fe_values.get_function_values (fe_function, data.function_values); if (update_flags & update_gradients) - fe_values.get_function_gradients (fe_function, function_grads); - - // FIXME - for (unsigned int q = 0; q < n_q_points; q++) - for (unsigned int c = 0; c < n_components; c++) - { - data.function_values[q][c] = function_values[q][c]; - data.function_grads[q][c] = function_grads[q][c]; - } + fe_values.get_function_gradients (fe_function, data.function_grads); difference(cell->active_cell_index()) = - integrate_difference_inner (exact_solution, norm, weight, - update_flags, exponent, - n_components, data); + integrate_difference_inner (exact_solution, norm, weight, + update_flags, exponent, + n_components, data); } else // the cell is a ghost cell or is artificial. write a zero into the @@ -6781,7 +6790,7 @@ namespace VectorTools point_gradient (const DoFHandler &dof, const InVector &fe_function, const Point &point, - std::vector > &gradients) + std::vector > &gradients) { point_gradient (StaticMappingQ1::mapping, @@ -6797,7 +6806,7 @@ namespace VectorTools point_gradient (const hp::DoFHandler &dof, const InVector &fe_function, const Point &point, - std::vector > &gradients) + std::vector > &gradients) { point_gradient(hp::StaticMappingQ1::mapping_collection, dof, @@ -6808,7 +6817,7 @@ namespace VectorTools template - Tensor<1, spacedim> + Tensor<1, spacedim, typename InVector::value_type> point_gradient (const DoFHandler &dof, const InVector &fe_function, const Point &point) @@ -6821,7 +6830,7 @@ namespace VectorTools template - Tensor<1, spacedim> + Tensor<1, spacedim, typename InVector::value_type> point_gradient (const hp::DoFHandler &dof, const InVector &fe_function, const Point &point) @@ -6839,7 +6848,7 @@ namespace VectorTools const DoFHandler &dof, const InVector &fe_function, const Point &point, - std::vector > &gradient) + std::vector > &gradient) { const FiniteElement &fe = dof.get_fe(); @@ -6868,7 +6877,7 @@ namespace VectorTools // the given fe_function at this point typedef typename InVector::value_type Number; std::vector > > - u_gradient(1, std::vector > (fe.n_components())); + u_gradient(1, std::vector > (fe.n_components())); fe_values.get_function_gradients(fe_function, u_gradient); gradient = u_gradient[0]; @@ -6881,7 +6890,7 @@ namespace VectorTools const hp::DoFHandler &dof, const InVector &fe_function, const Point &point, - std::vector > &gradient) + std::vector > &gradient) { typedef typename InVector::value_type Number; const hp::FECollection &fe = dof.get_fe(); @@ -6910,7 +6919,7 @@ namespace VectorTools // the given fe_function at this point typedef typename InVector::value_type Number; std::vector > > - u_gradient(1, std::vector > (fe.n_components())); + u_gradient(1, std::vector > (fe.n_components())); fe_values.get_function_gradients(fe_function, u_gradient); gradient = u_gradient[0]; @@ -6918,7 +6927,7 @@ namespace VectorTools template - Tensor<1, spacedim> + Tensor<1, spacedim, typename InVector::value_type> point_gradient (const Mapping &mapping, const DoFHandler &dof, const InVector &fe_function, @@ -6927,15 +6936,16 @@ namespace VectorTools Assert(dof.get_fe().n_components() == 1, ExcMessage ("Finite element is not scalar as is necessary for this function")); - std::vector > gradient(1); + std::vector > gradient(1); point_gradient (mapping, dof, fe_function, point, gradient); return gradient[0]; } + template - Tensor<1, spacedim> + Tensor<1, spacedim, typename InVector::value_type> point_gradient (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const InVector &fe_function, @@ -6944,7 +6954,7 @@ namespace VectorTools Assert(dof.get_fe().n_components() == 1, ExcMessage ("Finite element is not scalar as is necessary for this function")); - std::vector > gradient(1); + std::vector > gradient(1); point_gradient (mapping, dof, fe_function, point, gradient); return gradient[0]; @@ -7003,7 +7013,7 @@ namespace VectorTools namespace internal { template - void set_possibly_complex_number(Number &n, const double &r, const double &i) + void set_possibly_complex_number(Number &n, const double &r, const double &) { n = r; } diff --git a/source/numerics/vector_tools_point_gradient.inst.in b/source/numerics/vector_tools_point_gradient.inst.in index 0844b7ce62..63417b4693 100644 --- a/source/numerics/vector_tools_point_gradient.inst.in +++ b/source/numerics/vector_tools_point_gradient.inst.in @@ -24,10 +24,10 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const hp::DoFHandler&, const VEC&, const Point&, - std::vector >&); + std::vector >&); template - Tensor<1,deal_II_space_dimension> point_gradient ( + Tensor<1,deal_II_space_dimension,VEC::value_type> point_gradient ( const hp::DoFHandler&, const VEC&, const Point&); @@ -38,10 +38,10 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const hp::DoFHandler&, const VEC&, const Point&, - std::vector >&); + std::vector >&); template - Tensor<1,deal_II_space_dimension> point_gradient ( + Tensor<1,deal_II_space_dimension,VEC::value_type> point_gradient ( const hp::MappingCollection&, const hp::DoFHandler&, const VEC&, @@ -52,10 +52,10 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const DoFHandler&, const VEC&, const Point&, - std::vector >&); + std::vector >&); template - Tensor<1,deal_II_space_dimension> point_gradient ( + Tensor<1,deal_II_space_dimension,VEC::value_type> point_gradient ( const DoFHandler&, const VEC&, const Point&); @@ -66,10 +66,10 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const DoFHandler&, const VEC&, const Point&, - std::vector >&); + std::vector >&); template - Tensor<1,deal_II_space_dimension> point_gradient ( + Tensor<1,deal_II_space_dimension,VEC::value_type> point_gradient ( const Mapping&, const DoFHandler&, const VEC&, -- 2.39.5