#include <deal.II/fe/mapping_q_generic.h>
+#include <deal.II/matrix_free/evaluation_flags.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/tensor_product_kernels.h>
* values on the element as returned by `cell->get_dof_values(global_vector,
* solution_values)`.
*
- * @param[out] values Array where the values of the evaluation of the
- * solution interpolated to the given points is stored. The array must have
- * the same size as the `unit_points` argument if values should be
- * computed. If empty, values will not be computed.
- *
- * @param[out] gradients Array where the gradients of the evaluation of the
- * solution interpolated to the given points is stored. The gradients are
- * expressed in real space via the Mapping defined in the constructor. The
- * array must have the same size as the `unit_points` argument if gradients
- * should be computed. If empty, gradients will not be computed.
+ * @param[in] flags Flags specifying which quantities should be evaluated
+ * at the points.
*/
void
evaluate(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const ArrayView<const Point<dim>> &unit_points,
- const ArrayView<const double> & solution_values,
- const ArrayView<value_type> & values,
- const ArrayView<gradient_type> & gradients = {});
+ const ArrayView<const Point<dim>> & unit_points,
+ const ArrayView<const double> & solution_values,
+ const EvaluationFlags::EvaluationFlags &flags);
+
+ /**
+ * Return the value at quadrature point number @p q_point after a call to
+ * FEPointEvaluation::evaluate() with EvaluationFlags::value set, or
+ * the value that has been stored there with a call to
+ * FEPointEvaluation::submit_value(). If the object is vector-valued, a
+ * vector-valued return argument is given. Note that when
+ * vectorization is enabled, values from several points are grouped together.
+ */
+ const value_type &
+ get_value(const unsigned int q_point) const;
+
+ /**
+ * Return the gradient at quadrature point number @p q_point after a call to
+ * FEPointEvaluation::evaluate() with EvaluationFlags::gradient set, or
+ * the gradient that has been stored there with a call to
+ * FEPointEvaluation::submit_gradient(). If the object is vector-valued, a
+ * vector-valued return argument is given. Note that when
+ * vectorization is enabled, values from several points are grouped together.
+ */
+ const gradient_type &
+ get_gradient(const unsigned int q_point) const;
private:
/**
*/
std::vector<value_type> solution_renumbered;
+ /**
+ * Temporary array to store the values at the points.
+ */
+ std::vector<value_type> values;
+
+ /**
+ * Temporary array to store the gradients at the points.
+ */
+ std::vector<gradient_type> gradients;
+
/**
* Number of unknowns per component, i.e., number of unique basis functions,
* for the chosen FiniteElement (or base element).
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const ArrayView<const Point<dim>> & unit_points,
const ArrayView<const double> & solution_values,
- const ArrayView<value_type> & values,
- const ArrayView<gradient_type> & gradients)
+ const EvaluationFlags::EvaluationFlags & evaluation_flag)
{
AssertDimension(solution_values.size(), fe->dofs_per_cell);
- if ((values.size() > 0 || gradients.size() > 0) && !poly.empty())
+ if (((evaluation_flag & EvaluationFlags::values) ||
+ (evaluation_flag & EvaluationFlags::gradients)) &&
+ !poly.empty())
{
// fast path with tensor product evaluation
if (solution_renumbered.size() != dofs_per_component)
comp,
solution_renumbered[i]);
- if (values.size() > 0)
- AssertDimension(unit_points.size(), values.size());
- if (gradients.size() > 0)
- AssertDimension(unit_points.size(), gradients.size());
+ if (evaluation_flag & EvaluationFlags::values)
+ values.resize(unit_points.size());
+ if (evaluation_flag & EvaluationFlags::gradients)
+ gradients.resize(unit_points.size());
const std::size_t n_points = unit_points.size();
const std::size_t n_lanes = VectorizedArray<double>::size();
poly, solution_renumbered, vectorized_points, poly.size() == 2);
// convert back to standard format
- if (values.size() > 0)
+ if (evaluation_flag & EvaluationFlags::values)
for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
internal::FEPointEvaluation::EvaluatorTypeTraits<
dim,
n_components>::set_value(val_and_grad.first, j, values[i + j]);
- if (gradients.size() > 0)
+ if (evaluation_flag & EvaluationFlags::gradients)
for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
internal::FEPointEvaluation::
EvaluatorTypeTraits<dim, n_components>::set_gradient(
}
// let mapping compute the transformation
- if (gradients.size() > 0)
+ if (evaluation_flag & EvaluationFlags::gradients)
{
Assert(mapping_q_generic != nullptr, ExcInternalError());
mapping_q_generic->transform_variable(
cell,
mapping_covariant,
unit_points,
- ArrayView<const gradient_type>(gradients.begin(), gradients.size()),
- gradients);
+ ArrayView<const gradient_type>(gradients.data(), gradients.size()),
+ ArrayView<gradient_type>(gradients.data(), gradients.size()));
}
}
- else if (values.size() > 0 || gradients.size() > 0)
+ else if ((evaluation_flag & EvaluationFlags::values) ||
+ (evaluation_flag & EvaluationFlags::gradients))
{
// slow path with FEValues
const UpdateFlags flags =
- (values.size() > 0 ? update_values : update_default) |
- (gradients.size() > 0 ? update_gradients : update_default);
+ ((evaluation_flag & EvaluationFlags::values) ? update_values :
+ update_default) |
+ ((evaluation_flag & EvaluationFlags::gradients) ? update_gradients :
+ update_default);
FEValues<dim, spacedim> fe_values(
*mapping,
*fe,
flags);
fe_values.reinit(cell);
- if (values.size() > 0)
+ if (evaluation_flag & EvaluationFlags::values)
{
- AssertDimension(unit_points.size(), values.size());
+ values.resize(unit_points.size());
std::fill(values.begin(), values.end(), value_type());
for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
{
}
}
- if (gradients.size() > 0)
+ if (evaluation_flag & EvaluationFlags::gradients)
{
- AssertDimension(unit_points.size(), gradients.size());
+ gradients.resize(unit_points.size());
std::fill(gradients.begin(), gradients.end(), gradient_type());
for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
{
}
}
+
+
+template <int n_components, int dim, int spacedim>
+inline const typename FEPointEvaluation<n_components, dim, spacedim>::value_type
+ &
+ FEPointEvaluation<n_components, dim, spacedim>::get_value(
+ const unsigned int q_point) const
+{
+ AssertIndexRange(q_point, values.size());
+ return values[q_point];
+}
+
+
+
+template <int n_components, int dim, int spacedim>
+inline const typename FEPointEvaluation<n_components, dim, spacedim>::
+ gradient_type &
+ FEPointEvaluation<n_components, dim, spacedim>::get_gradient(
+ const unsigned int q_point) const
+{
+ AssertIndexRange(q_point, gradients.size());
+ return gradients[q_point];
+}
+
DEAL_II_NAMESPACE_CLOSE
#endif
std::vector<double> function_values(unit_points.size());
std::vector<Tensor<1, dim>> function_gradients(unit_points.size());
- std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
- std::vector<double> function_values_1(unit_points.size());
- std::vector<Tensor<1, dim>> function_gradients_1(unit_points.size());
-
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
evaluator.evaluate(cell,
unit_points,
solution_values,
- function_values_1,
- function_gradients_1);
+ EvaluationFlags::values | EvaluationFlags::gradients);
deallog << "Cell with center " << cell->center(true) << std::endl;
- for (unsigned int i = 0; i < function_values_1.size(); ++i)
+ for (unsigned int i = 0; i < function_values.size(); ++i)
deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i])
- << ": " << function_values_1[i] << " error value "
- << function_values[i] - function_values_1[i] << " error grad "
- << (function_gradients_1[i] - function_gradients[i]).norm()
+ << ": " << evaluator.get_value(i) << " error value "
+ << function_values[i] - evaluator.get_value(i) << " error grad "
+ << (evaluator.get_gradient(i) - function_gradients[i]).norm()
<< std::endl;
deallog << std::endl;
}
std::vector<double> function_values(unit_points.size());
std::vector<Tensor<1, dim>> function_gradients(unit_points.size());
- std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
- std::vector<double> function_values_1(unit_points.size());
- std::vector<Tensor<1, dim>> function_gradients_1(unit_points.size());
-
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
evaluator.evaluate(cell,
unit_points,
solution_values,
- function_values_1,
- function_gradients_1);
+ EvaluationFlags::values | EvaluationFlags::gradients);
deallog << "Cell with center " << cell->center(true) << std::endl;
- for (unsigned int i = 0; i < function_values_1.size(); ++i)
+ for (unsigned int i = 0; i < function_values.size(); ++i)
deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i])
- << ": " << function_values_1[i] << " error value "
- << function_values[i] - function_values_1[i] << " error grad "
- << (function_gradients_1[i] - function_gradients[i]).norm()
+ << ": " << evaluator.get_value(i) << " error value "
+ << function_values[i] - evaluator.get_value(i) << " error grad "
+ << (evaluator.get_gradient(i) - function_gradients[i]).norm()
<< std::endl;
deallog << std::endl;
}
std::vector<Tensor<1, dim>> function_values(unit_points.size());
std::vector<Tensor<2, dim>> function_gradients(unit_points.size());
- std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
- std::vector<Tensor<1, dim>> function_values_1(unit_points.size());
- std::vector<Tensor<2, dim>> function_gradients_1(unit_points.size());
- FEValuesExtractors::Vector extractor(0);
+ FEValuesExtractors::Vector extractor(0);
for (const auto &cell : dof_handler.active_cell_iterators())
{
evaluator.evaluate(cell,
unit_points,
solution_values,
- function_values_1,
- function_gradients_1);
+ EvaluationFlags::values | EvaluationFlags::gradients);
deallog << "Cell with center " << cell->center(true) << std::endl;
- for (unsigned int i = 0; i < function_values_1.size(); ++i)
+ for (unsigned int i = 0; i < function_values.size(); ++i)
deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i])
- << ": " << function_values_1[i] << " error value "
- << (function_values[i] - function_values_1[i]).norm()
+ << ": " << evaluator.get_value(i) << " error value "
+ << (function_values[i] - evaluator.get_value(i)).norm()
<< " error grad "
- << (function_gradients_1[i] - function_gradients[i]).norm()
+ << (evaluator.get_gradient(i) - function_gradients[i]).norm()
<< std::endl;
deallog << std::endl;
}
std::vector<double> function_values_scalar(unit_points.size());
std::vector<Tensor<1, dim>> function_gradients_scalar(unit_points.size());
- std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
-
- std::vector<Tensor<1, dim>> function_values_1(unit_points.size());
- std::vector<Tensor<2, dim>> function_gradients_1(unit_points.size());
- std::vector<double> function_values_scalar_1(unit_points.size());
- std::vector<Tensor<1, dim>> function_gradients_scalar_1(unit_points.size());
-
FEValuesExtractors::Vector extractor(0);
FEValuesExtractors::Scalar extractor_s(dim);
evaluator.evaluate(cell,
unit_points,
solution_values,
- function_values_1,
- function_gradients_1);
+ EvaluationFlags::values | EvaluationFlags::gradients);
evaluator_scalar.evaluate(cell,
unit_points,
solution_values,
- function_values_scalar_1,
- function_gradients_scalar_1);
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
deallog << "Cell with center " << cell->center(true) << std::endl;
- for (unsigned int i = 0; i < function_values_1.size(); ++i)
+ for (unsigned int i = 0; i < function_values.size(); ++i)
deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i])
- << ": " << function_values_1[i] << " error value "
- << (function_values[i] - function_values_1[i]).norm()
+ << ": " << evaluator.get_value(i) << " error value "
+ << (function_values[i] - evaluator.get_value(i)).norm()
<< " error grad "
- << (function_gradients_1[i] - function_gradients[i]).norm()
+ << (evaluator.get_gradient(i) - function_gradients[i]).norm()
<< " error scalar value "
- << function_values_scalar[i] - function_values_scalar_1[i]
+ << function_values_scalar[i] - evaluator_scalar.get_value(i)
<< " error scalar grad "
<< (function_gradients_scalar[i] -
- function_gradients_scalar_1[i])
+ evaluator_scalar.get_gradient(i))
.norm()
<< std::endl;
deallog << std::endl;
std::vector<double> function_values(unit_points.size());
std::vector<Tensor<1, dim>> function_gradients(unit_points.size());
- std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
- std::vector<double> function_values_1(unit_points.size());
- std::vector<Tensor<1, dim>> function_gradients_1(unit_points.size());
-
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
evaluator.evaluate(cell,
unit_points,
solution_values,
- function_values_1,
- function_gradients_1);
+ EvaluationFlags::values | EvaluationFlags::gradients);
deallog << "Cell with center " << cell->center(true) << std::endl;
- for (unsigned int i = 0; i < function_values_1.size(); ++i)
+ for (unsigned int i = 0; i < function_values.size(); ++i)
deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i])
- << ": " << function_values_1[i] << " error value "
- << function_values[i] - function_values_1[i] << " error grad "
- << (function_gradients_1[i] - function_gradients[i]).norm()
+ << ": " << evaluator.get_value(i) << " error value "
+ << function_values[i] - evaluator.get_value(i) << " error grad "
+ << (evaluator.get_gradient(i) - function_gradients[i]).norm()
<< std::endl;
deallog << std::endl;
}
std::vector<double> function_values(unit_points.size());
std::vector<Tensor<1, dim>> function_gradients(unit_points.size());
- std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
- std::vector<double> function_values_1(unit_points.size());
- std::vector<Tensor<1, dim>> function_gradients_1(unit_points.size());
-
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
evaluator.evaluate(cell,
unit_points,
solution_values,
- function_values_1,
- function_gradients_1);
+ EvaluationFlags::values | EvaluationFlags::gradients);
deallog << "Cell with center " << cell->center(true) << std::endl;
- for (unsigned int i = 0; i < function_values_1.size(); ++i)
+ for (unsigned int i = 0; i < function_values.size(); ++i)
deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i])
- << ": " << function_values_1[i] << " error value "
- << function_values[i] - function_values_1[i] << " error grad "
- << (function_gradients_1[i] - function_gradients[i]).norm()
+ << ": " << evaluator.get_value(i) << " error value "
+ << function_values[i] - evaluator.get_value(i) << " error grad "
+ << (evaluator.get_gradient(i) - function_gradients[i]).norm()
<< std::endl;
deallog << std::endl;
}
std::vector<Tensor<1, dim>> function_values(unit_points.size());
std::vector<Tensor<2, dim>> function_gradients(unit_points.size());
- std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
- std::vector<Tensor<1, dim>> function_values_1(unit_points.size());
- std::vector<Tensor<2, dim>> function_gradients_1(unit_points.size());
- FEValuesExtractors::Vector extractor(0);
+ FEValuesExtractors::Vector extractor(0);
for (const auto &cell : dof_handler.active_cell_iterators())
{
evaluator.evaluate(cell,
unit_points,
solution_values,
- function_values_1,
- function_gradients_1);
+ EvaluationFlags::values | EvaluationFlags::gradients);
deallog << "Cell with center " << cell->center(true) << std::endl;
- for (unsigned int i = 0; i < function_values_1.size(); ++i)
+ for (unsigned int i = 0; i < function_values.size(); ++i)
deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i])
- << ": " << function_values_1[i] << " error value "
- << (function_values[i] - function_values_1[i]).norm()
+ << ": " << evaluator.get_value(i) << " error value "
+ << (function_values[i] - evaluator.get_value(i)).norm()
<< " error grad "
- << (function_gradients_1[i] - function_gradients[i]).norm()
+ << (evaluator.get_gradient(i) - function_gradients[i]).norm()
<< std::endl;
deallog << std::endl;
}
std::vector<double> function_values_scalar(unit_points.size());
std::vector<Tensor<1, dim>> function_gradients_scalar(unit_points.size());
- std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
-
- std::vector<Tensor<1, dim>> function_values_1(unit_points.size());
- std::vector<Tensor<2, dim>> function_gradients_1(unit_points.size());
- std::vector<double> function_values_scalar_1(unit_points.size());
- std::vector<Tensor<1, dim>> function_gradients_scalar_1(unit_points.size());
-
FEValuesExtractors::Vector extractor(0);
FEValuesExtractors::Scalar extractor_s(dim);
evaluator.evaluate(cell,
unit_points,
solution_values,
- function_values_1,
- function_gradients_1);
+ EvaluationFlags::values | EvaluationFlags::gradients);
evaluator_scalar.evaluate(cell,
unit_points,
solution_values,
- function_values_scalar_1,
- function_gradients_scalar_1);
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
deallog << "Cell with center " << cell->center(true) << std::endl;
- for (unsigned int i = 0; i < function_values_1.size(); ++i)
+ for (unsigned int i = 0; i < function_values.size(); ++i)
deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i])
- << ": " << function_values_1[i] << " error value "
- << (function_values[i] - function_values_1[i]).norm()
+ << ": " << evaluator.get_value(i) << " error value "
+ << (function_values[i] - evaluator.get_value(i)).norm()
<< " error grad "
- << (function_gradients_1[i] - function_gradients[i]).norm()
+ << (evaluator.get_gradient(i) - function_gradients[i]).norm()
<< " error scalar value "
- << function_values_scalar[i] - function_values_scalar_1[i]
+ << function_values_scalar[i] - evaluator_scalar.get_value(i)
<< " error scalar grad "
<< (function_gradients_scalar[i] -
- function_gradients_scalar_1[i])
+ evaluator_scalar.get_gradient(i))
.norm()
<< std::endl;
deallog << std::endl;
std::vector<std::vector<Tensor<1, dim>>> function_gradients(
unit_points.size(), std::vector<Tensor<1, dim>>(dim + 1));
- std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
- std::vector<Tensor<1, dim + 1>> function_values_1(unit_points.size());
- std::vector<Tensor<1, dim + 1, Tensor<1, dim>>> function_gradients_1(
- unit_points.size());
-
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
evaluator.evaluate(cell,
unit_points,
solution_values,
- function_values_1,
- function_gradients_1);
+ EvaluationFlags::values | EvaluationFlags::gradients);
deallog << "Cell with center " << cell->center(true) << std::endl;
- for (unsigned int i = 0; i < function_values_1.size(); ++i)
+ for (unsigned int i = 0; i < function_values.size(); ++i)
{
deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i])
- << ": " << function_values_1[i] << " error value ";
+ << ": " << evaluator.get_value(i) << " error value ";
double error = 0;
for (unsigned int d = 0; d < dim + 1; ++d)
- error += std::abs(function_values[i][d] - function_values_1[i][d]);
+ error +=
+ std::abs(function_values[i][d] - evaluator.get_value(i)[d]);
deallog << error << " error grad ";
error = 0;
for (unsigned int d = 0; d < dim + 1; ++d)
error +=
- (function_gradients[i][d] - function_gradients_1[i][d]).norm();
+ (function_gradients[i][d] - evaluator.get_gradient(i)[d]).norm();
deallog << error << std::endl;
}
deallog << std::endl;
projects / dealii.git / commitdiff