From 3b684d9266691877c4bdd27c6cbba8e6edbe72fc Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Tue, 23 Feb 2016 14:33:12 +0100 Subject: [PATCH] fixed a bug in FEFieldFunction::vector_gradient() introduced in 0dc4d47 --- .../numerics/fe_field_function.templates.h | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/include/deal.II/numerics/fe_field_function.templates.h b/include/deal.II/numerics/fe_field_function.templates.h index 025d9cf59f..aea385d4b6 100644 --- a/include/deal.II/numerics/fe_field_function.templates.h +++ b/include/deal.II/numerics/fe_field_function.templates.h @@ -114,6 +114,7 @@ namespace Functions vector_gradient (const Point &p, std::vector > &gradients) const { + typedef typename VectorType::value_type number; Assert (gradients.size() == n_components, ExcDimensionMismatch(gradients.size(), n_components)); typename DoFHandlerType::active_cell_iterator cell = cell_hint.get(); @@ -140,7 +141,25 @@ namespace Functions FEValues fe_v(mapping, cell->get_fe(), quad, update_gradients); fe_v.reinit(cell); - fe_v.get_function_gradients(data_vector, gradients); + + if (n_components == 1) + { + // the size of the @p gradients coincidentally coincides + // with the number of quadrature points we evaluate the function at. + fe_v.get_function_gradients(data_vector, gradients); + } + else + { + // Unfortunately we still need a temporary argument as we want to + // evaluate a gradient of a (generally) multicomponent function at + // a single quadrature point. Note that the first std::vector<> is related + // to the number of quadrature points (always one here), whereas the second + // to the number of components. + std::vector< std::vector > > vgrads + (1, std::vector >(n_components) ); + fe_v.get_function_gradients(data_vector, vgrads); + gradients = vgrads[0]; + } } -- 2.39.5