From 519d8b9b6ded8c357068136ba6d1cef3a1312882 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Mon, 22 Feb 2016 13:37:20 +0100 Subject: [PATCH] extend VectorTools to complex-valued PETSc vector --- doc/news/changes.h | 5 +- include/deal.II/numerics/vector_tools.h | 58 ++++++------ .../deal.II/numerics/vector_tools.templates.h | 90 ++++++++++++------- .../numerics/vector_tools_interpolate.inst.in | 8 +- .../numerics/vector_tools_mean_value.inst.in | 4 +- .../numerics/vector_tools_point_value.inst.in | 24 ++--- 6 files changed, 109 insertions(+), 80 deletions(-) diff --git a/doc/news/changes.h b/doc/news/changes.h index 6bff97718f..1f84ee4d2d 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -38,8 +38,9 @@ inconvenience this causes.

    -
  1. Changed: ConstraintMatrix::distribute_local_to_global() now requires - matching data types. This is done to correctly handle complex-valued case. +
  2. Changed: ConstraintMatrix::distribute_local_to_global() and numerous + functions in VectorTools namespace now require matching data types. + This is done to correctly handle complex-valued case.
    (Denis Davydov, 2016/02/22)
    3. diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 635ee356d2..64f046f0c4 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -403,7 +403,7 @@ namespace VectorTools template class DoFHandlerType> void interpolate (const Mapping &mapping, const DoFHandlerType &dof, - const Function &function, + const Function &function, VectorType &vec); /** @@ -412,7 +412,7 @@ namespace VectorTools */ template void interpolate (const DoFHandlerType &dof, - const Function &function, + const Function &function, VectorType &vec); /** @@ -489,7 +489,7 @@ namespace VectorTools interpolate_based_on_material_id (const Mapping &mapping, const DoFHandlerType &dof_handler, - const std::map *> &function_map, + const std::map *> &function_map, VectorType &dst, const ComponentMask &component_mask = ComponentMask()); @@ -592,7 +592,7 @@ namespace VectorTools const DoFHandler &dof, const ConstraintMatrix &constraints, const Quadrature &quadrature, - const Function &function, + const Function &function, VectorType &vec, const bool enforce_zero_boundary = false, const Quadrature &q_boundary = (dim > 1 ? @@ -608,7 +608,7 @@ namespace VectorTools void project (const DoFHandler &dof, const ConstraintMatrix &constraints, const Quadrature &quadrature, - const Function &function, + const Function &function, VectorType &vec, const bool enforce_zero_boundary = false, const Quadrature &q_boundary = (dim > 1 ? @@ -625,7 +625,7 @@ namespace VectorTools const hp::DoFHandler &dof, const ConstraintMatrix &constraints, const hp::QCollection &quadrature, - const Function &function, + const Function &function, VectorType &vec, const bool enforce_zero_boundary = false, const hp::QCollection &q_boundary = hp::QCollection(dim > 1 ? @@ -641,7 +641,7 @@ namespace VectorTools void project (const hp::DoFHandler &dof, const ConstraintMatrix &constraints, const hp::QCollection &quadrature, - const Function &function, + const Function &function, VectorType &vec, const bool enforce_zero_boundary = false, const hp::QCollection &q_boundary = hp::QCollection(dim > 1 ? @@ -1989,8 +1989,8 @@ namespace VectorTools template void point_difference (const DoFHandler &dof, const VectorType &fe_function, - const Function &exact_solution, - Vector &difference, + const Function &exact_solution, + Vector &difference, const Point &point); /** @@ -2009,8 +2009,8 @@ namespace VectorTools void point_difference (const Mapping &mapping, const DoFHandler &dof, const VectorType &fe_function, - const Function &exact_solution, - Vector &difference, + const Function &exact_solution, + Vector &difference, const Point &point); /** @@ -2029,7 +2029,7 @@ namespace VectorTools point_value (const DoFHandler &dof, const VectorType &fe_function, const Point &point, - Vector &value); + Vector &value); /** * Same as above for hp. @@ -2042,7 +2042,7 @@ namespace VectorTools point_value (const hp::DoFHandler &dof, const VectorType &fe_function, const Point &point, - Vector &value); + Vector &value); /** * Evaluate a scalar finite element function defined by the given DoFHandler @@ -2060,7 +2060,7 @@ namespace VectorTools * exception of type VectorTools::ExcPointNotAvailableHere is thrown. */ template - double + typename VectorType::value_type point_value (const DoFHandler &dof, const VectorType &fe_function, const Point &point); @@ -2072,7 +2072,7 @@ namespace VectorTools * exception of type VectorTools::ExcPointNotAvailableHere is thrown. */ template - double + typename VectorType::value_type point_value (const hp::DoFHandler &dof, const VectorType &fe_function, const Point &point); @@ -2094,7 +2094,7 @@ namespace VectorTools const DoFHandler &dof, const VectorType &fe_function, const Point &point, - Vector &value); + Vector &value); /** * Same as above for hp. @@ -2108,7 +2108,7 @@ namespace VectorTools const hp::DoFHandler &dof, const VectorType &fe_function, const Point &point, - Vector &value); + Vector &value); /** * Evaluate a scalar finite element function defined by the given DoFHandler @@ -2122,7 +2122,7 @@ namespace VectorTools * exception of type VectorTools::ExcPointNotAvailableHere is thrown. */ template - double + typename VectorType::value_type point_value (const Mapping &mapping, const DoFHandler &dof, const VectorType &fe_function, @@ -2135,7 +2135,7 @@ namespace VectorTools * exception of type VectorTools::ExcPointNotAvailableHere is thrown. */ template - double + typename VectorType::value_type point_value (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const VectorType &fe_function, @@ -2346,21 +2346,23 @@ namespace VectorTools * mean value and subtract it right inside the evaluation routine. */ template - double compute_mean_value (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &quadrature, - const VectorType &v, - const unsigned int component); + typename VectorType::value_type + compute_mean_value (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &quadrature, + const VectorType &v, + const unsigned int component); /** * Calls the other compute_mean_value() function, see above, with * mapping=MappingQGeneric@(1). */ template - double compute_mean_value (const DoFHandler &dof, - const Quadrature &quadrature, - const VectorType &v, - const unsigned int component); + typename VectorType::value_type + compute_mean_value (const DoFHandler &dof, + const Quadrature &quadrature, + const VectorType &v, + const unsigned int component); //@} /** * Geometrical interpolation diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index e6889a36b4..6ed3e4d05b 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -77,9 +77,10 @@ namespace VectorTools template class DoFHandlerType> void interpolate (const Mapping &mapping, const DoFHandlerType &dof, - const Function &function, + const Function &function, VectorType &vec) { + typedef typename VectorType::value_type number; Assert (vec.size() == dof.n_dofs(), ExcDimensionMismatch (vec.size(), dof.n_dofs())); Assert (dof.get_fe().n_components() == function.n_components, @@ -187,8 +188,8 @@ namespace VectorTools // have two versions, one for system fe // and one for scalar ones, to take the // more efficient one respectively - std::vector > function_values_scalar(fe.size()); - std::vector > > function_values_system(fe.size()); + std::vector > function_values_scalar(fe.size()); + std::vector > > function_values_system(fe.size()); // Make a quadrature rule from support points // to feed it into FEValues @@ -234,7 +235,7 @@ namespace VectorTools // all components function_values_system[fe_index] .resize (n_rep_points[fe_index], - Vector (fe[fe_index].n_components())); + Vector (fe[fe_index].n_components())); function.vector_value_list (rep_points, function_values_system[fe_index]); // distribute the function @@ -273,7 +274,7 @@ namespace VectorTools template void interpolate (const DoFHandlerType &dof, - const Function &function, + const Function &function, VectorType &vec) { interpolate(StaticMappingQ1 cell_data_1(dof_1.get_fe().dofs_per_cell); - Vector cell_data_2(dof_2.get_fe().dofs_per_cell); + typedef typename OutVector::value_type number; + Vector cell_data_1(dof_1.get_fe().dofs_per_cell); + Vector cell_data_2(dof_2.get_fe().dofs_per_cell); std::vector touch_count (dof_2.n_dofs(), 0); //TODO: check on datatype... kinda strange (UK) std::vector local_dof_indices (dof_2.get_fe().dofs_per_cell); @@ -6102,6 +6104,25 @@ namespace VectorTools std::vector >(n_components)); } + namespace + { + 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. Panick 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(); + } + } + // avoid compiling inner function for many vector types when we always // really do the same thing by putting the main work into this helper @@ -6219,6 +6240,7 @@ namespace VectorTools } double diff = 0; + Number diff_mean = 0; // First work on function values: switch (norm) @@ -6227,10 +6249,10 @@ namespace VectorTools // Compute values in quadrature points and integrate for (unsigned int q=0; q(data.psi_values[q](k)*data.psi_values[q](k)), + static_cast(numbers::NumberTraits::abs_square(data.psi_values[q](k))), exponent/2.) * data.weight_vectors[q](k); diff += sum * fe_values.JxW(q); } @@ -6260,7 +6282,7 @@ namespace VectorTools { double sum = 0; for (unsigned int k=0; k::abs_square(data.psi_values[q](k)) * data.weight_vectors[q](k); diff += sum * fe_values.JxW(q); } @@ -6299,7 +6321,7 @@ namespace VectorTools double sum = 0; for (unsigned int k=0; k(data.psi_grads[q][k]*data.psi_grads[q][k]), + data.psi_grads[q][k].norm_square(), exponent/2.) * data.weight_vectors[q](k); diff += sum * fe_values.JxW(q); } @@ -6312,7 +6334,7 @@ namespace VectorTools { double sum = 0; for (unsigned int k=0; k::abs_square(sum) * fe_values.JxW(q); } diff = std::sqrt(diff); break; @@ -6354,6 +6376,9 @@ namespace VectorTools break; } + if (norm == mean) + diff = mean_to_double(diff_mean); + // append result of this cell to the end of the vector AssertIsFinite(diff); return diff; @@ -6554,8 +6579,8 @@ namespace VectorTools void point_difference (const DoFHandler &dof, const VectorType &fe_function, - const Function &exact_function, - Vector &difference, + const Function &exact_function, + Vector &difference, const Point &point) { point_difference(StaticMappingQ1::mapping, @@ -6572,8 +6597,8 @@ namespace VectorTools point_difference (const Mapping &mapping, const DoFHandler &dof, const VectorType &fe_function, - const Function &exact_function, - Vector &difference, + const Function &exact_function, + Vector &difference, const Point &point) { typedef typename VectorType::value_type Number; @@ -6618,7 +6643,7 @@ namespace VectorTools point_value (const DoFHandler &dof, const VectorType &fe_function, const Point &point, - Vector &value) + Vector &value) { point_value (StaticMappingQ1::mapping, @@ -6634,7 +6659,7 @@ namespace VectorTools point_value (const hp::DoFHandler &dof, const VectorType &fe_function, const Point &point, - Vector &value) + Vector &value) { point_value(hp::StaticMappingQ1::mapping_collection, dof, @@ -6645,7 +6670,7 @@ namespace VectorTools template - double + typename VectorType::value_type point_value (const DoFHandler &dof, const VectorType &fe_function, const Point &point) @@ -6658,7 +6683,7 @@ namespace VectorTools template - double + typename VectorType::value_type point_value (const hp::DoFHandler &dof, const VectorType &fe_function, const Point &point) @@ -6676,7 +6701,7 @@ namespace VectorTools const DoFHandler &dof, const VectorType &fe_function, const Point &point, - Vector &value) + Vector &value) { typedef typename VectorType::value_type Number; const FiniteElement &fe = dof.get_fe(); @@ -6717,7 +6742,7 @@ namespace VectorTools const hp::DoFHandler &dof, const VectorType &fe_function, const Point &point, - Vector &value) + Vector &value) { typedef typename VectorType::value_type Number; const hp::FECollection &fe = dof.get_fe(); @@ -6752,7 +6777,7 @@ namespace VectorTools template - double + typename VectorType::value_type point_value (const Mapping &mapping, const DoFHandler &dof, const VectorType &fe_function, @@ -6761,7 +6786,7 @@ namespace VectorTools Assert(dof.get_fe().n_components() == 1, ExcMessage ("Finite element is not scalar as is necessary for this function")); - Vector value(1); + Vector value(1); point_value(mapping, dof, fe_function, point, value); return value(0); @@ -6769,7 +6794,7 @@ namespace VectorTools template - double + typename VectorType::value_type point_value (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const VectorType &fe_function, @@ -6778,7 +6803,7 @@ namespace VectorTools Assert(dof.get_fe().n_components() == 1, ExcMessage ("Finite element is not scalar as is necessary for this function")); - Vector value(1); + Vector value(1); point_value(mapping, dof, fe_function, point, value); return value(0); @@ -6992,7 +7017,8 @@ namespace VectorTools for (unsigned int i=0; i - double + typename VectorType::value_type compute_mean_value (const Mapping &mapping, const DoFHandler &dof, const Quadrature &quadrature, @@ -7097,7 +7123,7 @@ namespace VectorTools template - double + typename VectorType::value_type compute_mean_value (const DoFHandler &dof, const Quadrature &quadrature, const VectorType &v, diff --git a/source/numerics/vector_tools_interpolate.inst.in b/source/numerics/vector_tools_interpolate.inst.in index bd96acf66f..454b944ccf 100644 --- a/source/numerics/vector_tools_interpolate.inst.in +++ b/source/numerics/vector_tools_interpolate.inst.in @@ -23,13 +23,13 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens void interpolate (const Mapping &, const DoFHandler&, - const Function&, + const Function&, VEC&); template void interpolate (const DoFHandler&, - const Function&, + const Function&, VEC&); template @@ -64,12 +64,12 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens void interpolate (const Mapping&, const hp::DoFHandler&, - const Function&, + const Function&, VEC&); template void interpolate (const hp::DoFHandler&, - const Function&, + const Function&, VEC&); template diff --git a/source/numerics/vector_tools_mean_value.inst.in b/source/numerics/vector_tools_mean_value.inst.in index 6892ce087c..06d304eb70 100644 --- a/source/numerics/vector_tools_mean_value.inst.in +++ b/source/numerics/vector_tools_mean_value.inst.in @@ -20,7 +20,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens namespace VectorTools \{ template - double compute_mean_value + VEC::value_type compute_mean_value (const Mapping&, const DoFHandler&, const Quadrature&, @@ -28,7 +28,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const unsigned int); template - double compute_mean_value + VEC::value_type compute_mean_value (const DoFHandler&, const Quadrature&, const VEC&, diff --git a/source/numerics/vector_tools_point_value.inst.in b/source/numerics/vector_tools_point_value.inst.in index c5ee68d2aa..4c77fab747 100644 --- a/source/numerics/vector_tools_point_value.inst.in +++ b/source/numerics/vector_tools_point_value.inst.in @@ -27,10 +27,10 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const hp::DoFHandler&, const VEC&, const Point&, - Vector&); + Vector&); template - double point_value ( + VEC::value_type point_value ( const hp::DoFHandler&, const VEC&, const Point&); @@ -41,10 +41,10 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const hp::DoFHandler&, const VEC&, const Point&, - Vector&); + Vector&); template - double point_value ( + VEC::value_type point_value ( const hp::MappingCollection&, const hp::DoFHandler&, const VEC&, @@ -54,8 +54,8 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens void point_difference ( const DoFHandler&, const VEC&, - const Function&, - Vector&, + const Function&, + Vector&, const Point&); template @@ -63,8 +63,8 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const Mapping&, const DoFHandler&, const VEC&, - const Function&, - Vector&, + const Function&, + Vector&, const Point&); template @@ -72,10 +72,10 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const DoFHandler&, const VEC&, const Point&, - Vector&); + Vector&); template - double point_value ( + VEC::value_type point_value ( const DoFHandler&, const VEC&, const Point&); @@ -86,10 +86,10 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens const DoFHandler&, const VEC&, const Point&, - Vector&); + Vector&); template - double point_value ( + VEC::value_type point_value ( const Mapping&, const DoFHandler&, const VEC&, -- 2.39.5