void
setup(const unsigned int first_selected_component);
+ /**
+ * Shared functionality of all @p reinit() functions. Resizes data fields and
+ * precomputes the @p shapes vector, holding the evaluation of 1D basis
+ * functions of tensor product polynomials, if necessary.
+ */
+ void
+ do_reinit();
+
/**
* Number of quadrature points of the current cell/face.
*/
* Bool indicating if class is reinitialized and data vectors a resized.
*/
bool is_reinitialized;
+
+ /**
+ * Vector containing tensor product shape functions evaluated (during
+ * reinit()) at the vectorized unit points.
+ */
+ AlignedVector<dealii::ndarray<VectorizedArray<Number>, 2, dim>> shapes;
};
// ----------------------- template and inline function ----------------------
AssertIndexRange(first_selected_component + n_components,
fe->n_components() + 1);
+ shapes.reserve(100);
+
bool same_base_element = true;
unsigned int base_element_number = 0;
component_in_base_element = 0;
fe_values->reinit(cell);
}
- const_cast<unsigned int &>(n_q_points) = unit_points.size();
-
- if (update_flags & update_values)
- values.resize(n_q_points, numbers::signaling_nan<value_type>());
- if (update_flags & update_gradients)
- gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
-
- is_reinitialized = true;
+ do_reinit();
}
current_cell_index = cell_index;
current_face_number = numbers::invalid_unsigned_int;
- const_cast<unsigned int &>(n_q_points) =
- mapping_info->get_unit_points(current_cell_index, current_face_number)
- .size();
-
- if (update_flags & update_values)
- values.resize(n_q_points, numbers::signaling_nan<value_type>());
- if (update_flags & update_gradients)
- gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
-
- is_reinitialized = true;
+ do_reinit();
}
current_cell_index = cell_index;
current_face_number = face_number;
- const_cast<unsigned int &>(n_q_points) =
- mapping_info->get_unit_points(current_cell_index, current_face_number)
- .size();
+ do_reinit();
+}
+
+
+
+template <int n_components, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components, dim, spacedim, Number>::do_reinit()
+{
+ const auto unit_points =
+ mapping_info->get_unit_points(current_cell_index, current_face_number);
+
+ const_cast<unsigned int &>(n_q_points) = unit_points.size();
if (update_flags & update_values)
values.resize(n_q_points, numbers::signaling_nan<value_type>());
if (update_flags & update_gradients)
gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
+ if (!polynomials_are_hat_functions)
+ {
+ const std::size_t n_points = unit_points.size();
+ const std::size_t n_lanes = VectorizedArray<Number>::size();
+ const std::size_t n_batches =
+ n_points / n_lanes + (n_points % n_lanes > 0 ? 1 : 0);
+ const std::size_t n_shapes = poly.size();
+ shapes.resize_fast(n_batches * n_shapes);
+ for (unsigned int i = 0, qb = 0; i < n_points; i += n_lanes, ++qb)
+ {
+ // convert to vectorized format
+ Point<dim, VectorizedArray<Number>> vectorized_points;
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ for (unsigned int d = 0; d < dim; ++d)
+ vectorized_points[d][j] = unit_points[i + j][d];
+
+ auto view =
+ make_array_view(shapes.begin() + qb * n_shapes,
+ shapes.begin() + (qb * n_shapes + n_shapes));
+
+ internal::compute_values_of_array(view, poly, vectorized_points);
+ }
+ }
+
is_reinitialized = true;
}
const ArrayView<const Number> & solution_values,
const EvaluationFlags::EvaluationFlags &evaluation_flag)
{
+ if (!is_reinitialized)
+ reinit(numbers::invalid_unsigned_int);
+
if (n_q_points == 0)
return;
fast_path)
{
// fast path with tensor product evaluation
- const auto unit_points =
- mapping_info->get_unit_points(current_cell_index, current_face_number);
-
- // we need to call reinit() here if we reuse the same MappingInfo object
- // for several FEPointEvaluation objects to resize the data fields
- if (!is_reinitialized || n_q_points != unit_points.size())
- reinit(numbers::invalid_unsigned_int);
-
if (solution_renumbered.size() != dofs_per_component)
solution_renumbered.resize(dofs_per_component);
for (unsigned int comp = 0; comp < n_components; ++comp)
unit_gradients.resize(n_q_points,
numbers::signaling_nan<gradient_type>());
+ const auto unit_points =
+ mapping_info->get_unit_points(current_cell_index, current_face_number);
+ const auto &mapping_data =
+ mapping_info->get_mapping_data(current_cell_index, current_face_number);
+
const std::size_t n_points = unit_points.size();
const std::size_t n_lanes = VectorizedArray<Number>::size();
- for (unsigned int i = 0; i < n_points; i += n_lanes)
+ for (unsigned int i = 0, qb = 0; i < n_points; i += n_lanes, ++qb)
{
- // convert to vectorized format
- Point<dim, VectorizedArray<Number>> vectorized_points;
- for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
- for (unsigned int d = 0; d < dim; ++d)
- vectorized_points[d][j] = unit_points[i + j][d];
-
// compute
- const auto val_and_grad =
- internal::evaluate_tensor_product_value_and_gradient(
- poly,
- solution_renumbered,
- vectorized_points,
- polynomials_are_hat_functions);
+ const unsigned int n_shapes = poly.size();
+ const auto val_and_grad = [&]() {
+ if (polynomials_are_hat_functions)
+ {
+ // convert to vectorized format
+ Point<dim, VectorizedArray<Number>> vectorized_points;
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ for (unsigned int d = 0; d < dim; ++d)
+ vectorized_points[d][j] = unit_points[i + j][d];
+
+ return internal::
+ evaluate_tensor_product_value_and_gradient_linear(
+ poly, solution_renumbered, vectorized_points);
+ }
+ else
+ return internal::
+ evaluate_tensor_product_value_and_gradient_shapes<
+ dim,
+ value_type,
+ VectorizedArray<Number>>(
+ make_array_view(shapes.begin() + qb * n_shapes,
+ shapes.begin() + (qb * n_shapes + n_shapes)),
+ poly.size(),
+ solution_renumbered);
+ }();
// convert back to standard format
if (evaluation_flag & EvaluationFlags::values)
Number>::set_gradient(val_and_grad.second,
j,
unit_gradients[i + j]);
- const auto &mapping_data =
- mapping_info->get_mapping_data(current_cell_index,
- current_face_number);
gradients[i + j] = apply_transformation(
mapping_data.inverse_jacobians[i + j].transpose(),
unit_gradients[i + j]);
fast_path)
{
// fast path with tensor product integration
- const auto unit_points =
- mapping_info->get_unit_points(current_cell_index, current_face_number);
-
- // we need to call reinit() here if we reuse the same MappingInfo object
- // for several FEPointEvaluation objects to resize the data fields
- if (!is_reinitialized || n_q_points != unit_points.size())
- reinit(numbers::invalid_unsigned_int);
if (integration_flags & EvaluationFlags::values)
AssertIndexRange(n_q_points, values.size() + 1);
n_components,
VectorizedArray<Number>>::value_type());
+ const auto unit_points =
+ mapping_info->get_unit_points(current_cell_index, current_face_number);
+ const auto &mapping_data =
+ mapping_info->get_mapping_data(current_cell_index, current_face_number);
+
const std::size_t n_points = unit_points.size();
const std::size_t n_lanes = VectorizedArray<Number>::size();
- for (unsigned int i = 0; i < n_points; i += n_lanes)
+ for (unsigned int i = 0, qb = 0; i < n_points; i += n_lanes, ++qb)
{
- // convert to vectorized format
- Point<dim, VectorizedArray<Number>> vectorized_points;
- for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
- for (unsigned int d = 0; d < dim; ++d)
- vectorized_points[d][j] = unit_points[i + j][d];
-
typename internal::ProductTypeNoPoint<value_type,
VectorizedArray<Number>>::type
value = {};
if (integration_flags & EvaluationFlags::gradients)
for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
{
- const auto &mapping_data =
- mapping_info->get_mapping_data(current_cell_index,
- current_face_number);
gradients[i + j] =
apply_transformation(mapping_data.inverse_jacobians[i + j],
gradients[i + j]);
}
// compute
- internal::integrate_add_tensor_product_value_and_gradient(
- poly,
- value,
- gradient,
- vectorized_points,
- solution_renumbered_vectorized);
+ const unsigned int n_shapes = poly.size();
+ if (polynomials_are_hat_functions)
+ {
+ // convert to vectorized format
+ Point<dim, VectorizedArray<Number>> vectorized_points;
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ for (unsigned int d = 0; d < dim; ++d)
+ vectorized_points[d][j] = unit_points[i + j][d];
+
+ internal::integrate_add_tensor_product_value_and_gradient_linear(
+ poly,
+ value,
+ gradient,
+ solution_renumbered_vectorized,
+ vectorized_points);
+ }
+ else
+ internal::integrate_add_tensor_product_value_and_gradient_shapes<
+ dim,
+ VectorizedArray<Number>,
+ typename internal::
+ ProductTypeNoPoint<value_type, VectorizedArray<Number>>::type>(
+ make_array_view(shapes.begin() + qb * n_shapes,
+ shapes.begin() + (qb * n_shapes + n_shapes)),
+ n_shapes,
+ value,
+ gradient,
+ solution_renumbered_vectorized);
}
// add between the lanes and write into the result
+ /**
+ * Computes the values and derivatives of the 1d polynomials @p poly at the
+ * specified point @p p and stores it in @p shapes.
+ */
+ template <int dim, typename Number>
+ inline void
+ compute_values_of_array(
+ ArrayView<dealii::ndarray<Number, 2, dim>> & shapes,
+ const std::vector<Polynomials::Polynomial<double>> &poly,
+ const Point<dim, Number> & p)
+ {
+ const int n_shapes = poly.size();
+
+ // Evaluate 1d polynomials and their derivatives
+ std::array<Number, dim> point;
+ for (unsigned int d = 0; d < dim; ++d)
+ point[d] = p[d];
+ for (int i = 0; i < n_shapes; ++i)
+ poly[i].values_of_array(point, 1, &shapes[i][0]);
+ }
+
+
+
/**
* Interpolate inner dimensions of tensor product shape functions.
*/
template <int dim, int length, typename Number2, typename Number>
inline std::array<typename ProductTypeNoPoint<Number, Number2>::type, 3>
- do_interpolate_xy(const std::vector<Number> & values,
- const std::vector<unsigned int> & renumber,
- const dealii::ndarray<Number2, 200, 2, dim> &shapes,
+ do_interpolate_xy(const std::vector<Number> & values,
+ const std::vector<unsigned int> & renumber,
+ const ArrayView<dealii::ndarray<Number2, 2, dim>> &shapes,
const int n_shapes_runtime,
int & i)
{
/**
- * Compute the polynomial interpolation of a tensor product shape function
- * $\varphi_i$ given a vector of coefficients $u_i$ in the form
- * $u_h(\mathbf{x}) = \sum_{i=1}^{k^d} \varphi_i(\mathbf{x}) u_i$. The shape
- * functions $\varphi_i(\mathbf{x}) =
- * \prod_{d=1}^{\text{dim}}\varphi_{i_d}^\text{1d}(x_d)$ represent a tensor
- * product. The function returns a pair with the value of the interpolation
- * as the first component and the gradient in reference coordinates as the
- * second component. Note that for compound types (e.g. the `values` field
- * begin a Point<spacedim> argument), the components of the gradient are
- * sorted as Tensor<1, dim, Tensor<1, spacedim>> with the derivatives
- * as the first index; this is a consequence of the generic arguments in the
- * function.
- *
- * @param poly The underlying one-dimensional polynomial basis
- * $\{\varphi^{1d}_{i_1}\}$ given as a vector of polynomials.
- *
- * @param values The expansion coefficients $u_i$ of type `Number` in
- * the polynomial interpolation. The coefficients can be simply `double`
- * variables but e.g. also Point<spacedim> in case they define arithmetic
- * operations with the type `Number2`.
- *
- * @param p The position in reference coordinates where the interpolation
- * should be evaluated.
- *
- * @param d_linear Flag to specify whether a d-linear (linear in 1d,
- * bi-linear in 2d, tri-linear in 3d) interpolation should be made, which
- * allows to unroll loops and considerably speed up evaluation.
- *
- * @param renumber Optional parameter to specify a renumbering in the
- * coefficient vector, assuming that `values[renumber[i]]` returns
- * the lexicographic (tensor product) entry of the coefficients. If the
- * vector is entry, the values are assumed to be sorted lexicographically.
+ * Interpolates the values and gradients into the points specified in
+ * @p compute_values_of_array() with help of the precomputed @p shapes.
*/
template <int dim, typename Number, typename Number2>
inline std::pair<
typename ProductTypeNoPoint<Number, Number2>::type,
Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
- evaluate_tensor_product_value_and_gradient(
- const std::vector<Polynomials::Polynomial<double>> &poly,
- const std::vector<Number> & values,
- const Point<dim, Number2> & p,
- const bool d_linear = false,
- const std::vector<unsigned int> & renumber = {})
+ evaluate_tensor_product_value_and_gradient_shapes(
+ const ArrayView<dealii::ndarray<Number2, 2, dim>> &shapes,
+ const int n_shapes,
+ const std::vector<Number> & values,
+ const std::vector<unsigned int> & renumber = {})
{
static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
- // use `int` type for this variable and the loops below to inform the
- // compiler that the loops below will never overflow, which allows it to
- // generate more optimized code for the variable loop bounds in the
- // present context
- const int n_shapes = poly.size();
AssertDimension(Utilities::pow(n_shapes, dim), values.size());
Assert(renumber.empty() || renumber.size() == values.size(),
ExcDimensionMismatch(renumber.size(), values.size()));
- // shortcut for linear interpolation to speed up evaluation
- if (d_linear)
- {
- AssertDimension(poly.size(), 2);
- for (unsigned int i = 0; i < renumber.size(); ++i)
- AssertDimension(renumber[i], i);
-
- if (dim == 1)
- {
- Tensor<1, dim, Number3> derivative;
- derivative[0] = values[1] - values[0];
- return std::make_pair((1. - p[0]) * values[0] + p[0] * values[1],
- derivative);
- }
- else if (dim == 2)
- {
- const Number2 x0 = 1. - p[0], x1 = p[0];
- const Number3 tmp0 = x0 * values[0] + x1 * values[1];
- const Number3 tmp1 = x0 * values[2] + x1 * values[3];
- const Number3 mapped = (1. - p[1]) * tmp0 + p[1] * tmp1;
- Tensor<1, dim, Number3> derivative;
- derivative[0] = (1. - p[1]) * (values[1] - values[0]) +
- p[1] * (values[3] - values[2]);
- derivative[1] = tmp1 - tmp0;
- return std::make_pair(mapped, derivative);
- }
- else if (dim == 3)
- {
- const Number2 x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
- z0 = 1. - p[2], z1 = p[2];
- const Number3 tmp0 = x0 * values[0] + x1 * values[1];
- const Number3 tmp1 = x0 * values[2] + x1 * values[3];
- const Number3 tmpy0 = y0 * tmp0 + y1 * tmp1;
- const Number3 tmp2 = x0 * values[4] + x1 * values[5];
- const Number3 tmp3 = x0 * values[6] + x1 * values[7];
- const Number3 tmpy1 = y0 * tmp2 + y1 * tmp3;
- const Number3 mapped = z0 * tmpy0 + z1 * tmpy1;
- Tensor<1, dim, Number3> derivative;
- derivative[2] = tmpy1 - tmpy0;
- derivative[1] = z0 * (tmp1 - tmp0) + z1 * (tmp3 - tmp2);
- derivative[0] =
- z0 *
- (y0 * (values[1] - values[0]) + y1 * (values[3] - values[2])) +
- z1 *
- (y0 * (values[5] - values[4]) + y1 * (values[7] - values[6]));
- return std::make_pair(mapped, derivative);
- }
- }
-
- AssertIndexRange(n_shapes, 200);
- dealii::ndarray<Number2, 200, 2, dim> shapes;
-
- // Evaluate 1d polynomials and their derivatives
- std::array<Number2, dim> point;
- for (unsigned int d = 0; d < dim; ++d)
- point[d] = p[d];
- for (int i = 0; i < n_shapes; ++i)
- poly[i].values_of_array(point, 1, &shapes[i][0]);
-
// Go through the tensor product of shape functions and interpolate
// with optimal algorithm
std::pair<Number3, Tensor<1, dim, Number3>> result = {};
+ /**
+ * Specializes @p evaluate_tensor_product_value_and_gradient() for linear
+ * polynomials which massively reduces the necessary instructions.
+ */
+ template <int dim, typename Number, typename Number2>
+ inline std::pair<
+ typename ProductTypeNoPoint<Number, Number2>::type,
+ Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
+ evaluate_tensor_product_value_and_gradient_linear(
+ const std::vector<Polynomials::Polynomial<double>> &poly,
+ const std::vector<Number> & values,
+ const Point<dim, Number2> & p,
+ const std::vector<unsigned int> & renumber = {})
+ {
+ (void)poly;
+ static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
+
+ using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+
+ AssertDimension(Utilities::pow(poly.size(), dim), values.size());
+ Assert(renumber.empty() || renumber.size() == values.size(),
+ ExcDimensionMismatch(renumber.size(), values.size()));
+
+ AssertDimension(poly.size(), 2);
+ for (unsigned int i = 0; i < renumber.size(); ++i)
+ AssertDimension(renumber[i], i);
+
+ if (dim == 1)
+ {
+ Tensor<1, dim, Number3> derivative;
+ derivative[0] = values[1] - values[0];
+ return std::make_pair((1. - p[0]) * values[0] + p[0] * values[1],
+ derivative);
+ }
+ else if (dim == 2)
+ {
+ const Number2 x0 = 1. - p[0], x1 = p[0];
+ const Number3 tmp0 = x0 * values[0] + x1 * values[1];
+ const Number3 tmp1 = x0 * values[2] + x1 * values[3];
+ const Number3 mapped = (1. - p[1]) * tmp0 + p[1] * tmp1;
+ Tensor<1, dim, Number3> derivative;
+ derivative[0] = (1. - p[1]) * (values[1] - values[0]) +
+ p[1] * (values[3] - values[2]);
+ derivative[1] = tmp1 - tmp0;
+ return std::make_pair(mapped, derivative);
+ }
+ else if (dim == 3)
+ {
+ const Number2 x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
+ z0 = 1. - p[2], z1 = p[2];
+ const Number3 tmp0 = x0 * values[0] + x1 * values[1];
+ const Number3 tmp1 = x0 * values[2] + x1 * values[3];
+ const Number3 tmpy0 = y0 * tmp0 + y1 * tmp1;
+ const Number3 tmp2 = x0 * values[4] + x1 * values[5];
+ const Number3 tmp3 = x0 * values[6] + x1 * values[7];
+ const Number3 tmpy1 = y0 * tmp2 + y1 * tmp3;
+ const Number3 mapped = z0 * tmpy0 + z1 * tmpy1;
+ Tensor<1, dim, Number3> derivative;
+ derivative[2] = tmpy1 - tmpy0;
+ derivative[1] = z0 * (tmp1 - tmp0) + z1 * (tmp3 - tmp2);
+ derivative[0] =
+ z0 * (y0 * (values[1] - values[0]) + y1 * (values[3] - values[2])) +
+ z1 * (y0 * (values[5] - values[4]) + y1 * (values[7] - values[6]));
+ return std::make_pair(mapped, derivative);
+ }
+
+ // work around a compile error: missing return statement
+ return std::make_pair(Number3(), Tensor<1, dim, Number3>());
+ }
+
+
+
+ /**
+ * Compute the polynomial interpolation of a tensor product shape function
+ * $\varphi_i$ given a vector of coefficients $u_i$ in the form
+ * $u_h(\mathbf{x}) = \sum_{i=1}^{k^d} \varphi_i(\mathbf{x}) u_i$. The shape
+ * functions $\varphi_i(\mathbf{x}) =
+ * \prod_{d=1}^{\text{dim}}\varphi_{i_d}^\text{1d}(x_d)$ represent a tensor
+ * product. The function returns a pair with the value of the interpolation
+ * as the first component and the gradient in reference coordinates as the
+ * second component. Note that for compound types (e.g. the `values` field
+ * begin a Point<spacedim> argument), the components of the gradient are
+ * sorted as Tensor<1, dim, Tensor<1, spacedim>> with the derivatives
+ * as the first index; this is a consequence of the generic arguments in the
+ * function.
+ *
+ * @param poly The underlying one-dimensional polynomial basis
+ * $\{\varphi^{1d}_{i_1}\}$ given as a vector of polynomials.
+ *
+ * @param values The expansion coefficients $u_i$ of type `Number` in
+ * the polynomial interpolation. The coefficients can be simply `double`
+ * variables but e.g. also Point<spacedim> in case they define arithmetic
+ * operations with the type `Number2`.
+ *
+ * @param p The position in reference coordinates where the interpolation
+ * should be evaluated.
+ *
+ * @param d_linear Flag to specify whether a d-linear (linear in 1d,
+ * bi-linear in 2d, tri-linear in 3d) interpolation should be made, which
+ * allows to unroll loops and considerably speed up evaluation.
+ *
+ * @param renumber Optional parameter to specify a renumbering in the
+ * coefficient vector, assuming that `values[renumber[i]]` returns
+ * the lexicographic (tensor product) entry of the coefficients. If the
+ * vector is entry, the values are assumed to be sorted lexicographically.
+ */
+ template <int dim, typename Number, typename Number2>
+ inline std::pair<
+ typename ProductTypeNoPoint<Number, Number2>::type,
+ Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
+ evaluate_tensor_product_value_and_gradient(
+ const std::vector<Polynomials::Polynomial<double>> &poly,
+ const std::vector<Number> & values,
+ const Point<dim, Number2> & p,
+ const bool d_linear = false,
+ const std::vector<unsigned int> & renumber = {})
+ {
+ if (d_linear)
+ {
+ return evaluate_tensor_product_value_and_gradient_linear(poly,
+ values,
+ p,
+ renumber);
+ }
+ else
+ {
+ std::array<dealii::ndarray<Number2, 2, dim>, 200> shapes;
+
+ auto view = make_array_view(shapes);
+
+ compute_values_of_array(view, poly, p);
+
+ return evaluate_tensor_product_value_and_gradient_shapes<dim,
+ Number,
+ Number2>(
+ view, poly.size(), values, renumber);
+ }
+ }
+
+
+
template <int dim, typename Number, typename Number2>
SymmetricTensor<2, dim, typename ProductTypeNoPoint<Number, Number2>::type>
evaluate_tensor_product_hessian(
*/
template <int dim, int length, typename Number2, typename Number>
inline void
- do_apply_test_functions_xy(AlignedVector<Number2> & values,
- const std::vector<unsigned int> &renumber,
- const dealii::ndarray<Number, 200, 2, dim> &shapes,
- const std::array<Number2, 3> &test_grads_value,
- const int n_shapes_runtime,
- int & i)
+ do_apply_test_functions_xy(
+ AlignedVector<Number2> & values,
+ const std::vector<unsigned int> & renumber,
+ const ArrayView<dealii::ndarray<Number, 2, dim>> &shapes,
+ const std::array<Number2, 3> & test_grads_value,
+ const int n_shapes_runtime,
+ int & i)
{
const int n_shapes = length > 0 ? length : n_shapes_runtime;
for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
*/
template <int dim, typename Number, typename Number2>
inline void
- integrate_add_tensor_product_value_and_gradient(
- const std::vector<Polynomials::Polynomial<double>> &poly,
- const Number2 & value,
- const Tensor<1, dim, Number2> & gradient,
- const Point<dim, Number> & p,
- AlignedVector<Number2> & values,
- const std::vector<unsigned int> & renumber = {})
+ integrate_add_tensor_product_value_and_gradient_shapes(
+ const ArrayView<dealii::ndarray<Number, 2, dim>> &shapes,
+ const int n_shapes,
+ const Number2 & value,
+ const Tensor<1, dim, Number2> & gradient,
+ AlignedVector<Number2> & values,
+ const std::vector<unsigned int> & renumber = {})
{
static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
// as in evaluate, use `int` type to produce better code in this context
- const int n_shapes = poly.size();
AssertDimension(Utilities::pow(n_shapes, dim), values.size());
Assert(renumber.empty() || renumber.size() == values.size(),
ExcDimensionMismatch(renumber.size(), values.size()));
- AssertIndexRange(n_shapes, 200);
- dealii::ndarray<Number, 200, 2, dim> shapes;
-
- // Evaluate 1d polynomials and their derivatives
- std::array<Number, dim> point;
- for (unsigned int d = 0; d < dim; ++d)
- point[d] = p[d];
- for (int i = 0; i < n_shapes; ++i)
- poly[i].values_of_array(point, 1, &shapes[i][0]);
-
// Implement the transpose of the function above
std::array<Number2, 3> test_grads_value;
for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
+ /**
+ * Specializes @p evaluate_tensor_product_value_and_gradient() for linear
+ * polynomials which massively reduces the necessary instructions.
+ */
+ template <int dim, typename Number, typename Number2>
+ inline void
+ integrate_add_tensor_product_value_and_gradient_linear(
+ const std::vector<Polynomials::Polynomial<double>> &poly,
+ const Number2 & value,
+ const Tensor<1, dim, Number2> & gradient,
+ AlignedVector<Number2> & values,
+ const Point<dim, Number> & p,
+ const std::vector<unsigned int> & renumber = {})
+ {
+ (void)poly;
+ static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
+
+ AssertDimension(Utilities::pow(poly.size(), dim), values.size());
+ Assert(renumber.empty() || renumber.size() == values.size(),
+ ExcDimensionMismatch(renumber.size(), values.size()));
+
+ AssertDimension(poly.size(), 2);
+ for (unsigned int i = 0; i < renumber.size(); ++i)
+ AssertDimension(renumber[i], i);
+
+ if (dim == 1)
+ {
+ const auto x0 = 1. - p[0], x1 = p[0];
+
+ values[0] = value * x0 - gradient[0];
+ values[1] = value * x1 + gradient[0];
+ }
+ else if (dim == 2)
+ {
+ const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1];
+
+ const auto test_value_y0 = value * y0 - gradient[1];
+ const auto test_grad_xy0 = gradient[0] * y0;
+ const auto test_value_y1 = value * y1 + gradient[1];
+ const auto test_grad_xy1 = gradient[0] * y1;
+
+ values[0] += x0 * test_value_y0 - test_grad_xy0;
+ values[1] += x1 * test_value_y0 + test_grad_xy0;
+ values[2] += x0 * test_value_y1 - test_grad_xy1;
+ values[3] += x1 * test_value_y1 + test_grad_xy1;
+ }
+ else if (dim == 3)
+ {
+ const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
+ z0 = 1. - p[2], z1 = p[2];
+
+ const auto test_value_z0 = value * z0 - gradient[2];
+ const auto test_grad_x0 = gradient[0] * z0;
+ const auto test_grad_y0 = gradient[1] * z0;
+ const auto test_value_z1 = value * z1 + gradient[2];
+ const auto test_grad_x1 = gradient[0] * z1;
+ const auto test_grad_y1 = gradient[1] * z1;
+
+ const auto test_value_y00 = test_value_z0 * y0 - test_grad_y0;
+ const auto test_grad_xy00 = test_grad_x0 * y0;
+ const auto test_value_y01 = test_value_z0 * y1 + test_grad_y0;
+ const auto test_grad_xy01 = test_grad_x0 * y1;
+ const auto test_value_y10 = test_value_z1 * y0 - test_grad_y1;
+ const auto test_grad_xy10 = test_grad_x1 * y0;
+ const auto test_value_y11 = test_value_z1 * y1 + test_grad_y1;
+ const auto test_grad_xy11 = test_grad_x1 * y1;
+
+ values[0] += x0 * test_value_y00 - test_grad_xy00;
+ values[1] += x1 * test_value_y00 + test_grad_xy00;
+ values[2] += x0 * test_value_y01 - test_grad_xy01;
+ values[3] += x1 * test_value_y01 + test_grad_xy01;
+ values[4] += x0 * test_value_y10 - test_grad_xy10;
+ values[5] += x1 * test_value_y10 + test_grad_xy10;
+ values[6] += x0 * test_value_y11 - test_grad_xy11;
+ values[7] += x1 * test_value_y11 + test_grad_xy11;
+ }
+ }
+
+
+
template <int dim, int n_points_1d_template, typename Number>
inline void
weight_fe_q_dofs_by_entity(const Number * weights,