--- /dev/null
+Changed: The integrate() function of FEPointEvaluation now multiplies the submitted
+entities internally with a JxW value to represent a proper integration like in FEEvaluation.
+The old behavior (where the entities are only multiplied with the test functions and summed up)
+is available through the new function test_and_sum().
+<br>
+(Maximilian Bergbauer, 2023/05/29)
evaluate(const ArrayView<const ScalarNumber> & solution_values,
const EvaluationFlags::EvaluationFlags &evaluation_flags);
+ /**
+ * 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 multiplied be
+ * the Jacobian determinant times the quadrature weight (JxW).
+ *
+ * @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 in `solution_values` not touched
+ * by this class will be set to zero.
+ *
+ * @param[in] integration_flags Flags specifying which quantities should be
+ * integrated at the points.
+ *
+ */
+ void
+ integrate(const ArrayView<ScalarNumber> & solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags);
+
/**
* 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.
+ * factor (in contrast to the integrate function of this class). This allows
+ * the class to naturally embed point information (e.g. particles) into a
+ * finite element formulation.
*
* @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.
+ * selected by the present class, the entries in `solution_values` not touched
+ * by this class will be set to zero.
*
* @param[in] integration_flags Flags specifying which quantities should be
* integrated at the points.
*
*/
void
- integrate(const ArrayView<ScalarNumber> & solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags);
+ test_and_sum(const ArrayView<ScalarNumber> & solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags);
/**
* Return the value at quadrature point number @p point_index after a call to
/**
* Fast path of the integrate function.
*/
+ template <bool do_JxW>
void
integrate_fast(const ArrayView<ScalarNumber> & solution_values,
const EvaluationFlags::EvaluationFlags &integration_flags);
/**
* Slow path of the integrate function using FEValues.
*/
+ template <bool do_JxW>
void
integrate_slow(const ArrayView<ScalarNumber> & solution_values,
const EvaluationFlags::EvaluationFlags &integration_flags);
+ /**
+ * Implementation of the integrate/test_and_sum function.
+ */
+ template <bool do_JxW>
+ void
+ do_integrate(const ArrayView<ScalarNumber> & solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags);
/**
* Number of quadrature batches of the current cell/face.
, n_q_points_scalar(numbers::invalid_unsigned_int)
, mapping(&mapping)
, fe(&fe)
+ , JxW_ptr(nullptr)
, use_face_path(false)
, update_flags(update_flags)
, mapping_info_on_the_fly(
, n_q_points_scalar(numbers::invalid_unsigned_int)
, mapping(&mapping_info.get_mapping())
, fe(&fe)
+ , JxW_ptr(nullptr)
, use_face_path(false)
, update_flags(mapping_info.get_update_flags())
, mapping_info(&mapping_info)
template <int n_components_, int dim, int spacedim, typename Number>
+template <bool do_JxW>
inline void
FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_fast(
const ArrayView<ScalarNumber> & solution_values,
for (unsigned int v = 0;
v < stride && (stride == 1 || q + v < n_q_points_scalar);
++v)
- ETT::get_value(value, v, values[qb * stride + v]);
+ {
+ const unsigned int offset = qb * stride + v;
+ if (do_JxW)
+ values[offset] *= JxW_ptr[offset];
+ ETT::get_value(value, v, values[offset]);
+ }
}
if (integration_flags & EvaluationFlags::gradients)
{
++v)
{
const unsigned int offset = qb * stride + v;
+ if (do_JxW)
+ gradients[offset] *= JxW_ptr[offset];
ETT::get_gradient(
gradient,
v,
template <int n_components_, int dim, int spacedim, typename Number>
+template <bool do_JxW>
inline void
FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_slow(
const ArrayView<ScalarNumber> & solution_values,
v < n_lanes_user_interface && q + v < n_points;
++v)
solution_values[i] += fe_values->shape_value(i, q + v) *
- ETT::access(values[qb], v, d);
+ ETT::access(values[qb], v, d) *
+ (do_JxW ? fe_values->JxW(q + v) : 1.);
else if (nonzero_shape_function_component[i][d])
for (unsigned int qb = 0, q = 0; q < n_points;
++qb, q += n_lanes_user_interface)
++v)
solution_values[i] +=
fe_values->shape_value_component(i, q + v, d) *
- ETT::access(values[qb], v, d);
+ ETT::access(values[qb], v, d) *
+ (do_JxW ? fe_values->JxW(q + v) : 1.);
}
}
v < n_lanes_user_interface && q + v < n_points;
++v)
solution_values[i] += fe_values->shape_grad(i, q + v) *
- ETT::access(gradients[qb], v, d);
+ ETT::access(gradients[qb], v, d) *
+ (do_JxW ? fe_values->JxW(q + v) : 1.);
else if (nonzero_shape_function_component[i][d])
for (unsigned int qb = 0, q = 0; q < n_points;
++qb, q += n_lanes_user_interface)
++v)
solution_values[i] +=
fe_values->shape_grad_component(i, q + v, d) *
- ETT::access(gradients[qb], v, d);
+ ETT::access(gradients[qb], v, d) *
+ (do_JxW ? fe_values->JxW(q + v) : 1.);
}
}
}
template <int n_components_, int dim, int spacedim, typename Number>
+template <bool do_JxW>
void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
+FEPointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
const ArrayView<ScalarNumber> & solution_values,
const EvaluationFlags::EvaluationFlags &integration_flags)
{
return;
}
+ Assert(
+ !do_JxW || JxW_ptr != nullptr,
+ ExcMessage(
+ "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
+
AssertDimension(solution_values.size(), fe->dofs_per_cell);
if (fast_path)
- integrate_fast(solution_values, integration_flags);
+ integrate_fast<do_JxW>(solution_values, integration_flags);
else
- integrate_slow(solution_values, integration_flags);
+ integrate_slow<do_JxW>(solution_values, integration_flags);
+}
+
+
+
+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)
+{
+ do_integrate<true>(solution_values, integration_flags);
+}
+
+
+
+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)
+{
+ do_integrate<false>(solution_values, integration_flags);
}
evaluator.submit_value(
values[q + cell_data.reference_point_ptrs[cell]], q);
- evaluator.integrate(solution_values, EvaluationFlags::values);
+ evaluator.test_and_sum(solution_values, EvaluationFlags::values);
// resolve constraints and scatter
internal::VectorDistributorLocalToGlobal<Number, VectorizedArrayType>
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << factor_float * i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator_scalar.submit_gradient(evaluator_scalar.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
deallog << std::endl;
- evaluator_scalar.integrate(solution_values,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
+ evaluator_scalar.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator_scalar.submit_gradient(evaluator_scalar.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
deallog << std::endl;
- evaluator_scalar.integrate(solution_values,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
+ evaluator_scalar.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << i << ' ';
evaluator2.submit_gradient(evaluator2.get_gradient(i), i);
}
- evaluator1.integrate(solution_values,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
- evaluator2.integrate(solution_values2,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
+ evaluator1.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
+ evaluator2.test_and_sum(solution_values2,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (unsigned int i = 0; i < solution_values.size(); ++i)
deallog << solution_values[i] - solution_values2[i] << ' ';
evaluator2.submit_gradient(evaluator2.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
- evaluator2.integrate(solution_values2,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
+ evaluator2.test_and_sum(solution_values2,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (unsigned int i = 0; i < solution_values.size(); ++i)
deallog << solution_values[i] - solution_values2[i] << ' ';
evaluator2.submit_gradient(evaluator2.get_gradient(i), i);
}
- evaluator1.integrate(solution_values,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
- evaluator2.integrate(solution_values2,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
+ evaluator1.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
+ evaluator2.test_and_sum(solution_values2,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (unsigned int i = 0; i < solution_values.size(); ++i)
deallog << solution_values[i] - solution_values2[i] << ' ';
evaluator_move.submit_gradient(evaluator_move.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << factor_float * i << ' ';
- evaluator_move.integrate(solution_values,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
+ evaluator_move.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
deallog << std::endl;
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values,
- EvaluationFlags::values | EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (const auto i : solution_values)
deallog << factor_float * i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values_out,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values_out,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (double i : solution_values_out)
deallog << i << ' ';
evaluator.submit_gradient(evaluator.get_gradient(i), i);
}
- evaluator.integrate(solution_values_out,
- EvaluationFlags::values |
- EvaluationFlags::gradients);
+ evaluator.test_and_sum(solution_values_out,
+ EvaluationFlags::values |
+ EvaluationFlags::gradients);
for (double i : solution_values_out)
deallog << i << ' ';
for (const auto q : fe_point.quadrature_point_indices())
{
- fe_point.submit_value(fe_point.JxW(q) * fe_point.get_value(q), q);
- fe_point.submit_gradient(fe_point.JxW(q) * fe_point.get_gradient(q),
- q);
+ fe_point.submit_value(fe_point.get_value(q), q);
+ fe_point.submit_gradient(fe_point.get_gradient(q), q);
}
fe_point.integrate(solution_values_out,
for (const auto q : fe_point_faces_m.quadrature_point_indices())
{
- fe_point_faces_m.submit_value(fe_point_faces_m.JxW(q) *
- (fe_point_faces_m.get_value(q) -
- fe_point_faces_p.get_value(q)),
+ fe_point_faces_m.submit_value(fe_point_faces_m.get_value(q) -
+ fe_point_faces_p.get_value(q),
q);
fe_point_faces_m.submit_gradient(
- fe_point_faces_m.JxW(q) * (fe_point_faces_m.get_gradient(q) -
- fe_point_faces_p.get_gradient(q)),
+ fe_point_faces_m.get_gradient(q) -
+ fe_point_faces_p.get_gradient(q),
q);
}
for (unsigned int q = 0; q < unit_points.size(); ++q)
phi_force.submit_value(force_JxW[q], q);
- phi_force.integrate(buffer, EvaluationFlags::values);
+ phi_force.test_and_sum(buffer, EvaluationFlags::values);
constraints.distribute_local_to_global(buffer,
local_dof_indices,
// integrate_scatter force
{
- phi_force.integrate(buffer_dim, EvaluationFlags::values);
+ phi_force.test_and_sum(buffer_dim, EvaluationFlags::values);
constraints.distribute_local_to_global(buffer_dim,
local_dof_indices_dim,