* @param[in] integration_flags Flags specifying which quantities should be
* integrated at the points.
*
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
*/
template <std::size_t stride_view>
void
integrate(const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags);
-
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
/**
* This function multiplies the quantities passed in by previous
* @param[in] integration_flags Flags specifying which quantities should be
* integrated at the points.
*
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
*/
void
integrate(const ArrayView<ScalarNumber> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags)
- {
- integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
- solution_values.size()),
- integration_flags);
- }
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
/**
* This function multiplies the quantities passed in by previous
* @param[in] integration_flags Flags specifying which quantities should be
* integrated at the points.
*
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
*/
template <std::size_t stride_view>
void
test_and_sum(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags);
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
/**
* This function multiplies the quantities passed in by previous
* @param[in] integration_flags Flags specifying which quantities should be
* integrated at the points.
*
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
*/
void
test_and_sum(const ArrayView<ScalarNumber> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags)
- {
- test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
- solution_values.size()),
- integration_flags);
- }
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
/**
* Return the value at quadrature point number @p point_index after a call to
finish_integrate_fast(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
const EvaluationFlags::EvaluationFlags &integration_flags,
- vectorized_value_type *solution_values_vectorized_linear);
+ vectorized_value_type *solution_values_vectorized_linear,
+ const bool sum_into_values);
/**
* Fast path of the integrate function.
void
integrate_fast(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags);
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
/**
* Slow path of the integrate function using FEValues.
void
integrate_slow(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags);
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
/**
* Implementation of the integrate/test_and_sum function.
void
do_integrate(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags);
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
/**
* Number of quadrature batches of the current cell/face.
FEPointEvaluation<n_components_, dim, spacedim, Number>::finish_integrate_fast(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
const EvaluationFlags::EvaluationFlags &integration_flags,
- vectorized_value_type *solution_values_vectorized_linear)
+ vectorized_value_type *solution_values_vectorized_linear,
+ const bool sum_into_values)
{
- if (fe->n_components() > n_components)
+ if (!sum_into_values && fe->n_components() > n_components)
for (unsigned int i = 0; i < solution_values.size(); ++i)
solution_values[i] = 0;
if (is_linear || renumber.empty())
{
for (unsigned int i = 0; i < dofs_per_comp; ++i)
- solution_values[i + offset] = output[i];
+ if (sum_into_values)
+ solution_values[i + offset] += output[i];
+ else
+ solution_values[i + offset] = output[i];
}
else
{
const unsigned int *renumber_ptr = renumber.data() + offset;
for (unsigned int i = 0; i < dofs_per_comp; ++i)
- solution_values[renumber_ptr[i]] = output[i];
+ if (sum_into_values)
+ solution_values[renumber_ptr[i]] += output[i];
+ else
+ solution_values[renumber_ptr[i]] = output[i];
}
}
else
if (is_linear || renumber.empty())
{
for (unsigned int i = 0; i < dofs_per_comp; ++i)
- solution_values[i + offset] =
- ETT::sum_value(comp,
- is_linear ?
- *(solution_values_vectorized_linear + i) :
- solution_renumbered_vectorized[i]);
+ if (sum_into_values)
+ solution_values[i + offset] +=
+ ETT::sum_value(comp,
+ is_linear ?
+ *(solution_values_vectorized_linear + i) :
+ solution_renumbered_vectorized[i]);
+ else
+ solution_values[i + offset] =
+ ETT::sum_value(comp,
+ is_linear ?
+ *(solution_values_vectorized_linear + i) :
+ solution_renumbered_vectorized[i]);
}
else
{
const unsigned int *renumber_ptr = renumber.data() + offset;
for (unsigned int i = 0; i < dofs_per_comp; ++i)
- solution_values[renumber_ptr[i]] =
- ETT::sum_value(comp, solution_renumbered_vectorized[i]);
+ if (sum_into_values)
+ solution_values[renumber_ptr[i]] +=
+ ETT::sum_value(comp, solution_renumbered_vectorized[i]);
+ else
+ solution_values[renumber_ptr[i]] =
+ ETT::sum_value(comp, solution_renumbered_vectorized[i]);
}
}
}
inline void
FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_fast(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags)
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
{
// zero out lanes of incomplete last quadrature point batch
if constexpr (stride == 1)
finish_integrate_fast<is_face_path, is_linear>(
solution_values,
integration_flags,
- solution_values_vectorized_linear.data());
+ solution_values_vectorized_linear.data(),
+ sum_into_values);
}
inline void
FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_slow(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags)
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
{
// slow path with FEValues
Assert(fe_values.get() != nullptr,
ExcMessage(
"Not initialized. Please call FEPointEvaluation::reinit()!"));
- for (unsigned int i = 0; i < solution_values.size(); ++i)
- solution_values[i] = 0;
+ if (!sum_into_values)
+ for (unsigned int i = 0; i < solution_values.size(); ++i)
+ solution_values[i] = 0;
const std::size_t n_points = fe_values->get_quadrature().size();
void
FEPointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags)
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
{
if (!is_reinitialized)
reinit();
(integration_flags &
EvaluationFlags::gradients))) // no integration flags
{
- for (unsigned int i = 0; i < solution_values.size(); ++i)
- solution_values[i] = 0;
+ if (!sum_into_values)
+ for (unsigned int i = 0; i < solution_values.size(); ++i)
+ solution_values[i] = 0;
return;
}
{
if (use_linear_path)
integrate_fast<do_JxW, true, true>(solution_values,
- integration_flags);
+ integration_flags,
+ sum_into_values);
else
integrate_fast<do_JxW, true, false>(solution_values,
- integration_flags);
+ integration_flags,
+ sum_into_values);
}
else
{
if (use_linear_path)
integrate_fast<do_JxW, false, true>(solution_values,
- integration_flags);
+ integration_flags,
+ sum_into_values);
else
integrate_fast<do_JxW, false, false>(solution_values,
- integration_flags);
+ integration_flags,
+ sum_into_values);
}
}
else
- integrate_slow<do_JxW>(solution_values, integration_flags);
+ integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
}
void
FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags)
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
+{
+ do_integrate<true>(solution_values, integration_flags, sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
+ const ArrayView<ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
{
- do_integrate<true>(solution_values, integration_flags);
+ integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
+ solution_values.size()),
+ integration_flags,
+ sum_into_values);
}
void
FEPointEvaluation<n_components_, dim, spacedim, Number>::test_and_sum(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags)
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
+{
+ do_integrate<false>(solution_values, integration_flags, sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::test_and_sum(
+ const ArrayView<ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
{
- do_integrate<false>(solution_values, integration_flags);
+ test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
+ solution_values.size()),
+ integration_flags,
+ sum_into_values);
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check FEPointEvaluation integrate_add()
+
+#include <deal.II/base/function_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim, typename Number = double>
+void
+test(const unsigned int degree)
+{
+ using namespace dealii;
+ Triangulation<dim> tria;
+
+ if (dim > 1)
+ GridGenerator::hyper_shell(tria, Point<dim>(), 0.5, 1, 6);
+ else
+ GridGenerator::subdivided_hyper_cube(tria, 2, 0, 1);
+
+ MappingQ<dim> mapping(degree);
+ deallog << "Mapping of degree " << degree << std::endl;
+
+ std::vector<Point<dim>> unit_points;
+ for (unsigned int i = 0; i < Utilities::pow(2, dim); ++i)
+ {
+ Point<dim> p;
+ for (unsigned int d = 0; d < dim; ++d)
+ p[d] = static_cast<Number>(i) / 17. + 0.015625 * d;
+ unit_points.push_back(p);
+ }
+
+ QGauss<dim> q_gauss(2);
+
+ Quadrature<dim> quad(unit_points, q_gauss.get_weights());
+
+ FE_Q<dim> fe(degree);
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+ Vector<Number> vector(dof_handler.n_dofs());
+
+ NonMatching::MappingInfo mapping_info(mapping,
+ update_values | update_gradients |
+ update_JxW_values);
+
+ std::vector<Quadrature<dim>> quad_vec_cell;
+ quad_vec_cell.reserve(tria.n_cells());
+
+ for (const auto &cell : tria.active_cell_iterators())
+ quad_vec_cell.push_back(quad);
+
+ mapping_info.reinit_cells(tria.active_cell_iterators(), quad_vec_cell);
+
+ FEPointEvaluation<1, dim, dim, Number> evaluator(mapping_info, fe);
+
+ Tensor<1, dim, Number> exponents;
+ exponents[0] = 1.;
+ VectorTools::interpolate(mapping,
+ dof_handler,
+ Functions::Monomial<dim, Number>(exponents),
+ vector);
+
+ std::vector<Number> solution_values(fe.dofs_per_cell);
+
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ cell->get_dof_values(vector,
+ solution_values.begin(),
+ solution_values.end());
+
+ evaluator.reinit(cell->active_cell_index());
+ evaluator.evaluate(solution_values,
+ EvaluationFlags::values | EvaluationFlags::gradients);
+
+ for (const unsigned int i : evaluator.quadrature_point_indices())
+ {
+ evaluator.submit_value(evaluator.get_value(i), i);
+ evaluator.submit_gradient(evaluator.get_gradient(i), i);
+ }
+
+ evaluator.integrate(solution_values,
+ EvaluationFlags::values | EvaluationFlags::gradients);
+
+ deallog << "Cell with center " << cell->center(true) << std::endl;
+ for (const auto i : solution_values)
+ deallog << i << ' ';
+ deallog << std::endl;
+
+ evaluator.integrate(solution_values,
+ EvaluationFlags::values | EvaluationFlags::gradients,
+ true);
+
+ for (const auto i : solution_values)
+ deallog << i << ' ';
+ deallog << std::endl;
+ }
+}
+
+
+
+int
+main()
+{
+ initlog();
+ deallog << std::setprecision(10);
+
+ test<2>(1);
+ test<2>(2);
+}
--- /dev/null
+
+DEAL::Mapping of degree 1
+DEAL::Cell with center 0.6495190528 0.3750000000
+DEAL::1.094128055 0.06405871628 -0.7129529194 -0.07186179469
+DEAL::2.188256110 0.1281174326 -1.425905839 -0.1437235894
+DEAL::Cell with center 4.592425497e-17 0.7500000000
+DEAL::0.5256716676 -0.3778726922 0.05767485942 -0.04324345223
+DEAL::1.051343335 -0.7557453844 0.1153497188 -0.08648690447
+DEAL::Cell with center -0.6495190528 0.3750000000
+DEAL::-0.5684563874 -0.4419314085 0.7706277789 0.02861834246
+DEAL::-1.136912775 -0.8838628170 1.541255558 0.05723668491
+DEAL::Cell with center -0.6495190528 -0.3750000000
+DEAL::-1.094128055 -0.06405871628 0.7129529194 0.07186179469
+DEAL::-2.188256110 -0.1281174326 1.425905839 0.1437235894
+DEAL::Cell with center -1.377727649e-16 -0.7500000000
+DEAL::-0.5256716676 0.3778726922 -0.05767485942 0.04324345223
+DEAL::-1.051343335 0.7557453844 -0.1153497188 0.08648690447
+DEAL::Cell with center 0.6495190528 -0.3750000000
+DEAL::0.5684563874 0.4419314085 -0.7706277789 -0.02861834246
+DEAL::1.136912775 0.8838628170 -1.541255558 -0.05723668491
+DEAL::Mapping of degree 2
+DEAL::Cell with center 0.6495190528 0.3750000000
+DEAL::2.482631076 -0.1564725804 0.4910617067 -0.02558423728 -2.474880528 0.1748480233 0.7176813542 0.1215956751 -0.8420572937
+DEAL::4.965262152 -0.3129451607 0.9821234133 -0.05116847457 -4.949761056 0.3496960465 1.435362708 0.2431913503 -1.684114587
+DEAL::Cell with center 4.592425497e-17 0.7500000000
+DEAL::1.903653500 0.1627020517 0.1365654459 -0.02631613009 -0.6782922627 0.1422251499 -0.8006550044 0.1439386250 -0.7778032374
+DEAL::3.807306999 0.3254041034 0.2731308918 -0.05263226019 -1.356584525 0.2844502999 -1.601310009 0.2878772501 -1.555606475
+DEAL::Cell with center -0.6495190528 0.3750000000
+DEAL::-0.5789775763 0.3191746321 -0.3544962608 -0.0007318928105 1.796588265 -0.03262287333 -1.518336359 0.02234294989 0.06425405628
+DEAL::-1.157955153 0.6383492641 -0.7089925215 -0.001463785621 3.593176530 -0.06524574667 -3.036672717 0.04468589979 0.1285081126
+DEAL::Cell with center -0.6495190528 -0.3750000000
+DEAL::-2.482631076 0.1564725804 -0.4910617067 0.02558423728 2.474880528 -0.1748480233 -0.7176813542 -0.1215956751 0.8420572937
+DEAL::-4.965262152 0.3129451607 -0.9821234133 0.05116847457 4.949761056 -0.3496960465 -1.435362708 -0.2431913503 1.684114587
+DEAL::Cell with center -1.377727649e-16 -0.7500000000
+DEAL::-1.903653500 -0.1627020517 -0.1365654459 0.02631613009 0.6782922627 -0.1422251499 0.8006550044 -0.1439386250 0.7778032374
+DEAL::-3.807306999 -0.3254041034 -0.2731308918 0.05263226019 1.356584525 -0.2844502999 1.601310009 -0.2878772501 1.555606475
+DEAL::Cell with center 0.6495190528 -0.3750000000
+DEAL::0.5789775763 -0.3191746321 0.3544962608 0.0007318928105 -1.796588265 0.03262287333 1.518336359 -0.02234294989 -0.06425405628
+DEAL::1.157955153 -0.6383492641 0.7089925215 0.001463785621 -3.593176530 0.06524574667 3.036672717 -0.04468589979 -0.1285081126
projects / dealii.git / commitdiff