From: Daniel Garcia-Sanchez Date: Wed, 9 Oct 2019 07:43:41 +0000 (+0200) Subject: Add std::complex support to integrate_difference() X-Git-Tag: v9.2.0-rc1~984^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=57e6ebf571633f7e44a5fb6e64831d494f9bf79d;p=dealii.git Add std::complex support to integrate_difference() The type of exact_solution was double and fe_function could have any numerical type (double, float, std::complex or std::complex). This would lead in certain cases to unnecessary casts, for example from double to float. In addition it was not possible to compare a std::complex exact_solution to a std::complex fe_function. Now the types of exact_solution and fe_function must be the same, as a result integrate_difference() can be used to compare a std::complex exact_solution to a std::complex fe_function. This change is not backward-compatible. The old version of the function has been deprecated. --- diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 986e3125ef..52ae5023a1 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -2516,6 +2516,82 @@ namespace VectorTools */ template void + integrate_difference( + const Mapping & mapping, + const DoFHandler & dof, + const InVector & fe_function, + const Function &exact_solution, + OutVector & difference, + const Quadrature & q, + const NormType & norm, + const Function * weight = nullptr, + const double exponent = 2.); + + /** + * Call the integrate_difference() function, see above, with + * mapping=MappingQGeneric@(1). + */ + template + void + integrate_difference( + const DoFHandler & dof, + const InVector & fe_function, + const Function &exact_solution, + OutVector & difference, + const Quadrature & q, + const NormType & norm, + const Function * weight = nullptr, + const double exponent = 2.); + + /** + * Same as above for hp. + */ + template + void + integrate_difference( + const hp::MappingCollection & mapping, + const hp::DoFHandler & dof, + const InVector & fe_function, + const Function &exact_solution, + OutVector & difference, + const hp::QCollection & q, + const NormType & norm, + const Function * weight = nullptr, + const double exponent = 2.); + + /** + * Call the integrate_difference() function, see above, with + * mapping=MappingQGeneric@(1). + */ + template + void + integrate_difference( + const hp::DoFHandler & dof, + const InVector & fe_function, + const Function &exact_solution, + OutVector & difference, + const hp::QCollection & q, + const NormType & norm, + const Function * weight = nullptr, + const double exponent = 2.); + + /** + * Compute the cellwise error of the finite element solution. Integrate the + * difference between a reference function which is given as a continuous + * function object, and a finite element function. The result of this + * function is the vector @p difference that contains one value per active + * cell $K$ of the triangulation. Each of the values of this vector $d$ + * equals + * @f{align*}{ + * d_K = \| u-u_h \|_X + * @f} + * where $X$ denotes the norm chosen and $u$ represents the exact solution. + * + * @deprecated Use integrate_difference(const Mapping &, const DoFHandler &, const InVector &, const Function &, OutVector &, const Quadrature &, const NormType &, const Function *, const double) instead. + */ + template + DEAL_II_DEPRECATED typename std::enable_if< + !std::is_same::value>::type integrate_difference(const Mapping & mapping, const DoFHandler & dof, const InVector & fe_function, @@ -2529,9 +2605,12 @@ namespace VectorTools /** * Call the integrate_difference() function, see above, with * mapping=MappingQGeneric@(1). + * + * @deprecated Use integrate_difference(const DoFHandler &, const InVector &, const Function &exact_solution, OutVector &, const Quadrature &, const NormType &, const Function *, const double) instead. */ template - void + DEAL_II_DEPRECATED typename std::enable_if< + !std::is_same::value>::type integrate_difference(const DoFHandler & dof, const InVector & fe_function, const Function &exact_solution, @@ -2543,9 +2622,12 @@ namespace VectorTools /** * Same as above for hp. + * + * @deprecated Use integrate_difference(const hp::MappingCollection &, const hp::DoFHandler &, const InVector &, const Function &, OutVector &, const hp::QCollection &, const NormType &, const Function *, const double) instead. */ template - void + DEAL_II_DEPRECATED typename std::enable_if< + !std::is_same::value>::type integrate_difference(const hp::MappingCollection &mapping, const hp::DoFHandler & dof, const InVector & fe_function, @@ -2559,9 +2641,12 @@ namespace VectorTools /** * Call the integrate_difference() function, see above, with * mapping=MappingQGeneric@(1). + * + * @deprecated Use integrate_difference(const hp::DoFHandler &, const InVector &, const Function &, OutVector &, const hp::QCollection &, const NormType &, const Function *, const double) instead. */ template - void + DEAL_II_DEPRECATED typename std::enable_if< + !std::is_same::value>::type integrate_difference(const hp::DoFHandler &dof, const InVector & fe_function, const Function & exact_solution, diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 61a9b3588e..35b085efd8 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -7819,10 +7819,11 @@ namespace VectorTools 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; + std::vector tmp_values; + std::vector> tmp_vector_values; + std::vector> tmp_gradients; + std::vector>> + tmp_vector_gradients; dealii::hp::FEValues x_fe_values; }; @@ -7859,55 +7860,420 @@ namespace VectorTools weight_values.resize(n_q_points); weight_vectors.resize(n_q_points, Vector(n_components)); - psi_values.resize(n_q_points, Vector(n_components)); - psi_grads.resize(n_q_points, - std::vector>(n_components)); - psi_scalar.resize(n_q_points); + psi_values.resize(n_q_points, 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_vector_values.resize(n_q_points, Vector(n_components)); + tmp_gradients.resize(n_q_points); + tmp_vector_gradients.resize( + n_q_points, std::vector>(n_components)); + } + + template + struct DEAL_II_DEPRECATED DeprecatedIDScratchData + { + DeprecatedIDScratchData( + const dealii::hp::MappingCollection &mapping, + const dealii::hp::FECollection & fe, + const dealii::hp::QCollection & q, + const UpdateFlags update_flags); + + DeprecatedIDScratchData(const DeprecatedIDScratchData &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_vector_values; + std::vector> tmp_gradients; + std::vector>> tmp_vector_gradients; + + dealii::hp::FEValues x_fe_values; + }; + + + template + DeprecatedIDScratchData::DeprecatedIDScratchData( + 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 + DeprecatedIDScratchData::DeprecatedIDScratchData( + const DeprecatedIDScratchData &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 + DeprecatedIDScratchData::resize_vectors( + const unsigned int n_q_points, + const unsigned int n_components) + { + function_values.resize(n_q_points, 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, Vector(n_components)); + + psi_values.resize(n_q_points, 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_vector_values.resize(n_q_points, Vector(n_components)); + tmp_gradients.resize(n_q_points); + tmp_vector_gradients.resize( + n_q_points, std::vector>(n_components)); + } + + namespace internal + { + template + double + mean_to_double(const number &mean_value) + { + return mean_value; + } + + template + double + mean_to_double(const std::complex &mean_value) + { + // we need to return double as a norm, but mean value is a complex + // number. Panic and return real-part while warning the user that + // he shall never do that. + Assert( + false, + ExcMessage( + "Mean value norm is not implemented for complex-valued vectors")); + return mean_value.real(); + } + } // namespace internal + + + // 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 != nullptr) + { + 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 < n_q_points; ++k) + data.weight_vectors[k] = data.weight_values[k]; + } + } + else + { + for (unsigned int k = 0; k < n_q_points; ++k) + data.weight_vectors[k] = 1.; + } + + + 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 + // 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.tmp_vector_values); + for (unsigned int i = 0; i < n_q_points; ++i) + data.psi_values[i] = data.tmp_vector_values[i]; + } + else + { + exact_solution.value_list(fe_values.get_quadrature_points(), + data.tmp_values); + for (unsigned int i = 0; i < n_q_points; ++i) + data.psi_values[i](0) = data.tmp_values[i]; + } + + // then subtract finite element fe_function + for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int i = 0; i < data.psi_values[q].size(); ++i) + data.psi_values[q][i] -= data.function_values[q][i]; + } + + // Do the same for gradients, if required + if (update_flags & update_gradients) + { + // try to be a little clever to avoid recursive virtual function + // calls when calling gradient_list for functions that are really + // scalar functions + if (fe_is_system) + { + exact_solution.vector_gradient_list( + fe_values.get_quadrature_points(), data.tmp_vector_gradients); + for (unsigned int i = 0; i < n_q_points; ++i) + for (unsigned int comp = 0; comp < data.psi_grads[i].size(); + ++comp) + data.psi_grads[i][comp] = data.tmp_vector_gradients[i][comp]; + } + else + { + exact_solution.gradient_list(fe_values.get_quadrature_points(), + data.tmp_gradients); + for (unsigned int i = 0; i < n_q_points; ++i) + data.psi_grads[i][0] = data.tmp_gradients[i]; + } + + // then subtract finite element function_grads. We need to be + // careful in the codimension one case, since there we only have + // tangential gradients in the finite element function, not the full + // gradient. This is taken care of, by subtracting the normal + // component of the gradient from the exact function. + if (update_flags & update_normal_vectors) + for (unsigned int k = 0; k < n_components; ++k) + for (unsigned int q = 0; q < n_q_points; ++q) + { + // compute (f.n) n + const typename ProductType::type f_dot_n = + data.psi_grads[q][k] * fe_values.normal_vector(q); + const Tensor<1, spacedim, Number> f_dot_n_times_n = + f_dot_n * 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 < n_components; ++k) + for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int d = 0; d < spacedim; ++d) + data.psi_grads[q][k][d] -= data.function_grads[q][k][d]; + } + + double diff = 0; + Number diff_mean = 0; + + // First work on function values: + switch (norm) + { + case mean: + // Compute values in quadrature points and integrate + for (unsigned int q = 0; q < n_q_points; ++q) + { + Number sum = 0; + for (unsigned int k = 0; k < n_components; ++k) + if (data.weight_vectors[q](k) != 0) + sum += data.psi_values[q](k) * data.weight_vectors[q](k); + diff_mean += sum * fe_values.JxW(q); + } + break; + + case Lp_norm: + case L1_norm: + case W1p_norm: + // Compute values in quadrature points and integrate + for (unsigned int q = 0; q < n_q_points; ++q) + { + double sum = 0; + for (unsigned int k = 0; k < n_components; ++k) + if (data.weight_vectors[q](k) != 0) + sum += std::pow(static_cast( + numbers::NumberTraits::abs_square( + data.psi_values[q](k))), + exponent / 2.) * + data.weight_vectors[q](k); + diff += sum * fe_values.JxW(q); + } + + // Compute the root only if no derivative values are added later + if (!(update_flags & update_gradients)) + diff = std::pow(diff, 1. / exponent); + break; + + case L2_norm: + case H1_norm: + // Compute values in quadrature points and integrate + for (unsigned int q = 0; q < n_q_points; ++q) + { + double sum = 0; + for (unsigned int k = 0; k < n_components; ++k) + if (data.weight_vectors[q](k) != 0) + sum += numbers::NumberTraits::abs_square( + data.psi_values[q](k)) * + data.weight_vectors[q](k); + diff += sum * fe_values.JxW(q); + } + // Compute the root only, if no derivative values are added later + if (norm == L2_norm) + diff = std::sqrt(diff); + break; + + case Linfty_norm: + case W1infty_norm: + for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int k = 0; k < n_components; ++k) + if (data.weight_vectors[q](k) != 0) + diff = std::max(diff, + double(std::abs(data.psi_values[q](k) * + data.weight_vectors[q](k)))); + break; + + case H1_seminorm: + case Hdiv_seminorm: + case W1p_seminorm: + case W1infty_seminorm: + // function values are not used for these norms + break; + + default: + Assert(false, ExcNotImplemented()); + break; + } + + // Now compute terms depending on derivatives: + switch (norm) + { + case W1p_seminorm: + case W1p_norm: + for (unsigned int q = 0; q < n_q_points; ++q) + { + double sum = 0; + for (unsigned int k = 0; k < n_components; ++k) + if (data.weight_vectors[q](k) != 0) + sum += std::pow(data.psi_grads[q][k].norm_square(), + exponent / 2.) * + data.weight_vectors[q](k); + diff += sum * fe_values.JxW(q); + } + diff = std::pow(diff, 1. / exponent); + break; + + case H1_seminorm: + case H1_norm: + for (unsigned int q = 0; q < n_q_points; ++q) + { + double sum = 0; + for (unsigned int k = 0; k < n_components; ++k) + if (data.weight_vectors[q](k) != 0) + sum += data.psi_grads[q][k].norm_square() * + data.weight_vectors[q](k); + diff += sum * fe_values.JxW(q); + } + diff = std::sqrt(diff); + break; + + case Hdiv_seminorm: + for (unsigned int q = 0; q < n_q_points; ++q) + { + unsigned int idx = 0; + if (weight != nullptr) + for (; idx < n_components; ++idx) + if (data.weight_vectors[0](idx) > 0) + break; + + Assert( + n_components >= idx + dim, + ExcMessage( + "You can only ask for the Hdiv norm for a finite element " + "with at least 'dim' components. In that case, this function " + "will find the index of the first non-zero weight and take " + "the divergence of the 'dim' components that follow it.")); + + Number sum = 0; + // take the trace of the derivatives scaled by the weight and + // square it + for (unsigned int k = idx; k < idx + dim; ++k) + if (data.weight_vectors[q](k) != 0) + sum += data.psi_grads[q][k][k - idx] * + std::sqrt(data.weight_vectors[q](k)); + diff += numbers::NumberTraits::abs_square(sum) * + fe_values.JxW(q); + } + diff = std::sqrt(diff); + break; - tmp_values.resize(n_q_points); - tmp_vector_values.resize(n_q_points, Vector(n_components)); - tmp_gradients.resize(n_q_points); - tmp_vector_gradients.resize( - n_q_points, std::vector>(n_components)); - } + case W1infty_seminorm: + case W1infty_norm: + { + double t = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int k = 0; k < n_components; ++k) + if (data.weight_vectors[q](k) != 0) + for (unsigned int d = 0; d < dim; ++d) + t = std::max(t, + double(std::abs(data.psi_grads[q][k][d]) * + data.weight_vectors[q](k))); - namespace internal - { - template - double - mean_to_double(const number &mean_value) - { - return mean_value; - } + // then add seminorm to norm if that had previously been computed + diff += t; + } + break; + default: + break; + } - template - double - mean_to_double(const std::complex &mean_value) - { - // we need to return double as a norm, but mean value is a complex - // number. Panic and return real-part while warning the user that - // he shall never do that. - Assert( - false, - ExcMessage( - "Mean value norm is not implemented for complex-valued vectors")); - return mean_value.real(); - } - } // namespace internal + if (norm == mean) + diff = internal::mean_to_double(diff_mean); + // append result of this cell to the end of the vector + AssertIsFinite(diff); + return diff; + } - // 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) + DEAL_II_DEPRECATED + typename std::enable_if::value, + double>::type + 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, + DeprecatedIDScratchData &data) { const bool fe_is_system = (n_components != 1); const dealii::FEValues &fe_values = @@ -8193,6 +8559,131 @@ namespace VectorTools typename DoFHandlerType, int spacedim> static void + do_integrate_difference( + const dealii::hp::MappingCollection & mapping, + const DoFHandlerType & dof, + const InVector & fe_function, + const Function &exact_solution, + OutVector & difference, + const dealii::hp::QCollection & q, + const NormType & norm, + const Function * weight, + const double exponent_1) + { + using Number = typename InVector::value_type; + // 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(0).n_components(); + + Assert(exact_solution.n_components == n_components, + ExcDimensionMismatch(exact_solution.n_components, n_components)); + + if (weight != nullptr) + { + Assert((weight->n_components == 1) || + (weight->n_components == n_components), + ExcDimensionMismatch(weight->n_components, n_components)); + } + + difference.reinit(dof.get_triangulation().n_active_cells()); + + switch (norm) + { + case L2_norm: + case H1_seminorm: + case H1_norm: + case Hdiv_seminorm: + exponent = 2.; + break; + + case L1_norm: + exponent = 1.; + break; + + default: + break; + } + + UpdateFlags update_flags = + UpdateFlags(update_quadrature_points | update_JxW_values); + switch (norm) + { + case H1_seminorm: + case Hdiv_seminorm: + case W1p_seminorm: + case W1infty_seminorm: + update_flags |= UpdateFlags(update_gradients); + if (spacedim == dim + 1) + update_flags |= UpdateFlags(update_normal_vectors); + + break; + + case H1_norm: + case W1p_norm: + case W1infty_norm: + update_flags |= UpdateFlags(update_gradients); + if (spacedim == dim + 1) + update_flags |= UpdateFlags(update_normal_vectors); + DEAL_II_FALLTHROUGH; + + default: + update_flags |= UpdateFlags(update_values); + break; + } + + const dealii::hp::FECollection &fe_collection = + dof.get_fe_collection(); + IDScratchData data(mapping, + fe_collection, + q, + update_flags); + + // loop over all cells + for (typename DoFHandlerType::active_cell_iterator cell = + dof.begin_active(); + cell != dof.end(); + ++cell) + if (cell->is_locally_owned()) + { + // initialize for this cell + data.x_fe_values.reinit(cell); + + const dealii::FEValues &fe_values = + data.x_fe_values.get_present_fe_values(); + const unsigned int n_q_points = fe_values.n_quadrature_points; + data.resize_vectors(n_q_points, n_components); + + if (update_flags & update_values) + fe_values.get_function_values(fe_function, data.function_values); + if (update_flags & update_gradients) + 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); + } + else + // the cell is a ghost cell or is artificial. write a zero into the + // corresponding value of the returned vector + difference(cell->active_cell_index()) = 0; + } + + template + DEAL_II_DEPRECATED static typename std::enable_if< + !std::is_same::value>::type do_integrate_difference( const dealii::hp::MappingCollection &mapping, const DoFHandlerType & dof, @@ -8270,10 +8761,10 @@ namespace VectorTools const dealii::hp::FECollection &fe_collection = dof.get_fe_collection(); - IDScratchData data(mapping, - fe_collection, - q, - update_flags); + DeprecatedIDScratchData data(mapping, + fe_collection, + q, + update_flags); // loop over all cells for (typename DoFHandlerType::active_cell_iterator cell = @@ -8317,6 +8808,32 @@ namespace VectorTools template void + integrate_difference( + const Mapping & mapping, + const DoFHandler & dof, + const InVector & fe_function, + const Function &exact_solution, + OutVector & difference, + const Quadrature & q, + const NormType & norm, + const Function * weight, + const double exponent) + { + internal ::do_integrate_difference(hp::MappingCollection( + mapping), + dof, + fe_function, + exact_solution, + difference, + hp::QCollection(q), + norm, + weight, + exponent); + } + + template + DEAL_II_DEPRECATED typename std::enable_if< + !std::is_same::value>::type integrate_difference(const Mapping & mapping, const DoFHandler &dof, const InVector & fe_function, @@ -8342,6 +8859,32 @@ namespace VectorTools template void + integrate_difference( + const DoFHandler & dof, + const InVector & fe_function, + const Function &exact_solution, + OutVector & difference, + const Quadrature & q, + const NormType & norm, + const Function * weight, + const double exponent) + { + internal ::do_integrate_difference( + hp::StaticMappingQ1::mapping_collection, + dof, + fe_function, + exact_solution, + difference, + hp::QCollection(q), + norm, + weight, + exponent); + } + + + template + DEAL_II_DEPRECATED typename std::enable_if< + !std::is_same::value>::type integrate_difference(const DoFHandler &dof, const InVector & fe_function, const Function & exact_solution, @@ -8367,6 +8910,32 @@ namespace VectorTools template void + integrate_difference( + const dealii::hp::MappingCollection & mapping, + const dealii::hp::DoFHandler & dof, + const InVector & fe_function, + const Function &exact_solution, + OutVector & difference, + const dealii::hp::QCollection & q, + const NormType & norm, + const Function * weight, + const double exponent) + { + internal ::do_integrate_difference(hp::MappingCollection( + mapping), + dof, + fe_function, + exact_solution, + difference, + q, + norm, + weight, + exponent); + } + + template + DEAL_II_DEPRECATED typename std::enable_if< + !std::is_same::value>::type integrate_difference( const dealii::hp::MappingCollection &mapping, const dealii::hp::DoFHandler & dof, @@ -8393,6 +8962,31 @@ namespace VectorTools template void + integrate_difference( + const dealii::hp::DoFHandler & dof, + const InVector & fe_function, + const Function &exact_solution, + OutVector & difference, + const dealii::hp::QCollection & q, + const NormType & norm, + const Function * weight, + const double exponent) + { + internal ::do_integrate_difference( + hp::StaticMappingQ1::mapping_collection, + dof, + fe_function, + exact_solution, + difference, + q, + norm, + weight, + exponent); + } + + template + DEAL_II_DEPRECATED typename std::enable_if< + !std::is_same::value>::type integrate_difference(const dealii::hp::DoFHandler &dof, const InVector & fe_function, const Function & exact_solution, diff --git a/source/numerics/vector_tools_integrate_difference.inst.in b/source/numerics/vector_tools_integrate_difference.inst.in index 6b05ae0d7b..62b17d2197 100644 --- a/source/numerics/vector_tools_integrate_difference.inst.in +++ b/source/numerics/vector_tools_integrate_difference.inst.in @@ -29,7 +29,7 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; const Mapping &, const DoFHandler &, const VEC &, - const Function &, + const Function &, Vector &, const Quadrature &, const NormType &, @@ -43,7 +43,7 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; deal_II_space_dimension>( const DoFHandler &, const VEC &, - const Function &, + const Function &, Vector &, const Quadrature &, const NormType &, @@ -58,7 +58,7 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; const Mapping &, const DoFHandler &, const VEC &, - const Function &, + const Function &, Vector &, const Quadrature &, const NormType &, @@ -72,7 +72,7 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; deal_II_space_dimension>( const DoFHandler &, const VEC &, - const Function &, + const Function &, Vector &, const Quadrature &, const NormType &, @@ -88,7 +88,7 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; &, const hp::DoFHandler &, const VEC &, - const Function &, + const Function &, Vector &, const hp::QCollection &, const NormType &, @@ -102,7 +102,7 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; deal_II_space_dimension>( const hp::DoFHandler &, const VEC &, - const Function &, + const Function &, Vector &, const hp::QCollection &, const NormType &, @@ -118,7 +118,7 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; &, const hp::DoFHandler &, const VEC &, - const Function &, + const Function &, Vector &, const hp::QCollection &, const NormType &, @@ -132,12 +132,660 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; deal_II_space_dimension>( const hp::DoFHandler &, const VEC &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + \} +#endif + } + +// All the functions in this for loop are deprecated +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + namespace VectorTools + \{ + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const Mapping &, + const DoFHandler &, + const Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const Mapping &, + const DoFHandler &, + const BlockVector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const Mapping &, + const DoFHandler &, + const LinearAlgebra::Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const Mapping &, + const DoFHandler &, + const LinearAlgebra::distributed::Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const Mapping &, + const DoFHandler &, + const LinearAlgebra::distributed::BlockVector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + +# ifdef DEAL_II_TRILINOS_WITH_TPETRA + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const Mapping &, + const DoFHandler &, + const LinearAlgebra::TpetraWrappers::Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); +# endif + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const DoFHandler &, + const Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const DoFHandler &, + const BlockVector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const DoFHandler &, + const LinearAlgebra::Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const DoFHandler &, + const LinearAlgebra::distributed::Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const DoFHandler &, + const LinearAlgebra::distributed::BlockVector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + +# ifdef DEAL_II_TRILINOS_WITH_TPETRA + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const DoFHandler &, + const LinearAlgebra::TpetraWrappers::Vector, + &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); +# endif + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const Mapping &, + const DoFHandler &, + const BlockVector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const Mapping &, + const DoFHandler &, + const LinearAlgebra::Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const Mapping &, + const DoFHandler &, + const LinearAlgebra::distributed::Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const Mapping &, + const DoFHandler &, + const LinearAlgebra::distributed::BlockVector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + +# ifdef DEAL_II_TRILINOS_WITH_TPETRA + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const Mapping &, + const DoFHandler &, + const LinearAlgebra::TpetraWrappers::Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); +# endif + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const DoFHandler &, + const BlockVector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const DoFHandler &, + const LinearAlgebra::Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const DoFHandler &, + const LinearAlgebra::distributed::Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const DoFHandler &, + const LinearAlgebra::distributed::BlockVector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); + +# ifdef DEAL_II_TRILINOS_WITH_TPETRA + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const DoFHandler &, + const LinearAlgebra::TpetraWrappers::Vector &, + const Function &, + Vector &, + const Quadrature &, + const NormType &, + const Function *, + const double); +# endif + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::MappingCollection + &, + const hp::DoFHandler &, + const BlockVector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::MappingCollection + &, + const hp::DoFHandler &, + const LinearAlgebra::Vector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::MappingCollection + &, + const hp::DoFHandler &, + const LinearAlgebra::distributed::Vector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::MappingCollection + &, + const hp::DoFHandler &, + const LinearAlgebra::distributed::BlockVector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + +# ifdef DEAL_II_TRILINOS_WITH_TPETRA + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::MappingCollection + &, + const hp::DoFHandler &, + const LinearAlgebra::TpetraWrappers::Vector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); +# endif + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::DoFHandler &, + const BlockVector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::DoFHandler &, + const LinearAlgebra::Vector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::DoFHandler &, + const LinearAlgebra::distributed::Vector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::DoFHandler &, + const LinearAlgebra::distributed::BlockVector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + +# ifdef DEAL_II_TRILINOS_WITH_TPETRA + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::DoFHandler &, + const LinearAlgebra::TpetraWrappers::Vector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); +# endif + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::MappingCollection + &, + const hp::DoFHandler &, + const BlockVector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::MappingCollection + &, + const hp::DoFHandler &, + const LinearAlgebra::Vector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::MappingCollection + &, + const hp::DoFHandler &, + const LinearAlgebra::distributed::Vector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::MappingCollection + &, + const hp::DoFHandler &, + const LinearAlgebra::distributed::BlockVector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + +# ifdef DEAL_II_TRILINOS_WITH_TPETRA + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::MappingCollection + &, + const hp::DoFHandler &, + const LinearAlgebra::TpetraWrappers::Vector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); +# endif + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::DoFHandler &, + const BlockVector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::DoFHandler &, + const LinearAlgebra::Vector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::DoFHandler &, + const LinearAlgebra::distributed::Vector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::DoFHandler &, + const LinearAlgebra::distributed::BlockVector &, + const Function &, + Vector &, + const hp::QCollection &, + const NormType &, + const Function *, + const double); + +# ifdef DEAL_II_TRILINOS_WITH_TPETRA + template void + integrate_difference, + Vector, + deal_II_space_dimension>( + const hp::DoFHandler &, + const LinearAlgebra::TpetraWrappers::Vector &, const Function &, Vector &, const hp::QCollection &, const NormType &, const Function *, const double); +# endif \} #endif