From 708767cb534b0fdf0beb9dc3d8a035820ced30a5 Mon Sep 17 00:00:00 2001 From: David Wells Date: Fri, 18 Aug 2023 11:59:34 -0400 Subject: [PATCH] Move the internal FEValuesViews functions into a new file. Most of the compilation work for FEValues is actually the inner functions used by FEValuesViews objects. Explicitly instantiating these functions in their own translation unit significantly lowers compilation time. Here are the numbers for release mode: master: - fe_values_views: 7:58 wall time, 4.5 GB max RSS feature: - fe_values_views: 0:32 wall time, 1.6 GB max RSS - fe_values_views_internal: 2:58 wall time, 2.8 GB max RSS I had to adjust the type signatures (on the first ArrayView argument) to avoid instantiating things for 'const float' rather than float et al. --- include/deal.II/fe/fe_values_views.h | 21 +- include/deal.II/fe/fe_values_views_internal.h | 250 ++++ source/fe/CMakeLists.txt | 2 + source/fe/fe_values_views.cc | 1280 +---------------- source/fe/fe_values_views_internal.cc | 1167 +++++++++++++++ source/fe/fe_values_views_internal.inst.in | 204 +++ 6 files changed, 1689 insertions(+), 1235 deletions(-) create mode 100644 include/deal.II/fe/fe_values_views_internal.h create mode 100644 source/fe/fe_values_views_internal.cc create mode 100644 source/fe/fe_values_views_internal.inst.in diff --git a/include/deal.II/fe/fe_values_views.h b/include/deal.II/fe/fe_values_views.h index 5c03c9d153..b8459e04cf 100644 --- a/include/deal.II/fe/fe_values_views.h +++ b/include/deal.II/fe/fe_values_views.h @@ -18,33 +18,18 @@ #include -#include #include -#include -#include -#include -#include +#include #include +#include -#include -#include - -#include #include #include -#include -#include - -#include -#include - -#include #include -#include -#include #include +#include DEAL_II_NAMESPACE_OPEN diff --git a/include/deal.II/fe/fe_values_views_internal.h b/include/deal.II/fe/fe_values_views_internal.h new file mode 100644 index 0000000000..1c53cad801 --- /dev/null +++ b/include/deal.II/fe/fe_values_views_internal.h @@ -0,0 +1,250 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_fe_values_views_internal_h +#define dealii_fe_values_views_internal_h + +#include + +#include +#include +#include +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace FEValuesViews +{ + /** + * Internal namespace for the utility functions used to actually compute + * values by FEValuesViews + */ + namespace internal + { + // Given values of degrees of freedom, evaluate the + // values/gradients/... at quadrature points + + // ------------------------- scalar functions -------------------------- + + /** + * Compute function values for Scalars. + */ + template + void + do_function_values( + const ArrayView &dof_values, + const Table<2, double> & shape_values, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector::type> &values); + + /** + * same code for gradient and Hessian, template argument 'order' to give + * the order of the derivative (= rank of gradient/Hessian tensor) + */ + template + void + do_function_derivatives( + const ArrayView & dof_values, + const Table<2, dealii::Tensor> &shape_derivatives, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &derivatives); + + /** + * Compute Laplacian values for Scalars. + */ + template + void + do_function_laplacians( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_laplacian_type> &laplacians); + + // ----------------------------- vector part --------------------------- + + /** + * Compute function values for Vectors. + */ + template + void + do_function_values( + const ArrayView &dof_values, + const Table<2, double> & shape_values, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &values); + + /** + * same code for gradient and Hessian, template argument 'order' to give + * the order of the derivative (= rank of gradient/Hessian tensor) + */ + template + void + do_function_derivatives( + const ArrayView & dof_values, + const Table<2, dealii::Tensor> &shape_derivatives, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &derivatives); + + /** + * Compute gradient values for Vectors. + */ + template + void + do_function_symmetric_gradients( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &symmetric_gradients); + + /** + * Compute divergence values for Vectors. + */ + template + void + do_function_divergences( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_divergence_type> &divergences); + + /** + * Compute curl values for Vectors. + */ + template + void + do_function_curls( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector::type>::type> &curls); + + /** + * Compute Laplacian values for Vectors. + */ + template + void + do_function_laplacians( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_laplacian_type> &laplacians); + + // ---------------------- symmetric tensor part ------------------------ + + /** + * Compute values for symmetric tensors. + */ + template + void + do_function_values( + const ArrayView & dof_values, + const dealii::Table<2, double> &shape_values, + const std::vector< + typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &values); + + /** + * Compute divergence values for symmetric tensors. + */ + template + void + do_function_divergences( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector< + typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_divergence_type> &divergences); + + // ---------------------- non-symmetric tensor part ------------------------ + + /** + * Compute values for nonsymmetric tensors. + */ + template + void + do_function_values( + const ArrayView & dof_values, + const dealii::Table<2, double> &shape_values, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &values); + + /** + * Compute divergence values for nonsymmetric tensors. + */ + template + void + do_function_divergences( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_divergence_type> &divergences); + + /** + * Compute gradient values for nonsymmetric tensors. + */ + template + void + do_function_gradients( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_gradient_type> &gradients); + } // end of namespace internal +} // namespace FEValuesViews + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index 57b6cc663e..361a293e35 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -70,6 +70,7 @@ set(_separate_src fe_values.cc fe_values_base.cc fe_values_views.cc + fe_values_views_internal.cc mapping_fe_field.cc mapping_fe_field_inst2.cc fe_tools.cc @@ -127,6 +128,7 @@ set(_inst fe_trace.inst.in fe_values_base.inst.in fe_values_views.inst.in + fe_values_views_internal.inst.in fe_values.inst.in fe_wedge_p.inst.in mapping_c1.inst.in diff --git a/source/fe/fe_values_views.cc b/source/fe/fe_values_views.cc index a620da229b..8326181698 100644 --- a/source/fe/fe_values_views.cc +++ b/source/fe/fe_values_views.cc @@ -14,95 +14,51 @@ // --------------------------------------------------------------------- #include -#include -#include #include -#include -#include -#include #include -#include - #include #include #include -#include - -#include -#include +#include #include -#include - -#include -#include -#include - DEAL_II_NAMESPACE_OPEN namespace internal { - template - inline std::vector - make_shape_function_to_row_table(const FiniteElement &fe) - { - std::vector shape_function_to_row_table( - fe.n_dofs_per_cell() * fe.n_components(), numbers::invalid_unsigned_int); - unsigned int row = 0; - for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) - { - // loop over all components that are nonzero for this particular - // shape function. if a component is zero then we leave the - // value in the table unchanged (at the invalid value) - // otherwise it is mapped to the next free entry - unsigned int nth_nonzero_component = 0; - for (unsigned int c = 0; c < fe.n_components(); ++c) - if (fe.get_nonzero_components(i)[c] == true) - { - shape_function_to_row_table[i * fe.n_components() + c] = - row + nth_nonzero_component; - ++nth_nonzero_component; - } - row += fe.n_nonzero_components(i); - } - - return shape_function_to_row_table; - } - namespace { - // Check to see if a DoF value is zero, implying that subsequent operations - // with the value have no effect. - template - struct CheckForZero - { - static bool - value(const Number &value) - { - return value == dealii::internal::NumberType::value(0.0); - } - }; - - // For auto-differentiable numbers, the fact that a DoF value is zero - // does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - // Note that we also want to avoid actually checking the value itself, - // since some AD numbers are not contextually convertible to booleans. - template - struct CheckForZero< - Number, - std::enable_if_t::value>> + template + inline std::vector + make_shape_function_to_row_table(const FiniteElement &fe) { - static bool - value(const Number & /*value*/) - { - return false; - } - }; + std::vector shape_function_to_row_table( + fe.n_dofs_per_cell() * fe.n_components(), + numbers::invalid_unsigned_int); + unsigned int row = 0; + for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) + { + // loop over all components that are nonzero for this particular + // shape function. if a component is zero then we leave the + // value in the table unchanged (at the invalid value) + // otherwise it is mapped to the next free entry + unsigned int nth_nonzero_component = 0; + for (unsigned int c = 0; c < fe.n_components(); ++c) + if (fe.get_nonzero_components(i)[c] == true) + { + shape_function_to_row_table[i * fe.n_components() + c] = + row + nth_nonzero_component; + ++nth_nonzero_component; + } + row += fe.n_nonzero_components(i); + } + + return shape_function_to_row_table; + } } // namespace } // namespace internal @@ -400,1116 +356,6 @@ namespace FEValuesViews - namespace internal - { - // Given values of degrees of freedom, evaluate the - // values/gradients/... at quadrature points - - // ------------------------- scalar functions -------------------------- - template - void - do_function_values( - const ArrayView &dof_values, - const Table<2, double> & shape_values, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector::type> &values) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = values.size(); - - std::fill(values.begin(), - values.end(), - dealii::internal::NumberType::value(0.0)); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component) - { - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is - // zero does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - const double *shape_value_ptr = - &shape_values(shape_function_data[shape_function].row_index, 0); - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - values[q_point] += value * (*shape_value_ptr++); - } - } - - - - // same code for gradient and Hessian, template argument 'order' to give - // the order of the derivative (= rank of gradient/Hessian tensor) - template - void - do_function_derivatives( - const ArrayView & dof_values, - const Table<2, dealii::Tensor> &shape_derivatives, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector< - typename ProductType>::type> - &derivatives) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = derivatives.size(); - - std::fill( - derivatives.begin(), - derivatives.end(), - typename ProductType>::type()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component) - { - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is - // zero does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - const dealii::Tensor *shape_derivative_ptr = - &shape_derivatives[shape_function_data[shape_function].row_index] - [0]; - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - derivatives[q_point] += value * (*shape_derivative_ptr++); - } - } - - - - template - void - do_function_laplacians( - const ArrayView & dof_values, - const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector:: - template solution_laplacian_type> &laplacians) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = laplacians.size(); - - std::fill( - laplacians.begin(), - laplacians.end(), - typename Scalar::template solution_laplacian_type()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component) - { - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is - // zero does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - const dealii::Tensor<2, spacedim> *shape_hessian_ptr = - &shape_hessians[shape_function_data[shape_function].row_index][0]; - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - laplacians[q_point] += value * trace(*shape_hessian_ptr++); - } - } - - - - // ----------------------------- vector part --------------------------- - - template - void - do_function_values( - const ArrayView &dof_values, - const Table<2, double> & shape_values, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector< - typename ProductType>::type> - &values) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = values.size(); - - std::fill( - values.begin(), - values.end(), - typename ProductType>::type()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = - shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is zero - // does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - if (snc != -1) - { - const unsigned int comp = shape_function_data[shape_function] - .single_nonzero_component_index; - const double *shape_value_ptr = &shape_values(snc, 0); - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - values[q_point][comp] += value * (*shape_value_ptr++); - } - else - for (unsigned int d = 0; d < spacedim; ++d) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[d]) - { - const double *shape_value_ptr = &shape_values( - shape_function_data[shape_function].row_index[d], 0); - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - values[q_point][d] += value * (*shape_value_ptr++); - } - } - } - - - - template - void - do_function_derivatives( - const ArrayView & dof_values, - const Table<2, dealii::Tensor> &shape_derivatives, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector< - typename ProductType>::type> - &derivatives) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = derivatives.size(); - - std::fill( - derivatives.begin(), - derivatives.end(), - typename ProductType>::type()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = - shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is zero - // does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - if (snc != -1) - { - const unsigned int comp = shape_function_data[shape_function] - .single_nonzero_component_index; - const dealii::Tensor *shape_derivative_ptr = - &shape_derivatives[snc][0]; - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - derivatives[q_point][comp] += value * (*shape_derivative_ptr++); - } - else - for (unsigned int d = 0; d < spacedim; ++d) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[d]) - { - const dealii::Tensor *shape_derivative_ptr = - &shape_derivatives[shape_function_data[shape_function] - .row_index[d]][0]; - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - derivatives[q_point][d] += - value * (*shape_derivative_ptr++); - } - } - } - - - - template - void - do_function_symmetric_gradients( - const ArrayView & dof_values, - const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector< - typename ProductType>::type> - &symmetric_gradients) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = symmetric_gradients.size(); - - std::fill( - symmetric_gradients.begin(), - symmetric_gradients.end(), - typename ProductType>::type()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = - shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is zero - // does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - if (snc != -1) - { - const unsigned int comp = shape_function_data[shape_function] - .single_nonzero_component_index; - const dealii::Tensor<1, spacedim> *shape_gradient_ptr = - &shape_gradients[snc][0]; - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - symmetric_gradients[q_point] += - value * dealii::SymmetricTensor<2, spacedim>( - symmetrize_single_row(comp, *shape_gradient_ptr++)); - } - else - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - { - typename ProductType>::type - grad; - for (unsigned int d = 0; d < spacedim; ++d) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[d]) - grad[d] = - value * - shape_gradients[shape_function_data[shape_function] - .row_index[d]][q_point]; - symmetric_gradients[q_point] += symmetrize(grad); - } - } - } - - - - template - void - do_function_divergences( - const ArrayView & dof_values, - const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector:: - template solution_divergence_type> &divergences) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = divergences.size(); - - std::fill( - divergences.begin(), - divergences.end(), - typename Vector::template solution_divergence_type()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = - shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is zero - // does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - if (snc != -1) - { - const unsigned int comp = shape_function_data[shape_function] - .single_nonzero_component_index; - const dealii::Tensor<1, spacedim> *shape_gradient_ptr = - &shape_gradients[snc][0]; - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - divergences[q_point] += value * (*shape_gradient_ptr++)[comp]; - } - else - for (unsigned int d = 0; d < spacedim; ++d) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[d]) - { - const dealii::Tensor<1, spacedim> *shape_gradient_ptr = - &shape_gradients[shape_function_data[shape_function] - .row_index[d]][0]; - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - divergences[q_point] += value * (*shape_gradient_ptr++)[d]; - } - } - } - - - - template - void - do_function_curls( - const ArrayView & dof_values, - const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector::type>::type> &curls) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = curls.size(); - - std::fill(curls.begin(), - curls.end(), - typename ProductType< - Number, - typename dealii::internal::CurlType::type>::type()); - - switch (spacedim) - { - case 1: - { - Assert(false, - ExcMessage( - "Computing the curl in 1d is not a useful operation")); - break; - } - - case 2: - { - for (unsigned int shape_function = 0; - shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = shape_function_data[shape_function] - .single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value - // is zero does not imply that its derivatives are zero as - // well. So we can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == - true) - continue; - - if (snc != -1) - { - const dealii::Tensor<1, spacedim> *shape_gradient_ptr = - &shape_gradients[snc][0]; - - Assert(shape_function_data[shape_function] - .single_nonzero_component >= 0, - ExcInternalError()); - // we're in 2d, so the formula for the curl is simple: - if (shape_function_data[shape_function] - .single_nonzero_component_index == 0) - for (unsigned int q_point = 0; - q_point < n_quadrature_points; - ++q_point) - curls[q_point][0] -= - value * (*shape_gradient_ptr++)[1]; - else - for (unsigned int q_point = 0; - q_point < n_quadrature_points; - ++q_point) - curls[q_point][0] += - value * (*shape_gradient_ptr++)[0]; - } - else - // we have multiple non-zero components in the shape - // functions. not all of them must necessarily be within the - // 2-component window this FEValuesViews::Vector object - // considers, however. - { - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[0]) - { - const dealii::Tensor<1, - spacedim> *shape_gradient_ptr = - &shape_gradients[shape_function_data[shape_function] - .row_index[0]][0]; - - for (unsigned int q_point = 0; - q_point < n_quadrature_points; - ++q_point) - curls[q_point][0] -= - value * (*shape_gradient_ptr++)[1]; - } - - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[1]) - { - const dealii::Tensor<1, - spacedim> *shape_gradient_ptr = - &shape_gradients[shape_function_data[shape_function] - .row_index[1]][0]; - - for (unsigned int q_point = 0; - q_point < n_quadrature_points; - ++q_point) - curls[q_point][0] += - value * (*shape_gradient_ptr++)[0]; - } - } - } - break; - } - - case 3: - { - for (unsigned int shape_function = 0; - shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = shape_function_data[shape_function] - .single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value - // is zero does not imply that its derivatives are zero as - // well. So we can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == - true) - continue; - - if (snc != -1) - { - const dealii::Tensor<1, spacedim> *shape_gradient_ptr = - &shape_gradients[snc][0]; - - switch (shape_function_data[shape_function] - .single_nonzero_component_index) - { - case 0: - { - for (unsigned int q_point = 0; - q_point < n_quadrature_points; - ++q_point) - { - curls[q_point][1] += - value * (*shape_gradient_ptr)[2]; - curls[q_point][2] -= - value * (*shape_gradient_ptr++)[1]; - } - - break; - } - - case 1: - { - for (unsigned int q_point = 0; - q_point < n_quadrature_points; - ++q_point) - { - curls[q_point][0] -= - value * (*shape_gradient_ptr)[2]; - curls[q_point][2] += - value * (*shape_gradient_ptr++)[0]; - } - - break; - } - - case 2: - { - for (unsigned int q_point = 0; - q_point < n_quadrature_points; - ++q_point) - { - curls[q_point][0] += - value * (*shape_gradient_ptr)[1]; - curls[q_point][1] -= - value * (*shape_gradient_ptr++)[0]; - } - break; - } - - default: - Assert(false, ExcInternalError()); - } - } - - else - // we have multiple non-zero components in the shape - // functions. not all of them must necessarily be within the - // 3-component window this FEValuesViews::Vector object - // considers, however. - { - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[0]) - { - const dealii::Tensor<1, - spacedim> *shape_gradient_ptr = - &shape_gradients[shape_function_data[shape_function] - .row_index[0]][0]; - - for (unsigned int q_point = 0; - q_point < n_quadrature_points; - ++q_point) - { - curls[q_point][1] += - value * (*shape_gradient_ptr)[2]; - curls[q_point][2] -= - value * (*shape_gradient_ptr++)[1]; - } - } - - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[1]) - { - const dealii::Tensor<1, - spacedim> *shape_gradient_ptr = - &shape_gradients[shape_function_data[shape_function] - .row_index[1]][0]; - - for (unsigned int q_point = 0; - q_point < n_quadrature_points; - ++q_point) - { - curls[q_point][0] -= - value * (*shape_gradient_ptr)[2]; - curls[q_point][2] += - value * (*shape_gradient_ptr++)[0]; - } - } - - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[2]) - { - const dealii::Tensor<1, - spacedim> *shape_gradient_ptr = - &shape_gradients[shape_function_data[shape_function] - .row_index[2]][0]; - - for (unsigned int q_point = 0; - q_point < n_quadrature_points; - ++q_point) - { - curls[q_point][0] += - value * (*shape_gradient_ptr)[1]; - curls[q_point][1] -= - value * (*shape_gradient_ptr++)[0]; - } - } - } - } - } - } - } - - - - template - void - do_function_laplacians( - const ArrayView & dof_values, - const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector:: - template solution_laplacian_type> &laplacians) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = laplacians.size(); - - std::fill( - laplacians.begin(), - laplacians.end(), - typename Vector::template solution_laplacian_type()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = - shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is zero - // does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - if (snc != -1) - { - const unsigned int comp = shape_function_data[shape_function] - .single_nonzero_component_index; - const dealii::Tensor<2, spacedim> *shape_hessian_ptr = - &shape_hessians[snc][0]; - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - laplacians[q_point][comp] += - value * trace(*shape_hessian_ptr++); - } - else - for (unsigned int d = 0; d < spacedim; ++d) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[d]) - { - const dealii::Tensor<2, spacedim> *shape_hessian_ptr = - &shape_hessians[shape_function_data[shape_function] - .row_index[d]][0]; - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - laplacians[q_point][d] += - value * trace(*shape_hessian_ptr++); - } - } - } - - - - // ---------------------- symmetric tensor part ------------------------ - - template - void - do_function_values( - const ArrayView & dof_values, - const dealii::Table<2, double> &shape_values, - const std::vector< - typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData> - &shape_function_data, - std::vector< - typename ProductType>::type> - &values) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = values.size(); - - std::fill( - values.begin(), - values.end(), - typename ProductType>::type()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = - shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is zero - // does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - if (snc != -1) - { - const TableIndices<2> comp = dealii:: - SymmetricTensor<2, spacedim>::unrolled_to_component_indices( - shape_function_data[shape_function] - .single_nonzero_component_index); - const double *shape_value_ptr = &shape_values(snc, 0); - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - values[q_point][comp] += value * (*shape_value_ptr++); - } - else - for (unsigned int d = 0; - d < - dealii::SymmetricTensor<2, spacedim>::n_independent_components; - ++d) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[d]) - { - const TableIndices<2> comp = - dealii::SymmetricTensor<2, spacedim>:: - unrolled_to_component_indices(d); - const double *shape_value_ptr = &shape_values( - shape_function_data[shape_function].row_index[d], 0); - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - values[q_point][comp] += value * (*shape_value_ptr++); - } - } - } - - - - template - void - do_function_divergences( - const ArrayView & dof_values, - const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, - const std::vector< - typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData> - &shape_function_data, - std::vector:: - template solution_divergence_type> &divergences) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = divergences.size(); - - std::fill(divergences.begin(), - divergences.end(), - typename SymmetricTensor<2, dim, spacedim>:: - template solution_divergence_type()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = - shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is zero - // does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - if (snc != -1) - { - const unsigned int comp = shape_function_data[shape_function] - .single_nonzero_component_index; - - const dealii::Tensor<1, spacedim> *shape_gradient_ptr = - &shape_gradients[snc][0]; - - const unsigned int ii = dealii::SymmetricTensor<2, spacedim>:: - unrolled_to_component_indices(comp)[0]; - const unsigned int jj = dealii::SymmetricTensor<2, spacedim>:: - unrolled_to_component_indices(comp)[1]; - - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point, ++shape_gradient_ptr) - { - divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj]; - - if (ii != jj) - divergences[q_point][jj] += - value * (*shape_gradient_ptr)[ii]; - } - } - else - { - for (unsigned int d = 0; - d < - dealii::SymmetricTensor<2, - spacedim>::n_independent_components; - ++d) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[d]) - { - Assert(false, ExcNotImplemented()); - - // the following implementation needs to be looked over -- I - // think it can't be right, because we are in a case where - // there is no single nonzero component - // - // the following is not implemented! we need to consider the - // interplay between multiple non-zero entries in shape - // function and the representation as a symmetric - // second-order tensor - const unsigned int comp = - shape_function_data[shape_function] - .single_nonzero_component_index; - - const dealii::Tensor<1, spacedim> *shape_gradient_ptr = - &shape_gradients[shape_function_data[shape_function] - .row_index[d]][0]; - for (unsigned int q_point = 0; - q_point < n_quadrature_points; - ++q_point, ++shape_gradient_ptr) - { - for (unsigned int j = 0; j < spacedim; ++j) - { - const unsigned int vector_component = - dealii::SymmetricTensor<2, spacedim>:: - component_to_unrolled_index( - TableIndices<2>(comp, j)); - divergences[q_point][vector_component] += - value * (*shape_gradient_ptr++)[j]; - } - } - } - } - } - } - - // ---------------------- non-symmetric tensor part ------------------------ - - template - void - do_function_values( - const ArrayView & dof_values, - const dealii::Table<2, double> &shape_values, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector< - typename ProductType>::type> - &values) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = values.size(); - - std::fill( - values.begin(), - values.end(), - typename ProductType>::type()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = - shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is zero - // does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - if (snc != -1) - { - const unsigned int comp = shape_function_data[shape_function] - .single_nonzero_component_index; - - const TableIndices<2> indices = - dealii::Tensor<2, spacedim>::unrolled_to_component_indices( - comp); - - const double *shape_value_ptr = &shape_values(snc, 0); - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - values[q_point][indices] += value * (*shape_value_ptr++); - } - else - for (unsigned int d = 0; d < dim * dim; ++d) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[d]) - { - const TableIndices<2> indices = - dealii::Tensor<2, spacedim>::unrolled_to_component_indices( - d); - - const double *shape_value_ptr = &shape_values( - shape_function_data[shape_function].row_index[d], 0); - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point) - values[q_point][indices] += value * (*shape_value_ptr++); - } - } - } - - - - template - void - do_function_divergences( - const ArrayView & dof_values, - const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector:: - template solution_divergence_type> &divergences) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = divergences.size(); - - std::fill( - divergences.begin(), - divergences.end(), - typename Tensor<2, dim, spacedim>::template solution_divergence_type< - Number>()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = - shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is zero - // does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - if (snc != -1) - { - const unsigned int comp = shape_function_data[shape_function] - .single_nonzero_component_index; - - const dealii::Tensor<1, spacedim> *shape_gradient_ptr = - &shape_gradients[snc][0]; - - const TableIndices<2> indices = - dealii::Tensor<2, spacedim>::unrolled_to_component_indices( - comp); - const unsigned int ii = indices[0]; - const unsigned int jj = indices[1]; - - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point, ++shape_gradient_ptr) - { - divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj]; - } - } - else - { - for (unsigned int d = 0; d < dim * dim; ++d) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[d]) - { - Assert(false, ExcNotImplemented()); - } - } - } - } - - - - template - void - do_function_gradients( - const ArrayView & dof_values, - const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, - const std::vector::ShapeFunctionData> - &shape_function_data, - std::vector:: - template solution_gradient_type> &gradients) - { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = gradients.size(); - - std::fill( - gradients.begin(), - gradients.end(), - typename Tensor<2, dim, spacedim>::template solution_gradient_type< - Number>()); - - for (unsigned int shape_function = 0; shape_function < dofs_per_cell; - ++shape_function) - { - const int snc = - shape_function_data[shape_function].single_nonzero_component; - - if (snc == -2) - // shape function is zero for the selected components - continue; - - const Number &value = dof_values[shape_function]; - // For auto-differentiable numbers, the fact that a DoF value is zero - // does not imply that its derivatives are zero as well. So we - // can't filter by value for these number types. - if (dealii::internal::CheckForZero::value(value) == true) - continue; - - if (snc != -1) - { - const unsigned int comp = shape_function_data[shape_function] - .single_nonzero_component_index; - - const dealii::Tensor<1, spacedim> *shape_gradient_ptr = - &shape_gradients[snc][0]; - - const TableIndices<2> indices = - dealii::Tensor<2, spacedim>::unrolled_to_component_indices( - comp); - const unsigned int ii = indices[0]; - const unsigned int jj = indices[1]; - - for (unsigned int q_point = 0; q_point < n_quadrature_points; - ++q_point, ++shape_gradient_ptr) - { - gradients[q_point][ii][jj] += value * (*shape_gradient_ptr); - } - } - else - { - for (unsigned int d = 0; d < dim * dim; ++d) - if (shape_function_data[shape_function] - .is_nonzero_shape_function_component[d]) - { - Assert(false, ExcNotImplemented()); - } - } - } - } - - } // end of namespace internal - - - template template void @@ -1529,7 +375,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_values( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_values, shape_function_data, values); @@ -1553,7 +399,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_values( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_values, shape_function_data, values); @@ -1581,7 +427,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_derivatives<1, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, gradients); @@ -1605,7 +451,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_derivatives<1, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, gradients); @@ -1633,7 +479,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_derivatives<2, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_hessians, shape_function_data, hessians); @@ -1657,7 +503,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_derivatives<2, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_hessians, shape_function_data, hessians); @@ -1685,7 +531,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_laplacians( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_hessians, shape_function_data, laplacians); @@ -1709,7 +555,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_laplacians( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_hessians, shape_function_data, laplacians); @@ -1738,7 +584,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_derivatives<3, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_3rd_derivatives, shape_function_data, third_derivatives); @@ -1763,7 +609,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_derivatives<3, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_3rd_derivatives, shape_function_data, third_derivatives); @@ -1789,7 +635,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_values( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_values, shape_function_data, values); @@ -1813,7 +659,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_values( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_values, shape_function_data, values); @@ -1841,7 +687,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_derivatives<1, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, gradients); @@ -1865,7 +711,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_derivatives<1, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, gradients); @@ -1894,7 +740,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_symmetric_gradients( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, symmetric_gradients); @@ -1919,7 +765,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_symmetric_gradients( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, symmetric_gradients); @@ -1948,7 +794,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_divergences( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, divergences); @@ -1972,7 +818,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_divergences( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, divergences); @@ -2000,7 +846,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_curls( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, curls); @@ -2024,7 +870,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_curls( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, curls); @@ -2052,7 +898,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_derivatives<2, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_hessians, shape_function_data, hessians); @@ -2076,7 +922,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_derivatives<2, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_hessians, shape_function_data, hessians); @@ -2109,7 +955,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_laplacians( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_hessians, shape_function_data, laplacians); @@ -2136,7 +982,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_laplacians( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_hessians, shape_function_data, laplacians); @@ -2165,7 +1011,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_derivatives<3, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_3rd_derivatives, shape_function_data, third_derivatives); @@ -2190,7 +1036,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_derivatives<3, dim, spacedim>( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_3rd_derivatives, shape_function_data, third_derivatives); @@ -2216,7 +1062,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_values( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_values, shape_function_data, values); @@ -2240,7 +1086,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_values( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_values, shape_function_data, values); @@ -2269,7 +1115,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_divergences( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, divergences); @@ -2294,7 +1140,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_divergences( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, divergences); @@ -2320,7 +1166,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_values( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_values, shape_function_data, values); @@ -2344,7 +1190,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_values( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_values, shape_function_data, values); @@ -2373,7 +1219,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_divergences( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, divergences); @@ -2397,7 +1243,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_divergences( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, divergences); @@ -2426,7 +1272,7 @@ namespace FEValuesViews fe_values->present_cell.get_interpolated_dof_values(fe_function, dof_values); internal::do_function_gradients( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, gradients); @@ -2450,7 +1296,7 @@ namespace FEValuesViews AssertDimension(dof_values.size(), fe_values->dofs_per_cell); internal::do_function_gradients( - make_array_view(dof_values.begin(), dof_values.end()), + make_const_array_view(dof_values), fe_values->finite_element_output.shape_gradients, shape_function_data, gradients); diff --git a/source/fe/fe_values_views_internal.cc b/source/fe/fe_values_views_internal.cc new file mode 100644 index 0000000000..4a038e2122 --- /dev/null +++ b/source/fe/fe_values_views_internal.cc @@ -0,0 +1,1167 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include +#include + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace FEValuesViews +{ + namespace internal + { + namespace + { + // Check to see if a DoF value is zero, implying that subsequent + // operations with the value have no effect. + template + struct CheckForZero + { + static bool + value(const Number &value) + { + return value == dealii::internal::NumberType::value(0.0); + } + }; + + // For auto-differentiable numbers, the fact that a DoF value is zero + // does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + // Note that we also want to avoid actually checking the value itself, + // since some AD numbers are not contextually convertible to booleans. + template + struct CheckForZero< + Number, + std::enable_if_t::value>> + { + static bool + value(const Number & /*value*/) + { + return false; + } + }; + } // namespace + + template + void + do_function_values( + const ArrayView &dof_values, + const Table<2, double> & shape_values, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector::type> &values) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = values.size(); + + std::fill(values.begin(), + values.end(), + dealii::internal::NumberType::value(0.0)); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component) + { + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is + // zero does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + const double *shape_value_ptr = + &shape_values(shape_function_data[shape_function].row_index, 0); + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + values[q_point] += value * (*shape_value_ptr++); + } + } + + + + template + void + do_function_derivatives( + const ArrayView & dof_values, + const Table<2, dealii::Tensor> &shape_derivatives, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &derivatives) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = derivatives.size(); + + std::fill( + derivatives.begin(), + derivatives.end(), + typename ProductType>::type()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component) + { + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is + // zero does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + const dealii::Tensor *shape_derivative_ptr = + &shape_derivatives[shape_function_data[shape_function].row_index] + [0]; + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + derivatives[q_point] += value * (*shape_derivative_ptr++); + } + } + + + + template + void + do_function_laplacians( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_laplacian_type> &laplacians) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = laplacians.size(); + + std::fill( + laplacians.begin(), + laplacians.end(), + typename Scalar::template solution_laplacian_type()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component) + { + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is + // zero does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + const dealii::Tensor<2, spacedim> *shape_hessian_ptr = + &shape_hessians[shape_function_data[shape_function].row_index][0]; + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + laplacians[q_point] += value * trace(*shape_hessian_ptr++); + } + } + + + + // ----------------------------- vector part --------------------------- + + template + void + do_function_values( + const ArrayView &dof_values, + const Table<2, double> & shape_values, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &values) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = values.size(); + + std::fill( + values.begin(), + values.end(), + typename ProductType>::type()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = + shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is zero + // does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const unsigned int comp = shape_function_data[shape_function] + .single_nonzero_component_index; + const double *shape_value_ptr = &shape_values(snc, 0); + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + values[q_point][comp] += value * (*shape_value_ptr++); + } + else + for (unsigned int d = 0; d < spacedim; ++d) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[d]) + { + const double *shape_value_ptr = &shape_values( + shape_function_data[shape_function].row_index[d], 0); + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + values[q_point][d] += value * (*shape_value_ptr++); + } + } + } + + + + template + void + do_function_derivatives( + const ArrayView & dof_values, + const Table<2, dealii::Tensor> &shape_derivatives, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &derivatives) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = derivatives.size(); + + std::fill( + derivatives.begin(), + derivatives.end(), + typename ProductType>::type()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = + shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is zero + // does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const unsigned int comp = shape_function_data[shape_function] + .single_nonzero_component_index; + const dealii::Tensor *shape_derivative_ptr = + &shape_derivatives[snc][0]; + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + derivatives[q_point][comp] += value * (*shape_derivative_ptr++); + } + else + for (unsigned int d = 0; d < spacedim; ++d) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[d]) + { + const dealii::Tensor *shape_derivative_ptr = + &shape_derivatives[shape_function_data[shape_function] + .row_index[d]][0]; + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + derivatives[q_point][d] += + value * (*shape_derivative_ptr++); + } + } + } + + + + template + void + do_function_symmetric_gradients( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &symmetric_gradients) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = symmetric_gradients.size(); + + std::fill( + symmetric_gradients.begin(), + symmetric_gradients.end(), + typename ProductType>::type()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = + shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is zero + // does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const unsigned int comp = shape_function_data[shape_function] + .single_nonzero_component_index; + const dealii::Tensor<1, spacedim> *shape_gradient_ptr = + &shape_gradients[snc][0]; + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + symmetric_gradients[q_point] += + value * dealii::SymmetricTensor<2, spacedim>( + symmetrize_single_row(comp, *shape_gradient_ptr++)); + } + else + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + { + typename ProductType>::type + grad; + for (unsigned int d = 0; d < spacedim; ++d) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[d]) + grad[d] = + value * + shape_gradients[shape_function_data[shape_function] + .row_index[d]][q_point]; + symmetric_gradients[q_point] += symmetrize(grad); + } + } + } + + + + template + void + do_function_divergences( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_divergence_type> &divergences) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = divergences.size(); + + std::fill( + divergences.begin(), + divergences.end(), + typename Vector::template solution_divergence_type()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = + shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is zero + // does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const unsigned int comp = shape_function_data[shape_function] + .single_nonzero_component_index; + const dealii::Tensor<1, spacedim> *shape_gradient_ptr = + &shape_gradients[snc][0]; + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + divergences[q_point] += value * (*shape_gradient_ptr++)[comp]; + } + else + for (unsigned int d = 0; d < spacedim; ++d) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[d]) + { + const dealii::Tensor<1, spacedim> *shape_gradient_ptr = + &shape_gradients[shape_function_data[shape_function] + .row_index[d]][0]; + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + divergences[q_point] += value * (*shape_gradient_ptr++)[d]; + } + } + } + + + + template + void + do_function_curls( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector::type>::type> &curls) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = curls.size(); + + std::fill(curls.begin(), + curls.end(), + typename ProductType< + Number, + typename dealii::internal::CurlType::type>::type()); + + switch (spacedim) + { + case 1: + { + Assert(false, + ExcMessage( + "Computing the curl in 1d is not a useful operation")); + break; + } + + case 2: + { + for (unsigned int shape_function = 0; + shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = shape_function_data[shape_function] + .single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value + // is zero does not imply that its derivatives are zero as + // well. So we can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const dealii::Tensor<1, spacedim> *shape_gradient_ptr = + &shape_gradients[snc][0]; + + Assert(shape_function_data[shape_function] + .single_nonzero_component >= 0, + ExcInternalError()); + // we're in 2d, so the formula for the curl is simple: + if (shape_function_data[shape_function] + .single_nonzero_component_index == 0) + for (unsigned int q_point = 0; + q_point < n_quadrature_points; + ++q_point) + curls[q_point][0] -= + value * (*shape_gradient_ptr++)[1]; + else + for (unsigned int q_point = 0; + q_point < n_quadrature_points; + ++q_point) + curls[q_point][0] += + value * (*shape_gradient_ptr++)[0]; + } + else + // we have multiple non-zero components in the shape + // functions. not all of them must necessarily be within the + // 2-component window this FEValuesViews::Vector object + // considers, however. + { + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[0]) + { + const dealii::Tensor<1, + spacedim> *shape_gradient_ptr = + &shape_gradients[shape_function_data[shape_function] + .row_index[0]][0]; + + for (unsigned int q_point = 0; + q_point < n_quadrature_points; + ++q_point) + curls[q_point][0] -= + value * (*shape_gradient_ptr++)[1]; + } + + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[1]) + { + const dealii::Tensor<1, + spacedim> *shape_gradient_ptr = + &shape_gradients[shape_function_data[shape_function] + .row_index[1]][0]; + + for (unsigned int q_point = 0; + q_point < n_quadrature_points; + ++q_point) + curls[q_point][0] += + value * (*shape_gradient_ptr++)[0]; + } + } + } + break; + } + + case 3: + { + for (unsigned int shape_function = 0; + shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = shape_function_data[shape_function] + .single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value + // is zero does not imply that its derivatives are zero as + // well. So we can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const dealii::Tensor<1, spacedim> *shape_gradient_ptr = + &shape_gradients[snc][0]; + + switch (shape_function_data[shape_function] + .single_nonzero_component_index) + { + case 0: + { + for (unsigned int q_point = 0; + q_point < n_quadrature_points; + ++q_point) + { + curls[q_point][1] += + value * (*shape_gradient_ptr)[2]; + curls[q_point][2] -= + value * (*shape_gradient_ptr++)[1]; + } + + break; + } + + case 1: + { + for (unsigned int q_point = 0; + q_point < n_quadrature_points; + ++q_point) + { + curls[q_point][0] -= + value * (*shape_gradient_ptr)[2]; + curls[q_point][2] += + value * (*shape_gradient_ptr++)[0]; + } + + break; + } + + case 2: + { + for (unsigned int q_point = 0; + q_point < n_quadrature_points; + ++q_point) + { + curls[q_point][0] += + value * (*shape_gradient_ptr)[1]; + curls[q_point][1] -= + value * (*shape_gradient_ptr++)[0]; + } + break; + } + + default: + Assert(false, ExcInternalError()); + } + } + + else + // we have multiple non-zero components in the shape + // functions. not all of them must necessarily be within the + // 3-component window this FEValuesViews::Vector object + // considers, however. + { + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[0]) + { + const dealii::Tensor<1, + spacedim> *shape_gradient_ptr = + &shape_gradients[shape_function_data[shape_function] + .row_index[0]][0]; + + for (unsigned int q_point = 0; + q_point < n_quadrature_points; + ++q_point) + { + curls[q_point][1] += + value * (*shape_gradient_ptr)[2]; + curls[q_point][2] -= + value * (*shape_gradient_ptr++)[1]; + } + } + + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[1]) + { + const dealii::Tensor<1, + spacedim> *shape_gradient_ptr = + &shape_gradients[shape_function_data[shape_function] + .row_index[1]][0]; + + for (unsigned int q_point = 0; + q_point < n_quadrature_points; + ++q_point) + { + curls[q_point][0] -= + value * (*shape_gradient_ptr)[2]; + curls[q_point][2] += + value * (*shape_gradient_ptr++)[0]; + } + } + + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[2]) + { + const dealii::Tensor<1, + spacedim> *shape_gradient_ptr = + &shape_gradients[shape_function_data[shape_function] + .row_index[2]][0]; + + for (unsigned int q_point = 0; + q_point < n_quadrature_points; + ++q_point) + { + curls[q_point][0] += + value * (*shape_gradient_ptr)[1]; + curls[q_point][1] -= + value * (*shape_gradient_ptr++)[0]; + } + } + } + } + } + } + } + + + + template + void + do_function_laplacians( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<2, spacedim>> &shape_hessians, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_laplacian_type> &laplacians) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = laplacians.size(); + + std::fill( + laplacians.begin(), + laplacians.end(), + typename Vector::template solution_laplacian_type()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = + shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is zero + // does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const unsigned int comp = shape_function_data[shape_function] + .single_nonzero_component_index; + const dealii::Tensor<2, spacedim> *shape_hessian_ptr = + &shape_hessians[snc][0]; + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + laplacians[q_point][comp] += + value * trace(*shape_hessian_ptr++); + } + else + for (unsigned int d = 0; d < spacedim; ++d) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[d]) + { + const dealii::Tensor<2, spacedim> *shape_hessian_ptr = + &shape_hessians[shape_function_data[shape_function] + .row_index[d]][0]; + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + laplacians[q_point][d] += + value * trace(*shape_hessian_ptr++); + } + } + } + + + + // ---------------------- symmetric tensor part ------------------------ + + template + void + do_function_values( + const ArrayView & dof_values, + const dealii::Table<2, double> &shape_values, + const std::vector< + typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &values) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = values.size(); + + std::fill( + values.begin(), + values.end(), + typename ProductType>::type()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = + shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is zero + // does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const TableIndices<2> comp = dealii:: + SymmetricTensor<2, spacedim>::unrolled_to_component_indices( + shape_function_data[shape_function] + .single_nonzero_component_index); + const double *shape_value_ptr = &shape_values(snc, 0); + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + values[q_point][comp] += value * (*shape_value_ptr++); + } + else + for (unsigned int d = 0; + d < + dealii::SymmetricTensor<2, spacedim>::n_independent_components; + ++d) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[d]) + { + const TableIndices<2> comp = + dealii::SymmetricTensor<2, spacedim>:: + unrolled_to_component_indices(d); + const double *shape_value_ptr = &shape_values( + shape_function_data[shape_function].row_index[d], 0); + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + values[q_point][comp] += value * (*shape_value_ptr++); + } + } + } + + + + template + void + do_function_divergences( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector< + typename SymmetricTensor<2, dim, spacedim>::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_divergence_type> &divergences) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = divergences.size(); + + std::fill(divergences.begin(), + divergences.end(), + typename SymmetricTensor<2, dim, spacedim>:: + template solution_divergence_type()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = + shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is zero + // does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const unsigned int comp = shape_function_data[shape_function] + .single_nonzero_component_index; + + const dealii::Tensor<1, spacedim> *shape_gradient_ptr = + &shape_gradients[snc][0]; + + const unsigned int ii = dealii::SymmetricTensor<2, spacedim>:: + unrolled_to_component_indices(comp)[0]; + const unsigned int jj = dealii::SymmetricTensor<2, spacedim>:: + unrolled_to_component_indices(comp)[1]; + + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point, ++shape_gradient_ptr) + { + divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj]; + + if (ii != jj) + divergences[q_point][jj] += + value * (*shape_gradient_ptr)[ii]; + } + } + else + { + for (unsigned int d = 0; + d < + dealii::SymmetricTensor<2, + spacedim>::n_independent_components; + ++d) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[d]) + { + Assert(false, ExcNotImplemented()); + + // the following implementation needs to be looked over -- I + // think it can't be right, because we are in a case where + // there is no single nonzero component + // + // the following is not implemented! we need to consider the + // interplay between multiple non-zero entries in shape + // function and the representation as a symmetric + // second-order tensor + const unsigned int comp = + shape_function_data[shape_function] + .single_nonzero_component_index; + + const dealii::Tensor<1, spacedim> *shape_gradient_ptr = + &shape_gradients[shape_function_data[shape_function] + .row_index[d]][0]; + for (unsigned int q_point = 0; + q_point < n_quadrature_points; + ++q_point, ++shape_gradient_ptr) + { + for (unsigned int j = 0; j < spacedim; ++j) + { + const unsigned int vector_component = + dealii::SymmetricTensor<2, spacedim>:: + component_to_unrolled_index( + TableIndices<2>(comp, j)); + divergences[q_point][vector_component] += + value * (*shape_gradient_ptr++)[j]; + } + } + } + } + } + } + + // ---------------------- non-symmetric tensor part ------------------------ + + template + void + do_function_values( + const ArrayView & dof_values, + const dealii::Table<2, double> &shape_values, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector< + typename ProductType>::type> + &values) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = values.size(); + + std::fill( + values.begin(), + values.end(), + typename ProductType>::type()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = + shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is zero + // does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const unsigned int comp = shape_function_data[shape_function] + .single_nonzero_component_index; + + const TableIndices<2> indices = + dealii::Tensor<2, spacedim>::unrolled_to_component_indices( + comp); + + const double *shape_value_ptr = &shape_values(snc, 0); + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + values[q_point][indices] += value * (*shape_value_ptr++); + } + else + for (unsigned int d = 0; d < dim * dim; ++d) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[d]) + { + const TableIndices<2> indices = + dealii::Tensor<2, spacedim>::unrolled_to_component_indices( + d); + + const double *shape_value_ptr = &shape_values( + shape_function_data[shape_function].row_index[d], 0); + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point) + values[q_point][indices] += value * (*shape_value_ptr++); + } + } + } + + + + template + void + do_function_divergences( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_divergence_type> &divergences) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = divergences.size(); + + std::fill( + divergences.begin(), + divergences.end(), + typename Tensor<2, dim, spacedim>::template solution_divergence_type< + Number>()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = + shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is zero + // does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const unsigned int comp = shape_function_data[shape_function] + .single_nonzero_component_index; + + const dealii::Tensor<1, spacedim> *shape_gradient_ptr = + &shape_gradients[snc][0]; + + const TableIndices<2> indices = + dealii::Tensor<2, spacedim>::unrolled_to_component_indices( + comp); + const unsigned int ii = indices[0]; + const unsigned int jj = indices[1]; + + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point, ++shape_gradient_ptr) + { + divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj]; + } + } + else + { + for (unsigned int d = 0; d < dim * dim; ++d) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[d]) + { + Assert(false, ExcNotImplemented()); + } + } + } + } + + + + template + void + do_function_gradients( + const ArrayView & dof_values, + const Table<2, dealii::Tensor<1, spacedim>> &shape_gradients, + const std::vector::ShapeFunctionData> + &shape_function_data, + std::vector:: + template solution_gradient_type> &gradients) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = gradients.size(); + + std::fill( + gradients.begin(), + gradients.end(), + typename Tensor<2, dim, spacedim>::template solution_gradient_type< + Number>()); + + for (unsigned int shape_function = 0; shape_function < dofs_per_cell; + ++shape_function) + { + const int snc = + shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const Number &value = dof_values[shape_function]; + // For auto-differentiable numbers, the fact that a DoF value is zero + // does not imply that its derivatives are zero as well. So we + // can't filter by value for these number types. + if (CheckForZero::value(value) == true) + continue; + + if (snc != -1) + { + const unsigned int comp = shape_function_data[shape_function] + .single_nonzero_component_index; + + const dealii::Tensor<1, spacedim> *shape_gradient_ptr = + &shape_gradients[snc][0]; + + const TableIndices<2> indices = + dealii::Tensor<2, spacedim>::unrolled_to_component_indices( + comp); + const unsigned int ii = indices[0]; + const unsigned int jj = indices[1]; + + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point, ++shape_gradient_ptr) + { + gradients[q_point][ii][jj] += value * (*shape_gradient_ptr); + } + } + else + { + for (unsigned int d = 0; d < dim * dim; ++d) + if (shape_function_data[shape_function] + .is_nonzero_shape_function_component[d]) + { + Assert(false, ExcNotImplemented()); + } + } + } + } + } // end of namespace internal +} // namespace FEValuesViews + + + +/*------------------------------- Explicit Instantiations -------------*/ + +#include "fe_values_views_internal.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe_values_views_internal.inst.in b/source/fe/fe_values_views_internal.inst.in new file mode 100644 index 0000000000..2f4cb3af9e --- /dev/null +++ b/source/fe/fe_values_views_internal.inst.in @@ -0,0 +1,204 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef DOXYGEN + +for (S : ALL_SCALAR_TYPES; deal_II_dimension : DIMENSIONS; + deal_II_space_dimension : SPACE_DIMENSIONS) + { +# if deal_II_dimension <= deal_II_space_dimension + namespace FEValuesViews + + { + namespace internal + + { + template void + do_function_values( + const ArrayView &, + const Table<2, double> &, + const std::vector::ShapeFunctionData> + &, + std::vector::type> &); + + + template void + do_function_laplacians( + const ArrayView &, + const Table<2, dealii::Tensor<2, deal_II_space_dimension>> &, + const std::vector::ShapeFunctionData> + &, + std::vector:: + template solution_laplacian_type> &); + + template void + do_function_values( + const ArrayView &, + const Table<2, double> &, + const std::vector::ShapeFunctionData> + &, + std::vector< + ProductType>::type> + &); + + template void + do_function_symmetric_gradients( + const ArrayView &, + const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &, + const std::vector::ShapeFunctionData> + &, + std::vector>::type> &); + + template void + do_function_divergences( + const ArrayView &, + const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &, + const std::vector::ShapeFunctionData> + &, + std::vector:: + template solution_divergence_type> &); + + template void + do_function_curls( + const ArrayView &, + const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &, + const std::vector::ShapeFunctionData> + &, + std::vector::type>::type> + &); + + template void + do_function_laplacians( + const ArrayView &, + const Table<2, dealii::Tensor<2, deal_II_space_dimension>> &, + const std::vector::ShapeFunctionData> + &, + std::vector:: + template solution_laplacian_type> &); + + template void + do_function_values( + const ArrayView &, + const dealii::Table<2, double> &, + const std::vector< + SymmetricTensor<2, deal_II_dimension, deal_II_space_dimension>:: + ShapeFunctionData> &, + std::vector>::type> &); + + template void + do_function_divergences( + const ArrayView &, + const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &, + const std::vector< + SymmetricTensor<2, deal_II_dimension, deal_II_space_dimension>:: + ShapeFunctionData> &, + std::vector< + SymmetricTensor<2, deal_II_dimension, deal_II_space_dimension>:: + template solution_divergence_type> &); + + template void + do_function_values( + const ArrayView &, + const dealii::Table<2, double> &, + const std::vector< + Tensor<2, deal_II_dimension, deal_II_space_dimension>:: + ShapeFunctionData> &, + std::vector< + ProductType>::type> + &); + + template void + do_function_divergences( + const ArrayView &, + const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &, + const std::vector< + Tensor<2, deal_II_dimension, deal_II_space_dimension>:: + ShapeFunctionData> &, + std::vector:: + template solution_divergence_type> &); + + template void + do_function_gradients( + const ArrayView &, + const Table<2, dealii::Tensor<1, deal_II_space_dimension>> &, + const std::vector< + Tensor<2, deal_II_dimension, deal_II_space_dimension>:: + ShapeFunctionData> &, + std::vector:: + template solution_gradient_type> &); + + \} + \} +# endif + } + +for (S : ALL_SCALAR_TYPES; deal_II_dimension : DIMENSIONS; ORDER : DIMENSIONS; + deal_II_space_dimension : SPACE_DIMENSIONS) + { +# if deal_II_dimension <= deal_II_space_dimension + namespace FEValuesViews + \{ + namespace internal + \{ + template void + do_function_derivatives( + const ArrayView &, + const Table<2, dealii::Tensor> &, + const std::vector::ShapeFunctionData> + &, + std::vector< + ProductType>::type> + &); + + template void + do_function_derivatives( + const ArrayView &, + const Table<2, dealii::Tensor> &, + const std::vector::ShapeFunctionData> + &, + std::vector>::type> &); + \} + \} +# endif + } + + +#endif -- 2.39.5