const unsigned int first_selected_component = 0);
/**
- * This is one of the main functions in this class and the one that does the
- * heavy lifting.
+ * This function interpolates the finite element solution, represented by
+ * `solution_values` on the given cell, to the `unit_points` passed to the
+ * function.
*
* @param[in] cell An iterator to the current cell in question
*
const EvaluationFlags::EvaluationFlags &evaluation_flags);
/**
- * This is one of the main functions in this class and the one that does the
- * heavy lifting.
+ * This function multiplies the quantities passed in by previous
+ * submit_value() or submit_gradient() calls by the value or gradient of the
+ * test functions, and performs summation over all given points. This is
+ * similar to the integration of a bilinear form in terms of the test
+ * function, with the difference that this formula does not include a `JxW`
+ * factor. This allows the class to naturally embed point information
+ * (e.g. particles) into a finite element formulation. Of course, by
+ * multiplication of a `JxW` information of the data given to
+ * submit_value(), the integration can also be represented by this class.
*
* @param[in] cell An iterator to the current cell in question
*
* @param[in] unit_points List of points in the reference locations of the
* current cell where the FiniteElement object should be evaluated
*
- * @param[out] solution_values This array is filled and can be used to
- * during `cell->set_dof_values(global_vector, solution_values)`.
+ * @param[out] solution_values This array will contain the result of the
+ * integral, which can be used to during
+ * `cell->set_dof_values(solution_values, global_vector)` or
+ * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
+ * that for multi-component systems where only some of the components are
+ * selected by the present class, the entries not touched by this class will
+ * be zeroed out.
*
* @param[in] integration_flags Flags specifying which quantities should be
* integrated at the points.
*/
std::vector<value_type> solution_renumbered;
+ /**
+ * Temporary array to store a vectorized version of the `solution_values`
+ * computed during `integrate()` in a format compatible with the tensor
+ * product evaluators. For vector-valued setups, this array uses a
+ * `Tensor<1, n_components, VectorizedArray<Number>>` format.
+ */
+ AlignedVector<typename internal::FEPointEvaluation::EvaluatorTypeTraits<
+ dim,
+ n_components,
+ VectorizedArray<Number>>::value_type>
+ solution_renumbered_vectorized;
+
/**
* Temporary array to store the values at the points.
*/
poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
poly[1].value(1.) == 1.);
}
- if (true /*TODO: as long as the fast path of integrate() is not working*/)
+ else
{
nonzero_shape_function_component.resize(fe.n_dofs_per_cell());
for (unsigned int d = 0; d < n_components; ++d)
mapping_q_generic->transform_variable(
cell,
mapping_covariant,
+ /* apply_from_left */ true,
unit_points,
ArrayView<const gradient_type>(gradients.data(), gradients.size()),
ArrayView<gradient_type>(gradients.data(), gradients.size()));
}
AssertDimension(solution_values.size(), fe->dofs_per_cell);
- if (false /*TODO*/ && (((integration_flags & EvaluationFlags::values) ||
- (integration_flags & EvaluationFlags::gradients)) &&
- !poly.empty()))
+ if (((integration_flags & EvaluationFlags::values) ||
+ (integration_flags & EvaluationFlags::gradients)) &&
+ !poly.empty())
{
- Assert(false, ExcNotImplemented());
// fast path with tensor product integration
- if (solution_renumbered.size() != dofs_per_component)
- solution_renumbered.resize(dofs_per_component);
-
- // let mapping compute the transformation
+ // let mapping apply the transformation
if (integration_flags & EvaluationFlags::gradients)
{
Assert(mapping_q_generic != nullptr, ExcInternalError());
mapping_q_generic->transform_variable(
cell,
mapping_covariant,
+ /* apply_from_left */ false,
unit_points,
ArrayView<const gradient_type>(gradients.data(), gradients.size()),
ArrayView<gradient_type>(gradients.data(), gradients.size()));
if (integration_flags & EvaluationFlags::gradients)
AssertIndexRange(unit_points.size(), gradients.size() + 1);
+ if (solution_renumbered_vectorized.size() != dofs_per_component)
+ solution_renumbered_vectorized.resize(dofs_per_component);
+ // zero content
+ solution_renumbered_vectorized.fill(
+ typename internal::FEPointEvaluation::EvaluatorTypeTraits<
+ dim,
+ n_components,
+ VectorizedArray<Number>>::value_type());
+
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)
typename internal::ProductTypeNoPoint<value_type,
VectorizedArray<Number>>::type
- value;
+ value = {};
Tensor<1,
dim,
typename internal::ProductTypeNoPoint<
VectorizedArray<Number>>::type>
gradient;
- // convert back to standard format
if (integration_flags & EvaluationFlags::values)
for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
internal::FEPointEvaluation::
gradient, j, gradients[i + j]);
// compute
- internal::integrate_tensor_product_value_and_gradient(
+ internal::integrate_add_tensor_product_value_and_gradient(
poly,
- solution_renumbered,
value,
gradient,
vectorized_points,
- poly.size() == 2);
+ solution_renumbered_vectorized);
}
+ // add between the lanes and write into the result
+ std::fill(solution_values.begin(), solution_values.end(), Number());
for (unsigned int comp = 0; comp < n_components; ++comp)
for (unsigned int i = 0; i < dofs_per_component; ++i)
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::write_value(
- solution_values[renumber[comp * dofs_per_component + i]],
- comp,
- solution_renumbered[i]);
+ {
+ VectorizedArray<Number> result;
+ internal::FEPointEvaluation::
+ EvaluatorTypeTraits<dim, n_components, VectorizedArray<Number>>::
+ write_value(result, comp, solution_renumbered_vectorized[i]);
+ for (unsigned int lane = n_lanes / 2; lane > 0; lane /= 2)
+ for (unsigned int j = 0; j < lane; ++j)
+ result[j] += result[lane + j];
+ solution_values[renumber[comp * dofs_per_component + i]] =
+ result[0];
+ }
}
else if ((integration_flags & EvaluationFlags::values) ||
(integration_flags & EvaluationFlags::gradients))
* As compared to the other transform functions that rely on pre-computed
* information of InternalDataBase, this function chooses the flexible
* evaluation path on the cell and points passed in to the current
- * function. The types `Number` and `Number2` of the input and output arrays
- * must be such that `Number2 = apply_transformation(DerivativeForm<1,
- * spacedim, dim>, Number)`.
+ * function.
+ *
+ * @param cell The cell where to evaluate the mapping
+ *
+ * @param kind Select the kind of the mapping to be applied; currently, this
+ * class only implements `mapping_covariant`.
+ *
+ * @param apply_from_left If true, the mapping is applied to a input
+ * vector from the left, the usual choice for the transformation of field
+ * variables. If false, the mapping is applied to the vector from the right,
+ * representing the transpose of the initial operation; this is the choice
+ * to be used for implementing the action by a test function in an
+ * integration step.
+ *
+ * @param unit_points The points in reference coordinates where the
+ * transformation should be computed and applied to the vector.
+ *
+ * @param input The array of vectors (e.g., gradients for
+ * `mapping_covariant`) to be transformed.
+ *
+ * @param output The array where the result will be stored. The types
+ * `Number` and `Number2` of the input and output arrays must be such that
+ * `Number2 = apply_transformation(DerivativeForm<1, spacedim, dim>,
+ * Number)`. In case the two number types match, this array can be the same
+ * as input array.
*/
template <typename Number, typename Number2>
void
transform_variable(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const MappingKind kind,
+ const bool apply_from_left,
const ArrayView<const Point<dim>> & unit_points,
const ArrayView<const Number> & input,
const ArrayView<Number2> & output) const;
MappingQGeneric<dim, spacedim>::transform_variable(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const MappingKind kind,
+ const bool apply_from_left,
const ArrayView<const Point<dim>> & unit_points,
const ArrayView<const Number> & input,
const ArrayView<Number2> & output) const
renumber_lexicographic_to_hierarchic)
.second;
- const DerivativeForm<1, spacedim, dim, VectorizedArray<double>> jac =
- grad.transpose().covariant_form();
- for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ if (apply_from_left)
{
- DerivativeForm<1, spacedim, dim> jac_j;
- for (unsigned int d = 0; d < spacedim; ++d)
- for (unsigned int e = 0; e < dim; ++e)
- jac_j[d][e] = jac[d][e][j];
- output[i + j] = apply_transformation(jac_j, input[i + j]);
+ const DerivativeForm<1, spacedim, dim, VectorizedArray<double>>
+ jac = grad.transpose().covariant_form();
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ {
+ DerivativeForm<1, spacedim, dim> jac_j;
+ for (unsigned int d = 0; d < spacedim; ++d)
+ for (unsigned int e = 0; e < dim; ++e)
+ jac_j[d][e] = jac[d][e][j];
+ output[i + j] = apply_transformation(jac_j, input[i + j]);
+ }
+ }
+ else
+ {
+ const DerivativeForm<1, dim, spacedim, VectorizedArray<double>>
+ jac = grad.covariant_form();
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ {
+ DerivativeForm<1, dim, spacedim> jac_j;
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int e = 0; e < spacedim; ++e)
+ jac_j[d][e] = jac[d][e][j];
+ output[i + j] = apply_transformation(jac_j, input[i + j]);
+ }
}
}
else
polynomial_degree == 1,
renumber_lexicographic_to_hierarchic)
.second;
- output[i] =
- apply_transformation(grad.transpose().covariant_form(), input[i]);
+ if (apply_from_left)
+ output[i] =
+ apply_transformation(grad.transpose().covariant_form(), input[i]);
+ else
+ output[i] = apply_transformation(grad.covariant_form(), input[i]);
}
}
*/
template <int dim, typename Number, typename Number2>
inline void
- integrate_tensor_product_value_and_gradient(
- const std::vector<Polynomials::Polynomial<double>> & poly,
- const std::vector<Number> & values,
- const typename ProductTypeNoPoint<Number, Number2>::type &value,
- const Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>
- & gradient,
- const Point<dim, Number2> & p,
- const bool d_linear = false,
- const std::vector<unsigned int> &renumber = {})
+ 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 = {})
{
- Assert(false, ExcNotImplemented());
-
- (void)poly;
- (void)values;
- (void)value;
- (void)gradient;
- (void)p;
- (void)d_linear;
- (void)renumber;
+ static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
+
+ const unsigned 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);
+ std::array<Number, 2 * dim * 200> shapes;
+
+ // Evaluate 1D polynomials and their derivatives
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int i = 0; i < n_shapes; ++i)
+ poly[i].value(p[d], 1, shapes.data() + 2 * (d * n_shapes + i));
+
+ // Implement the transpose of the function above
+ for (unsigned int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
+ {
+ const Number2 test_value_z =
+ dim > 2 ? (value * shapes[4 * n_shapes + 2 * i2] +
+ gradient[2] * shapes[4 * n_shapes + 2 * i2 + 1]) :
+ value;
+ const Number2 test_grad_x =
+ dim > 2 ? gradient[0] * shapes[4 * n_shapes + 2 * i2] : gradient[0];
+ const Number2 test_grad_y =
+ dim > 2 ? gradient[1] * shapes[4 * n_shapes + 2 * i2] :
+ (dim > 1 ? gradient[1] : Number2());
+ for (unsigned int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
+ {
+ const Number2 test_value_y =
+ dim > 1 ? (test_value_z * shapes[2 * n_shapes + 2 * i1] +
+ test_grad_y * shapes[2 * n_shapes + 2 * i1 + 1]) :
+ test_value_z;
+ const Number2 test_grad_xy =
+ dim > 1 ? test_grad_x * shapes[2 * n_shapes + 2 * i1] :
+ test_grad_x;
+ if (renumber.empty())
+ for (unsigned int i0 = 0; i0 < n_shapes; ++i0, ++i)
+ values[i] += shapes[2 * i0] * test_value_y +
+ shapes[2 * i0 + 1] * test_grad_xy;
+ else
+ for (unsigned int i0 = 0; i0 < n_shapes; ++i0, ++i)
+ values[renumber[i]] += shapes[2 * i0] * test_value_y +
+ shapes[2 * i0 + 1] * test_grad_xy;
+ }
+ }
}
projects / dealii.git / commitdiff