![]() |
deal.II version 9.7.1
|
#include <deal.II/matrix_free/evaluation_kernels.h>
This struct implements the action of the inverse mass matrix operation with user-provided coefficients at quadrature points (in contrast to CellwiseInverseMassMatrixImplBasic, which implicitly uses `1/(|J|xW)' as coefficient). */ template <int dim, typename Number> struct CellwiseInverseMassMatrixImplFlexible { using Number2 = typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
template <int fe_degree, int = 0> static bool run(const unsigned int n_desired_components, const FEEvaluationData<dim, Number, false> &fe_eval, const ArrayView<const Number> &inverse_coefficients, const bool dyadic_coefficients, const Number *in_array, Number *out_array) { const unsigned int given_degree = (fe_degree > -1) ? fe_degree : fe_eval.get_shape_info().data.front().fe_degree;
const unsigned int dofs_per_component = Utilities::pow(given_degree + 1, dim);
Assert(inverse_coefficients.size() > 0 && inverse_coefficients.size() % dofs_per_component == 0, ExcMessage( "Expected diagonal to be a multiple of scalar dof per cells"));
if (!dyadic_coefficients) { if (inverse_coefficients.size() != dofs_per_component) AssertDimension(n_desired_components * dofs_per_component, inverse_coefficients.size()); } else { AssertDimension(n_desired_components * n_desired_components * dofs_per_component, inverse_coefficients.size()); }
Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); Assert(fe_eval.get_shape_info().element_type <= MatrixFreeFunctions::tensor_symmetric_no_collocation, ExcNotImplemented());
EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree + 1, fe_degree + 1, Number, Number2> evaluator({}, {}, fe_eval.get_shape_info().data.front().inverse_shape_values_eo, given_degree + 1, given_degree + 1);
const Number *in = in_array; Number *out = out_array;
const Number *inv_coefficient = inverse_coefficients.data();
const unsigned int shift_coefficient = inverse_coefficients.size() > dofs_per_component ? dofs_per_component : 0;
const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components; const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1;
for (unsigned int d = 0; d < n_comp_outer; ++d) { for (unsigned int di = 0; di < n_comp_inner; ++di) { const Number *in_ = in + di * dofs_per_component; Number *out_ = out + di * dofs_per_component; evaluator.template hessians<0, true, false>(in_, out_); if (dim > 1) evaluator.template hessians<1, true, false>(out_, out_); if (dim > 2) evaluator.template hessians<2, true, false>(out_, out_); } if (dyadic_coefficients) { const auto n_coeff_components = n_desired_components * n_desired_components; if (n_desired_components == dim) { for (unsigned int q = 0; q < dofs_per_component; ++q) vmult<dim>(&inv_coefficient[q * n_coeff_components], &in[q], &out[q], dofs_per_component); } else { for (unsigned int q = 0; q < dofs_per_component; ++q) vmult<-1>(&inv_coefficient[q * n_coeff_components], &in[q], &out[q], dofs_per_component, n_desired_components); } } else for (unsigned int q = 0; q < dofs_per_component; ++q) out[q] *= inv_coefficient[q];
for (unsigned int di = 0; di < n_comp_inner; ++di) { Number *out_ = out + di * dofs_per_component; if (dim > 2) evaluator.template hessians<2, false, false>(out_, out_); if (dim > 1) evaluator.template hessians<1, false, false>(out_, out_); evaluator.template hessians<0, false, false>(out_, out_); }
in += dofs_per_component; out += dofs_per_component; inv_coefficient += shift_coefficient; }
return false; }
private: template <int n_components> static inline void vmult(const Number *inverse_coefficients, const Number *src, Number *dst, const unsigned int dofs_per_component, const unsigned int n_given_components = 0) { const unsigned int n_desired_components = (n_components > -1) ? n_components : n_given_components;
std::array<Number, dim + 2> tmp = {}; Assert(n_desired_components <= dim + 2, ExcMessage( "Number of components larger than dim+2 not supported."));
for (unsigned int d = 0; d < n_desired_components; ++d) tmp[d] = src[d * dofs_per_component];
for (unsigned int d1 = 0; d1 < n_desired_components; ++d1) { const Number *inv_coeff_row = &inverse_coefficients[d1 * n_desired_components]; Number sum = inv_coeff_row[0] * tmp[0]; for (unsigned int d2 = 1; d2 < n_desired_components; ++d2) sum += inv_coeff_row[d2] * tmp[d2]; dst[d1 * dofs_per_component] = sum; } } };
/** This struct implements the action of a projection of the values given at the quadrature points to the support points, using an FEEvaluationData argument. For the derivation, see comments in step-67.
Definition at line 2567 of file evaluation_kernels.h.
Static Public Member Functions | |
| template<int fe_degree, int n_q_points_1d> | |
| static bool | run (const unsigned int n_desired_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array) |
|
inlinestatic |
Definition at line 2571 of file evaluation_kernels.h.