DerivativeForm &
operator=(const Tensor<1, dim, Number> &);
+ /**
+ * Number conversion operator.
+ */
+ template <typename OtherNumber>
+ DerivativeForm &
+ operator=(const DerivativeForm<order, dim, spacedim, OtherNumber> &df);
+
/**
* Converts a DerivativeForm <order, dim, dim, Number> to Tensor<order+1, dim,
* Number>. In particular, if order == 1 and the derivative is the Jacobian of
+template <int order, int dim, int spacedim, typename Number>
+template <typename OtherNumber>
+inline DerivativeForm<order, dim, spacedim, Number> &
+DerivativeForm<order, dim, spacedim, Number>::operator=(
+ const DerivativeForm<order, dim, spacedim, OtherNumber> &df)
+{
+ for (unsigned int j = 0; j < spacedim; ++j)
+ (*this)[j] = df[j];
+ return *this;
+}
+
+
+
template <int order, int dim, int spacedim, typename Number>
inline Tensor<order, dim, Number> &
DerivativeForm<order, dim, spacedim, Number>::operator[](const unsigned int i)
template <typename T>
struct VectorizedArrayTrait
{
- using value_type = T;
+ /**
+ * Define scalar value type.
+ */
+ using value_type = T;
+
+ /**
+ * Define width of template type.
+ */
static constexpr std::size_t width = 1;
- static T &
- get(T &value, unsigned int c)
+ /**
+ * Define vectorized value type for internal vectorization.
+ */
+ using vectorized_value_type = VectorizedArray<T>;
+
+ /**
+ * Define a stride which defines how often the template type T fits into the
+ * vectorized_value_type. This is useful to write vectorized templated code
+ * where the internal computation is vectorized and the user interface is
+ * optionally scalar or also vectorized.
+ */
+ static constexpr std::size_t stride = vectorized_value_type::size();
+
+ /**
+ * Get a reference to scalar value (on lane 0).
+ */
+ static value_type &
+ get(value_type &value, unsigned int c)
{
- AssertDimension(c, 0);
+ AssertIndexRange(c, 1);
(void)c;
return value;
}
- static const T &
- get(const T &value, unsigned int c)
+ /**
+ * Get a read-only reference to scalar value (on lane 0).
+ */
+ static const value_type &
+ get(const value_type &value, unsigned int c)
{
- AssertDimension(c, 0);
+ AssertIndexRange(c, 1);
(void)c;
return value;
}
+
+ /**
+ * Get a reference to scalar value on lane c from a vectorized values field.
+ */
+ static value_type &
+ get_from_vectorized(vectorized_value_type &values, unsigned int c)
+ {
+ AssertIndexRange(c, stride);
+
+ return values[c];
+ }
+
+ /**
+ * Get a read-only reference to scalar value on lane c from a vectorized
+ * values field.
+ */
+ static const value_type &
+ get_from_vectorized(const vectorized_value_type &values, unsigned int c)
+ {
+ AssertIndexRange(c, stride);
+
+ return values[c];
+ }
};
template <typename T, std::size_t width_>
struct VectorizedArrayTrait<VectorizedArray<T, width_>>
{
- using value_type = T;
+ /**
+ * Define scalar value type.
+ */
+ using value_type = T;
+
+ /**
+ * Define width of template type.
+ */
static constexpr std::size_t width = width_;
- static T &
- get(VectorizedArray<T, width_> &values, unsigned int c)
+ /**
+ * Define vectorized value type for internal vectorization.
+ */
+ using vectorized_value_type = VectorizedArray<T, width_>;
+
+ /**
+ * Define a stride which defines how often the template type
+ * VectorizedArray<T, width_> fits into the vectorized value type. This is
+ * useful to write vectorized templated code where the internal computation
+ * is vectorized and the user interface is optionally scalar or also
+ * vectorized.
+ */
+ static constexpr std::size_t stride = 1;
+
+ /**
+ * Get a reference to scalar value on lane c.
+ */
+ static value_type &
+ get(vectorized_value_type &values, unsigned int c)
{
AssertIndexRange(c, width_);
return values[c];
}
- static const T &
- get(const VectorizedArray<T, width_> &values, unsigned int c)
+ /**
+ * Get a read-only reference to scalar value on lane c.
+ */
+ static const value_type &
+ get(const vectorized_value_type &values, unsigned int c)
{
AssertIndexRange(c, width_);
return values[c];
}
+
+ /**
+ * Get a reference to vectorized values from a vectorized values field.
+ */
+ static vectorized_value_type &
+ get_from_vectorized(vectorized_value_type &values, unsigned int c)
+ {
+ (void)c;
+ AssertIndexRange(c, stride);
+
+ return values;
+ }
+
+ /**
+ * Get a read-only reference to vectorized values from a vectorized values
+ * field.
+ */
+ static const vectorized_value_type &
+ get_from_vectorized(const vectorized_value_type &values, unsigned int c)
+ {
+ (void)c;
+ AssertIndexRange(c, stride);
+
+ return values;
+ }
};
} // namespace internal
{
template <int dim>
class FEImmersedSurfaceValues;
- template <int dim, int spacedim>
- class MappingInfo;
+ namespace internal
+ {
+ template <int dim, int spacedim>
+ class ComputeMappingDataHelper;
+ }
} // namespace NonMatching
friend class FEFaceValues<dim, spacedim>;
friend class FESubfaceValues<dim, spacedim>;
friend class NonMatching::FEImmersedSurfaceValues<dim>;
- friend class NonMatching::MappingInfo<dim, spacedim>;
+ friend class NonMatching::internal::ComputeMappingDataHelper<dim, spacedim>;
};
{
namespace FEPointEvaluation
{
- template <typename Number>
- struct VectorizedPointsTypeTraits
- {
- using value_type = VectorizedArray<Number>;
- static constexpr std::size_t stride = value_type::size();
- };
-
- template <typename Number, std::size_t width_>
- struct VectorizedPointsTypeTraits<VectorizedArray<Number, width_>>
- {
- using value_type = VectorizedArray<Number, width_>;
- static constexpr std::size_t stride = 1;
- };
-
+ DeclException1(
+ ExcFEPointEvaluationAccessToUninitializedMappingField,
+ std::string,
+ << "You are requesting information from an FEPointEvaluation "
+ << "object for which this kind of information has not been computed. "
+ << "What information these objects compute is determined by the update_* "
+ << "flags you pass to MappingInfo() in the Constructor. Here, "
+ << "the operation you are attempting requires the <" << arg1
+ << "> flag to be set, but it was apparently not specified "
+ << "upon initialization.");
/**
* Struct to distinguish between the value and gradient types of different
using ScalarNumber =
typename internal::VectorizedArrayTrait<Number>::value_type;
using VectorizedArrayType =
- typename internal::FEPointEvaluation::VectorizedPointsTypeTraits<
- Number>::value_type;
+ typename dealii::internal::VectorizedArrayTrait<
+ Number>::vectorized_value_type;
using value_type = Tensor<1, n_components, Number>;
using scalar_value_type = Tensor<1, n_components, ScalarNumber>;
using vectorized_value_type =
Tensor<1, dim, Tensor<1, n_components, VectorizedArrayType>>;
static void
- read_value(const Number vector_entry,
+ read_value(const ScalarNumber vector_entry,
const unsigned int component,
- value_type & result)
+ scalar_value_type &result)
{
AssertIndexRange(component, n_components);
result[component] = vector_entry;
}
static void
- write_value(Number & vector_entry,
- const unsigned int component,
- const value_type & result)
+ write_value(VectorizedArrayType & vector_entry,
+ const unsigned int component,
+ const vectorized_value_type &result)
{
AssertIndexRange(component, n_components);
vector_entry = result[component];
static void
set_gradient(const interface_vectorized_gradient_type &value,
const unsigned int vector_lane,
- scalar_gradient_type & result)
+ gradient_type & result)
{
for (unsigned int i = 0; i < n_components; ++i)
for (unsigned int d = 0; d < dim; ++d)
- result[i][d] = value[d][i][vector_lane];
- }
-
- static void
- set_gradient(const interface_vectorized_gradient_type &value,
- const unsigned int,
- vectorized_gradient_type &result)
- {
- for (unsigned int i = 0; i < n_components; ++i)
- for (unsigned int d = 0; d < dim; ++d)
- result[i][d] = value[d][i];
+ result[i][d] =
+ internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+ value[d][i], vector_lane);
}
static void
get_gradient(interface_vectorized_gradient_type &value,
const unsigned int vector_lane,
- const scalar_gradient_type & result)
+ const gradient_type & result)
{
for (unsigned int i = 0; i < n_components; ++i)
for (unsigned int d = 0; d < dim; ++d)
- value[d][i][vector_lane] = result[i][d];
- }
-
- static void
- get_gradient(interface_vectorized_gradient_type &value,
- const unsigned int,
- const vectorized_gradient_type &result)
- {
- for (unsigned int i = 0; i < n_components; ++i)
- for (unsigned int d = 0; d < dim; ++d)
- value[d][i] = result[i][d];
+ internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+ value[d][i], vector_lane) = result[i][d];
}
static void
using ScalarNumber =
typename internal::VectorizedArrayTrait<Number>::value_type;
using VectorizedArrayType =
- typename internal::FEPointEvaluation::VectorizedPointsTypeTraits<
- Number>::value_type;
+ typename dealii::internal::VectorizedArrayTrait<
+ Number>::vectorized_value_type;
using value_type = Number;
using scalar_value_type = ScalarNumber;
using vectorized_value_type = VectorizedArrayType;
using vectorized_gradient_type = Tensor<1, dim, VectorizedArrayType>;
static void
- read_value(const Number vector_entry,
+ read_value(const ScalarNumber vector_entry,
const unsigned int,
- value_type &result)
+ scalar_value_type &result)
{
result = vector_entry;
}
static void
- write_value(Number &vector_entry,
+ write_value(VectorizedArrayType &vector_entry,
const unsigned int,
- const value_type &result)
+ const vectorized_value_type &result)
{
vector_entry = result;
}
using ScalarNumber =
typename internal::VectorizedArrayTrait<Number>::value_type;
using VectorizedArrayType =
- typename internal::FEPointEvaluation::VectorizedPointsTypeTraits<
- Number>::value_type;
+ typename dealii::internal::VectorizedArrayTrait<
+ Number>::vectorized_value_type;
using value_type = Tensor<1, dim, Number>;
using scalar_value_type = Tensor<1, dim, ScalarNumber>;
using vectorized_value_type = Tensor<1, dim, VectorizedArrayType>;
Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>>;
static void
- read_value(const Number vector_entry,
+ read_value(const ScalarNumber vector_entry,
const unsigned int component,
- value_type & result)
+ scalar_value_type &result)
{
result[component] = vector_entry;
}
static void
- write_value(Number & vector_entry,
- const unsigned int component,
- const value_type & result)
+ write_value(VectorizedArrayType & vector_entry,
+ const unsigned int component,
+ const vectorized_value_type &result)
{
vector_entry = result[component];
}
static void
set_gradient(const interface_vectorized_gradient_type &value,
const unsigned int vector_lane,
- scalar_gradient_type & result)
+ gradient_type & result)
{
for (unsigned int i = 0; i < dim; ++i)
for (unsigned int d = 0; d < dim; ++d)
- result[i][d] = value[d][i][vector_lane];
- }
-
- static void
- set_gradient(const interface_vectorized_gradient_type &value,
- const unsigned int,
- vectorized_gradient_type &result)
- {
- for (unsigned int i = 0; i < dim; ++i)
- for (unsigned int d = 0; d < dim; ++d)
- result[i][d] = value[d][i];
+ result[i][d] =
+ internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+ value[d][i], vector_lane);
}
static void
get_gradient(interface_vectorized_gradient_type &value,
const unsigned int vector_lane,
- const scalar_gradient_type & result)
- {
- for (unsigned int i = 0; i < dim; ++i)
- for (unsigned int d = 0; d < dim; ++d)
- value[d][i][vector_lane] = result[i][d];
- }
-
- static void
- get_gradient(interface_vectorized_gradient_type &value,
- const unsigned int,
- const vectorized_gradient_type &result)
+ const gradient_type & result)
{
for (unsigned int i = 0; i < dim; ++i)
for (unsigned int d = 0; d < dim; ++d)
- value[d][i] = result[i][d];
+ internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+ value[d][i], vector_lane) = result[i][d];
}
static void
using ScalarNumber =
typename internal::VectorizedArrayTrait<Number>::value_type;
using VectorizedArrayType =
- typename internal::FEPointEvaluation::VectorizedPointsTypeTraits<
- Number>::value_type;
+ typename dealii::internal::VectorizedArrayTrait<
+ Number>::vectorized_value_type;
using value_type = Number;
using scalar_value_type = ScalarNumber;
using vectorized_value_type = VectorizedArrayType;
using vectorized_gradient_type = Tensor<1, 1, VectorizedArrayType>;
static void
- read_value(const Number vector_entry,
+ read_value(const ScalarNumber vector_entry,
const unsigned int,
- value_type &result)
+ scalar_value_type &result)
{
result = vector_entry;
}
static void
- write_value(Number &vector_entry,
+ write_value(VectorizedArrayType &vector_entry,
const unsigned int,
- const value_type &result)
+ const vectorized_value_type &result)
{
vector_entry = result;
}
using ScalarNumber =
typename internal::VectorizedArrayTrait<Number>::value_type;
- using VectorizedArrayType =
- typename internal::FEPointEvaluation::VectorizedPointsTypeTraits<
- Number>::value_type;
+ using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
+ Number>::vectorized_value_type;
using ETT = typename internal::FEPointEvaluation::
EvaluatorTypeTraits<dim, n_components, Number>;
using value_type = typename ETT::value_type;
* objects, this parameter allows to select a range of `n_components`
* components starting from this parameter.
*/
- FEPointEvaluation(NonMatching::MappingInfo<dim, spacedim> &mapping_info,
- const FiniteElement<dim> & fe,
- const unsigned int first_selected_component = 0);
+ FEPointEvaluation(
+ NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+ const FiniteElement<dim> & fe,
+ const unsigned int first_selected_component = 0);
/**
* Copy constructor.
quadrature_point_indices() const;
private:
+ static constexpr std::size_t n_lanes_user_interface =
+ internal::VectorizedArrayTrait<Number>::width;
+ static constexpr std::size_t n_lanes_internal =
+ internal::VectorizedArrayTrait<VectorizedArrayType>::width;
+ static constexpr std::size_t stride =
+ internal::VectorizedArrayTrait<Number>::stride;
+
/**
* Common setup function for both constructors. Does the setup for both fast
* and slow path.
integrate_slow(const ArrayView<ScalarNumber> & solution_values,
const EvaluationFlags::EvaluationFlags &integration_flags);
+
/**
- * Number of quadrature points of the current cell/face.
+ * Number of quadrature batches of the current cell/face.
*/
const unsigned int n_q_points;
/**
- * Number of active quadrature points of the last quadrature point batch.
+ * Number of quadrature points of the current cell/face.
*/
- const unsigned int n_filled_lanes_last_batch;
+ const unsigned int n_q_points_scalar;
/**
* Pointer to the Mapping object passed to the constructor.
*/
std::vector<gradient_type> gradients;
+ /**
+ * Pointer to first unit point batch of current cell/face from MappingInfo,
+ * set internally during do_reinit().
+ */
+ const Point<dim, VectorizedArrayType> *unit_point_ptr;
+
+ /**
+ * Pointer to real point of first quadrature point of current cell/face from
+ * MappingInfo, set internally during do_reinit().
+ */
+ const Point<spacedim, Number> *real_point_ptr;
+
+ /**
+ * Pointer to Jacobian of first quadrature point of current cell/face from
+ * MappingInfo, set internally during do_reinit().
+ */
+ const DerivativeForm<1, dim, spacedim, Number> *jacobian_ptr;
+
+ /**
+ * Pointer to inverse Jacobian of first quadrature point of current cell/face
+ * from MappingInfo, set internally during do_reinit().
+ */
+ const DerivativeForm<1, spacedim, dim, Number> *inverse_jacobian_ptr;
+
+ /**
+ * Pointer to normal vector of first quadrature point of current cell/face
+ * from MappingInfo, set internally during do_reinit().
+ */
+ const Tensor<1, spacedim, Number> *normal_ptr;
+
+ /**
+ * Pointer to Jacobian determinant times quadrature weight of first quadrature
+ * point of current cell/face from MappingInfo, set internally during
+ * do_reinit().
+ */
+ const Number *JxW_ptr;
+
/**
* Number of unknowns per component, i.e., number of unique basis functions,
* for the chosen FiniteElement (or base element).
/**
* Pointer to mapping info on the fly computed during reinit.
*/
- std::unique_ptr<NonMatching::MappingInfo<dim, spacedim>>
+ std::unique_ptr<NonMatching::MappingInfo<dim, spacedim, Number>>
mapping_info_on_the_fly;
/**
* Pointer to currently used mapping info (either on the fly or external
* precomputed).
*/
- SmartPointer<NonMatching::MappingInfo<dim, spacedim>> mapping_info;
+ SmartPointer<NonMatching::MappingInfo<dim, spacedim, Number>> mapping_info;
/**
* The current cell index to access mapping data from mapping info.
const UpdateFlags update_flags,
const unsigned int first_selected_component)
: n_q_points(numbers::invalid_unsigned_int)
- , n_filled_lanes_last_batch(numbers::invalid_unsigned_int)
+ , n_q_points_scalar(numbers::invalid_unsigned_int)
, mapping(&mapping)
, fe(&fe)
, update_flags(update_flags)
, mapping_info_on_the_fly(
- std::make_unique<NonMatching::MappingInfo<dim, spacedim>>(mapping,
- update_flags))
+ std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
+ mapping,
+ update_flags))
, mapping_info(mapping_info_on_the_fly.get())
, current_cell_index(numbers::invalid_unsigned_int)
, current_face_number(numbers::invalid_unsigned_int)
template <int n_components, int dim, int spacedim, typename Number>
FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
- NonMatching::MappingInfo<dim, spacedim> &mapping_info,
- const FiniteElement<dim> & fe,
- const unsigned int first_selected_component)
+ NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+ const FiniteElement<dim> & fe,
+ const unsigned int first_selected_component)
: n_q_points(numbers::invalid_unsigned_int)
- , n_filled_lanes_last_batch(numbers::invalid_unsigned_int)
+ , n_q_points_scalar(numbers::invalid_unsigned_int)
, mapping(&mapping_info.get_mapping())
, fe(&fe)
, update_flags(mapping_info.get_update_flags())
FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
FEPointEvaluation<n_components_, dim, spacedim, Number> &other) noexcept
: n_q_points(other.n_q_points)
- , n_filled_lanes_last_batch(numbers::invalid_unsigned_int)
+ , n_q_points_scalar(other.n_q_points_scalar)
, mapping(other.mapping)
, fe(other.fe)
, poly(other.poly)
, fe_values(other.fe_values)
, mapping_info_on_the_fly(
other.mapping_info_on_the_fly ?
- std::make_unique<NonMatching::MappingInfo<dim, spacedim>>(
+ std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
*mapping,
update_flags) :
nullptr)
FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
FEPointEvaluation<n_components_, dim, spacedim, Number> &&other) noexcept
: n_q_points(other.n_q_points)
- , n_filled_lanes_last_batch(numbers::invalid_unsigned_int)
+ , n_q_points_scalar(other.n_q_points_scalar)
, mapping(other.mapping)
, fe(other.fe)
, poly(other.poly)
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 unsigned int n_q_points_unvectorized = unit_points.size();
-
- if (std::is_same<ScalarNumber, Number>::value)
- {
- const_cast<unsigned int &>(n_q_points) = unit_points.size();
- }
- else
- {
- const unsigned int n_lanes =
- internal::VectorizedArrayTrait<VectorizedArrayType>::width;
- const_cast<unsigned int &>(n_filled_lanes_last_batch) =
- n_q_points_unvectorized % n_lanes;
- const_cast<unsigned int &>(n_q_points) =
- n_q_points_unvectorized / n_lanes;
- if (n_filled_lanes_last_batch > 0)
- ++const_cast<unsigned int &>(n_q_points);
- }
+ const_cast<unsigned int &>(n_q_points_scalar) =
+ mapping_info->get_n_q_points_unvectorized(current_cell_index,
+ current_face_number);
+
+ const_cast<unsigned int &>(n_q_points) =
+ n_q_points_scalar / n_lanes_user_interface +
+ (n_q_points_scalar % n_lanes_user_interface > 0 ? 1 : 0);
+
+ // set unit point pointer
+ const unsigned int unit_point_offset =
+ mapping_info->compute_unit_point_index_offset(current_cell_index,
+ current_face_number);
+ unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
+
+ // set data pointers
+ const UpdateFlags update_flags_mapping =
+ mapping_info->get_update_flags_mapping();
+ const unsigned int data_offset =
+ mapping_info->compute_data_index_offset(current_cell_index,
+ current_face_number);
+ if (update_flags_mapping & UpdateFlags::update_quadrature_points)
+ real_point_ptr = mapping_info->get_real_point(data_offset);
+ if (update_flags_mapping & UpdateFlags::update_jacobians)
+ jacobian_ptr = mapping_info->get_jacobian(data_offset);
+ if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
+ inverse_jacobian_ptr = mapping_info->get_inverse_jacobian(data_offset);
+ if (update_flags_mapping & UpdateFlags::update_normal_vectors)
+ normal_ptr = mapping_info->get_normal_vector(data_offset);
+ if (update_flags_mapping & UpdateFlags::update_JxW_values)
+ JxW_ptr = mapping_info->get_JxW(data_offset);
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)
+ if (fast_path && !polynomials_are_hat_functions)
{
- const std::size_t n_points = unit_points.size();
- const std::size_t n_lanes = VectorizedArrayType::size();
const std::size_t n_batches =
- n_points / n_lanes + (n_points % n_lanes > 0 ? 1 : 0);
+ n_q_points_scalar / n_lanes_internal +
+ (n_q_points_scalar % n_lanes_internal > 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, VectorizedArrayType> 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);
- }
+ for (unsigned int qb = 0; qb < n_batches; ++qb)
+ internal::compute_values_of_array(make_array_view(shapes,
+ qb * n_shapes,
+ n_shapes),
+ poly,
+ unit_point_ptr[qb]);
}
is_reinitialized = true;
solution_renumbered.resize(dofs_per_component);
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, ScalarNumber>::read_value(
- solution_values[renumber[(component_in_base_element + comp) *
- dofs_per_component +
- i]],
- comp,
- solution_renumbered[i]);
+ ETT::read_value(
+ solution_values[renumber[(component_in_base_element + comp) *
+ dofs_per_component +
+ i]],
+ comp,
+ solution_renumbered[i]);
// unit gradients are currently only implemented with the fast tensor
// path
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 =
- internal::VectorizedArrayTrait<VectorizedArrayType>::width;
- const std::size_t stride =
- internal::FEPointEvaluation::VectorizedPointsTypeTraits<Number>::stride;
-
// loop over quadrature batches qb / points q
- for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+ const unsigned int n_shapes = poly.size();
+ for (unsigned int qb = 0, q = 0; q < n_q_points_scalar;
+ ++qb, q += n_lanes_internal)
{
// compute
- const unsigned int n_shapes = poly.size();
- const auto val_and_grad = [&]() {
- if (polynomials_are_hat_functions)
- {
- // convert quadrature points to vectorized format
- Point<dim, VectorizedArrayType> vectorized_points;
- for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v)
- for (unsigned int d = 0; d < dim; ++d)
- vectorized_points[d][v] = unit_points[q + v][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<
+ const auto val_and_grad =
+ polynomials_are_hat_functions ?
+ internal::evaluate_tensor_product_value_and_gradient_linear(
+ poly, solution_renumbered, unit_point_ptr[qb]) :
+ internal::evaluate_tensor_product_value_and_gradient_shapes<
dim,
scalar_value_type,
- VectorizedArrayType>(make_array_view(shapes.begin() + qb * n_shapes,
- shapes.begin() +
- (qb * n_shapes + n_shapes)),
- poly.size(),
+ VectorizedArrayType>(make_array_view(shapes,
+ qb * n_shapes,
+ n_shapes),
+ n_shapes,
solution_renumbered);
- }();
if (evaluation_flags & EvaluationFlags::values)
{
- for (unsigned int v = 0; v < stride && q + v < n_points; ++v)
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::set_value(
- val_and_grad.first, v, values[qb * stride + v]);
+ for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v)
+ ETT::set_value(val_and_grad.first, v, values[qb * stride + v]);
}
if (evaluation_flags & EvaluationFlags::gradients)
{
Assert(update_flags & update_gradients ||
update_flags & update_inverse_jacobians,
ExcNotInitialized());
- if (std::is_same<Number, ScalarNumber>::value)
- {
- // convert back to standard format
- for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v)
- {
- internal::FEPointEvaluation::EvaluatorTypeTraits<
- dim,
- n_components,
- Number>::set_gradient(val_and_grad.second,
- v,
- unit_gradients[q + v]);
- gradients[q + v] = apply_transformation(
- mapping_data.inverse_jacobians[q + v].transpose(),
- unit_gradients[q + v]);
- }
- }
- else
+
+ for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v)
{
- // convert to vectorized format
- DerivativeForm<1, spacedim, dim, Number>
- vectorized_inverse_transposed_jacobians;
- for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v)
- for (unsigned int d = 0; d < dim; ++d)
- for (unsigned int s = 0; s < spacedim; ++s)
- internal::VectorizedArrayTrait<Number>::get(
- vectorized_inverse_transposed_jacobians[s][d], v) =
- mapping_data.inverse_jacobians[q + v].transpose()[s][d];
-
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::set_gradient(
- val_and_grad.second, 0, unit_gradients[qb]);
- gradients[qb] =
- apply_transformation(vectorized_inverse_transposed_jacobians,
- unit_gradients[qb]);
+ const unsigned int offset = qb * stride + v;
+ ETT::set_gradient(val_and_grad.second, v, unit_gradients[offset]);
+ gradients[offset] =
+ apply_transformation(inverse_jacobian_ptr[offset].transpose(),
+ unit_gradients[offset]);
}
}
}
"Not initialized. Please call FEPointEvaluation::reinit()!"));
const std::size_t n_points = fe_values->get_quadrature().size();
- const std::size_t n_lanes = internal::VectorizedArrayTrait<Number>::width;
if (evaluation_flags & EvaluationFlags::values)
{
for (unsigned int d = 0; d < n_components; ++d)
if (nonzero_shape_function_component[i][d] &&
(fe->is_primitive(i) || fe->is_primitive()))
- for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
for (unsigned int v = 0;
- v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
- n_lanes);
+ v < n_lanes_user_interface && q + v < n_points;
++v)
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::access(
- values[qb],
- v,
- d,
- fe_values->shape_value(i, q + v) * value);
+ ETT::access(values[qb],
+ v,
+ d,
+ fe_values->shape_value(i, q + v) * value);
else if (nonzero_shape_function_component[i][d])
- for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
for (unsigned int v = 0;
- v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
- n_lanes);
+ v < n_lanes_user_interface && q + v < n_points;
++v)
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::access(
- values[qb],
- v,
- d,
- fe_values->shape_value_component(i, q + v, d) * value);
+ ETT::access(values[qb],
+ v,
+ d,
+ fe_values->shape_value_component(i, q + v, d) *
+ value);
}
}
for (unsigned int d = 0; d < n_components; ++d)
if (nonzero_shape_function_component[i][d] &&
(fe->is_primitive(i) || fe->is_primitive()))
- for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
for (unsigned int v = 0;
- v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
- n_lanes);
+ v < n_lanes_user_interface && q + v < n_points;
++v)
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::access(
- gradients[qb],
- v,
- d,
- fe_values->shape_grad(i, q + v) * value);
+ ETT::access(gradients[qb],
+ v,
+ d,
+ fe_values->shape_grad(i, q + v) * value);
else if (nonzero_shape_function_component[i][d])
- for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
for (unsigned int v = 0;
- v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
- n_lanes);
+ v < n_lanes_user_interface && q + v < n_points;
++v)
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::access(
- gradients[qb],
- v,
- d,
- fe_values->shape_grad_component(i, q + v, d) * value);
+ ETT::access(gradients[qb],
+ v,
+ d,
+ fe_values->shape_grad_component(i, q + v, d) *
+ value);
}
}
}
FEPointEvaluation<n_components, dim, spacedim, Number>::integrate_fast(
const ArrayView<ScalarNumber> & solution_values,
const EvaluationFlags::EvaluationFlags &integration_flags)
-
{
// fast path with tensor product integration
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,
- VectorizedArrayType>::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 =
- internal::VectorizedArrayTrait<VectorizedArrayType>::width;
- const std::size_t stride =
- internal::FEPointEvaluation::VectorizedPointsTypeTraits<Number>::stride;
+ solution_renumbered_vectorized.fill(vectorized_value_type());
// loop over quadrature batches qb / points q
- for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+ const unsigned int n_shapes = poly.size();
+ for (unsigned int qb = 0, q = 0; q < n_q_points_scalar;
+ ++qb, q += n_lanes_internal)
{
const bool incomplete_last_batch =
- (qb == (n_q_points - 1)) && (n_filled_lanes_last_batch > 0);
+ q + n_lanes_user_interface > n_q_points_scalar;
- typename internal::ProductTypeNoPoint<value_type,
- VectorizedArrayType>::type value =
- {};
- Tensor<1,
- dim,
- typename internal::ProductTypeNoPoint<value_type,
- VectorizedArrayType>::type>
- gradient;
+ vectorized_value_type value = {};
+ Tensor<1, dim, vectorized_value_type> gradient;
if (integration_flags & EvaluationFlags::values)
{
// zero out lanes of incomplete last quadrature point batch
if (incomplete_last_batch)
{
- for (unsigned int v = n_filled_lanes_last_batch; v < n_lanes; ++v)
- internal::FEPointEvaluation::EvaluatorTypeTraits<
- dim,
- n_components,
- Number>::set_zero_value(values[qb], v);
+ const unsigned int n_filled_lanes_last_batch =
+ n_q_points_scalar % n_lanes_internal;
+ for (unsigned int v = n_filled_lanes_last_batch;
+ v < n_lanes_internal;
+ ++v)
+ ETT::set_zero_value(values[qb], v);
}
- for (unsigned int v = 0; v < stride && q + v < n_points; ++v)
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::get_value(
- value, v, values[qb * stride + v]);
+
+ for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v)
+ ETT::get_value(value, v, values[qb * stride + v]);
}
if (integration_flags & EvaluationFlags::gradients)
{
- if (std::is_same<Number, ScalarNumber>::value)
+ // zero out lanes of incomplete last quadrature point batch
+ if (incomplete_last_batch)
{
- for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v)
- {
- gradients[q + v] =
- apply_transformation(mapping_data.inverse_jacobians[q + v],
- gradients[q + v]);
- internal::FEPointEvaluation::EvaluatorTypeTraits<
- dim,
- n_components,
- Number>::get_gradient(gradient, v, gradients[q + v]);
- }
+ const unsigned int n_filled_lanes_last_batch =
+ n_q_points_scalar % n_lanes_internal;
+ for (unsigned int v = n_filled_lanes_last_batch;
+ v < n_lanes_internal;
+ ++v)
+ ETT::set_zero_gradient(gradients[qb], v);
}
- else
+
+ for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v)
{
- // convert to vectorized format
- DerivativeForm<1, spacedim, dim, Number>
- vectorized_inverse_jacobians;
- for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v)
- for (unsigned int d = 0; d < dim; ++d)
- for (unsigned int s = 0; s < spacedim; ++s)
- internal::VectorizedArrayTrait<Number>::get(
- vectorized_inverse_jacobians[s][d], v) =
- mapping_data.inverse_jacobians[q + v][s][d];
-
- // zero out lanes of incomplete last quadrature point batch
- if (incomplete_last_batch)
- {
- for (unsigned int v = n_filled_lanes_last_batch; v < n_lanes;
- ++v)
- internal::FEPointEvaluation::EvaluatorTypeTraits<
- dim,
- n_components,
- Number>::set_zero_gradient(gradients[qb], v);
- }
-
- gradients[qb] = apply_transformation(vectorized_inverse_jacobians,
- gradients[qb]);
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::get_gradient(
- gradient, 0, gradients[qb]);
+ const unsigned int offset = qb * stride + v;
+ ETT::get_gradient(
+ gradient,
+ v,
+ apply_transformation(inverse_jacobian_ptr[offset],
+ gradients[offset]));
}
}
// compute
- const unsigned int n_shapes = poly.size();
if (polynomials_are_hat_functions)
- {
- // convert to vectorized format
- Point<dim, VectorizedArrayType> vectorized_points;
- for (unsigned int j = 0; j < n_lanes && q + j < n_points; ++j)
- for (unsigned int d = 0; d < dim; ++d)
- vectorized_points[d][j] = unit_points[q + j][d];
-
- internal::integrate_add_tensor_product_value_and_gradient_linear(
- poly,
- value,
- gradient,
- solution_renumbered_vectorized,
- vectorized_points);
- }
+ internal::integrate_add_tensor_product_value_and_gradient_linear(
+ poly,
+ value,
+ gradient,
+ solution_renumbered_vectorized,
+ unit_point_ptr[qb]);
else
internal::integrate_add_tensor_product_value_and_gradient_shapes<
dim,
VectorizedArrayType,
- typename internal::ProductTypeNoPoint<value_type,
- VectorizedArrayType>::type>(
- make_array_view(shapes.begin() + qb * n_shapes,
- shapes.begin() + (qb * n_shapes + n_shapes)),
- n_shapes,
- value,
- gradient,
- solution_renumbered_vectorized);
+ vectorized_value_type>(make_array_view(shapes,
+ qb * n_shapes,
+ n_shapes),
+ n_shapes,
+ value,
+ gradient,
+ solution_renumbered_vectorized);
}
// add between the lanes and write into the result
for (unsigned int i = 0; i < dofs_per_component; ++i)
{
VectorizedArrayType result;
- internal::FEPointEvaluation::EvaluatorTypeTraits<
- dim,
- n_components,
- VectorizedArrayType>::write_value(result,
- comp,
- solution_renumbered_vectorized[i]);
- for (unsigned int lane = n_lanes / 2; lane > 0; lane /= 2)
+ ETT::write_value(result, comp, solution_renumbered_vectorized[i]);
+ for (unsigned int lane = n_lanes_internal / 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];
std::fill(solution_values.begin(), solution_values.end(), 0.0);
const std::size_t n_points = fe_values->get_quadrature().size();
- const std::size_t n_lanes = internal::VectorizedArrayTrait<Number>::width;
if (integration_flags & EvaluationFlags::values)
{
for (unsigned int d = 0; d < n_components; ++d)
if (nonzero_shape_function_component[i][d] &&
(fe->is_primitive(i) || fe->is_primitive()))
- for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
for (unsigned int v = 0;
- v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
- n_lanes);
+ v < n_lanes_user_interface && q + v < n_points;
++v)
- solution_values[i] +=
- fe_values->shape_value(i, q + v) *
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::access(
- values[qb], v, d);
+ solution_values[i] += fe_values->shape_value(i, q + v) *
+ ETT::access(values[qb], v, d);
else if (nonzero_shape_function_component[i][d])
- for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
for (unsigned int v = 0;
- v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
- n_lanes);
+ v < n_lanes_user_interface && q + v < n_points;
++v)
solution_values[i] +=
fe_values->shape_value_component(i, q + v, d) *
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::access(
- values[qb], v, d);
+ ETT::access(values[qb], v, d);
}
}
for (unsigned int d = 0; d < n_components; ++d)
if (nonzero_shape_function_component[i][d] &&
(fe->is_primitive(i) || fe->is_primitive()))
- for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
for (unsigned int v = 0;
- v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
- n_lanes);
+ v < n_lanes_user_interface && q + v < n_points;
++v)
- solution_values[i] +=
- fe_values->shape_grad(i, q + v) *
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::access(
- gradients[qb], v, d);
+ solution_values[i] += fe_values->shape_grad(i, q + v) *
+ ETT::access(gradients[qb], v, d);
else if (nonzero_shape_function_component[i][d])
- for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
for (unsigned int v = 0;
- v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
- n_lanes);
+ v < n_lanes_user_interface && q + v < n_points;
++v)
solution_values[i] +=
fe_values->shape_grad_component(i, q + v, d) *
- internal::FEPointEvaluation::
- EvaluatorTypeTraits<dim, n_components, Number>::access(
- gradients[qb], v, d);
+ ETT::access(gradients[qb], v, d);
}
}
}
FEPointEvaluation<n_components, dim, spacedim, Number>::jacobian(
const unsigned int point_index) const
{
- const auto &mapping_data =
- mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
- const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
- DerivativeForm<1, dim, spacedim, Number> vectorized_jacobian;
- for (unsigned int v = 0;
- v < n_lanes && point_index * n_lanes + v < mapping_data.jacobians.size();
- ++v)
- for (unsigned int d = 0; d < dim; ++d)
- for (unsigned int s = 0; s < spacedim; ++s)
- internal::VectorizedArrayTrait<Number>::get(vectorized_jacobian[d][s],
- v) =
- mapping_data.jacobians[point_index * n_lanes + v][d][s];
- return vectorized_jacobian;
+ AssertIndexRange(point_index, n_q_points);
+ Assert(jacobian_ptr != nullptr,
+ internal::FEPointEvaluation::
+ ExcFEPointEvaluationAccessToUninitializedMappingField(
+ "update_jacobians"));
+ return jacobian_ptr[point_index];
}
FEPointEvaluation<n_components, dim, spacedim, Number>::inverse_jacobian(
const unsigned int point_index) const
{
- const auto &mapping_data =
- mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
- const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
- DerivativeForm<1, dim, spacedim, Number> vectorized_inverse_jacobian;
- for (unsigned int v = 0;
- v < n_lanes &&
- point_index * n_lanes + v < mapping_data.inverse_jacobians.size();
- ++v)
- for (unsigned int d = 0; d < dim; ++d)
- for (unsigned int s = 0; s < spacedim; ++s)
- internal::VectorizedArrayTrait<Number>::get(
- vectorized_inverse_jacobian[d][s], v) =
- mapping_data.inverse_jacobians[point_index * n_lanes + v][d][s];
- return vectorized_inverse_jacobian;
+ AssertIndexRange(point_index, n_q_points);
+ Assert(inverse_jacobian_ptr != nullptr,
+ internal::FEPointEvaluation::
+ ExcFEPointEvaluationAccessToUninitializedMappingField(
+ "update_inverse_jacobians"));
+ return inverse_jacobian_ptr[point_index];
}
FEPointEvaluation<n_components, dim, spacedim, Number>::JxW(
const unsigned int point_index) const
{
- const auto &mapping_data =
- mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
- const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
- Number vectorized_JxW = 0;
- for (unsigned int v = 0; v < n_lanes && point_index * n_lanes + v <
- mapping_data.JxW_values.size();
- ++v)
- internal::VectorizedArrayTrait<Number>::get(vectorized_JxW, v) =
- mapping_data.JxW_values[point_index * n_lanes + v];
- return vectorized_JxW;
+ AssertIndexRange(point_index, n_q_points);
+ Assert(JxW_ptr != nullptr,
+ internal::FEPointEvaluation::
+ ExcFEPointEvaluationAccessToUninitializedMappingField(
+ "update_JxW_values"));
+ return JxW_ptr[point_index];
}
+
template <int n_components, int dim, int spacedim, typename Number>
inline Tensor<1, spacedim, Number>
FEPointEvaluation<n_components, dim, spacedim, Number>::normal_vector(
const unsigned int point_index) const
{
- const auto &mapping_data =
- mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
- const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
- Tensor<1, spacedim, Number> vectorized_normal_vectors;
- for (unsigned int v = 0; v < n_lanes && point_index * n_lanes + v <
- mapping_data.normal_vectors.size();
- ++v)
- for (unsigned int s = 0; s < spacedim; ++s)
- internal::VectorizedArrayTrait<Number>::get(vectorized_normal_vectors[s],
- v) =
- mapping_data.normal_vectors[point_index * n_lanes + v][s];
- return vectorized_normal_vectors;
+ AssertIndexRange(point_index, n_q_points);
+ Assert(normal_ptr != nullptr,
+ internal::FEPointEvaluation::
+ ExcFEPointEvaluationAccessToUninitializedMappingField(
+ "update_normal_vectors"));
+ return normal_ptr[point_index];
}
FEPointEvaluation<n_components, dim, spacedim, Number>::real_point(
const unsigned int point_index) const
{
- const auto &mapping_data =
- mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
- const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
- Point<spacedim, Number> vectorized_real_point;
- for (unsigned int v = 0;
- v < n_lanes &&
- point_index * n_lanes + v < mapping_data.quadrature_points.size();
- ++v)
- for (unsigned int s = 0; s < spacedim; ++s)
- internal::VectorizedArrayTrait<Number>::get(vectorized_real_point[s], v) =
- mapping_data.quadrature_points[point_index * n_lanes + v][s];
- return vectorized_real_point;
+ AssertIndexRange(point_index, n_q_points);
+ Assert(real_point_ptr != nullptr,
+ internal::FEPointEvaluation::
+ ExcFEPointEvaluationAccessToUninitializedMappingField(
+ "update_quadrature_points"));
+ return real_point_ptr[point_index];
}
FEPointEvaluation<n_components, dim, spacedim, Number>::unit_point(
const unsigned int point_index) const
{
- const auto unit_points =
- mapping_info->get_unit_points(current_cell_index, current_face_number);
- const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
- Point<spacedim, Number> vectorized_unit_point;
- for (unsigned int v = 0;
- v < n_lanes && point_index * n_lanes + v < unit_points.size();
- ++v)
- for (unsigned int d = 0; d < dim; ++d)
- internal::VectorizedArrayTrait<Number>::get(vectorized_unit_point[d], v) =
- unit_points[point_index * n_lanes + v][d];
- return vectorized_unit_point;
+ AssertIndexRange(point_index, n_q_points);
+ Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
+ Point<dim, Number> unit_point;
+ for (unsigned int d = 0; d < dim; ++d)
+ unit_point[d] = internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+ unit_point_ptr[point_index / stride][d], point_index % stride);
+ return unit_point;
}
template <int dim, typename Number>
inline void
compute_values_of_array(
- ArrayView<dealii::ndarray<Number, 2, dim>> & shapes,
+ ArrayView<dealii::ndarray<Number, 2, dim>> shapes,
const std::vector<Polynomials::Polynomial<double>> &poly,
const Point<dim, Number> & p)
{
namespace NonMatching
{
+ namespace internal
+ {
+ template <int dim, int spacedim = dim>
+ class ComputeMappingDataHelper
+ {
+ using MappingData =
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
+ spacedim>;
+
+ public:
+ static UpdateFlags
+ required_update_flags(
+ const SmartPointer<const Mapping<dim, spacedim>> mapping,
+ const UpdateFlags & update_flags)
+ {
+ return mapping->requires_update_flags(update_flags);
+ }
+
+ static void
+ compute_mapping_data_for_quadrature(
+ const SmartPointer<const Mapping<dim, spacedim>> mapping,
+ const UpdateFlags & update_flags_mapping,
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ CellSimilarity::Similarity &cell_similarity,
+ const Quadrature<dim> & quadrature,
+ std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+ internal_mapping_data,
+ MappingData &mapping_data)
+ {
+ mapping_data.initialize(quadrature.size(), update_flags_mapping);
+
+ // reuse internal_mapping_data for MappingQ to avoid memory allocations
+ if (const MappingQ<dim, spacedim> *mapping_q =
+ dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+ {
+ (void)mapping_q;
+ auto &data =
+ dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+ *internal_mapping_data);
+ data.initialize(update_flags_mapping,
+ quadrature,
+ quadrature.size());
+ }
+ else
+ {
+ internal_mapping_data =
+ mapping->get_data(update_flags_mapping, quadrature);
+ }
+
+ cell_similarity = mapping->fill_fe_values(cell,
+ cell_similarity,
+ quadrature,
+ *internal_mapping_data,
+ mapping_data);
+ }
+
+
+
+ static void
+ compute_mapping_data_for_immersed_surface_quadrature(
+ const SmartPointer<const Mapping<dim, spacedim>> mapping,
+ const UpdateFlags & update_flags_mapping,
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const ImmersedSurfaceQuadrature<dim> & quadrature,
+ std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+ internal_mapping_data,
+ MappingData &mapping_data)
+ {
+ mapping_data.initialize(quadrature.size(), update_flags_mapping);
+
+ // reuse internal_mapping_data for MappingQ to avoid memory allocations
+ if (dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+ {
+ auto &data =
+ dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+ *internal_mapping_data);
+ data.initialize(update_flags_mapping,
+ quadrature,
+ quadrature.size());
+ }
+ else
+ {
+ internal_mapping_data =
+ mapping->get_data(update_flags_mapping, quadrature);
+ }
+
+ mapping->fill_fe_immersed_surface_values(cell,
+ quadrature,
+ *internal_mapping_data,
+ mapping_data);
+ }
+
+
+
+ static void
+ compute_mapping_data_for_face_quadrature(
+ const SmartPointer<const Mapping<dim, spacedim>> mapping,
+ const UpdateFlags & update_flags_mapping,
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature<dim - 1> & quadrature,
+ std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+ internal_mapping_data,
+ MappingData &mapping_data)
+ {
+ mapping_data.initialize(quadrature.size(), update_flags_mapping);
+
+ // reuse internal_mapping_data for MappingQ to avoid memory allocations
+ if (const MappingQ<dim, spacedim> *mapping_q =
+ dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+ {
+ auto &data =
+ dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+ *internal_mapping_data);
+ data.initialize_face(update_flags_mapping,
+ QProjector<dim>::project_to_oriented_face(
+ ReferenceCells::get_hypercube<dim>(),
+ quadrature,
+ face_no,
+ cell->face_orientation(face_no),
+ cell->face_flip(face_no),
+ cell->face_rotation(face_no)),
+ quadrature.size());
+
+ mapping_q->fill_mapping_data_for_face_quadrature(
+ cell, face_no, quadrature, *internal_mapping_data, mapping_data);
+ }
+ else
+ {
+ auto internal_mapping_data =
+ mapping->get_face_data(update_flags_mapping,
+ hp::QCollection<dim - 1>(quadrature));
+
+ mapping->fill_fe_face_values(cell,
+ face_no,
+ hp::QCollection<dim - 1>(quadrature),
+ *internal_mapping_data,
+ mapping_data);
+ }
+ }
+ };
+ } // namespace internal
+
/**
* This class provides the mapping information computation and mapping data
* storage to be used together with FEPointEvaluation.
+ *
+ * MappingInfo is vectorized across quadrature points, which means data can be
+ * provided in vectorized format.
+ *
+ * Two different modes are available: partially vectorized and fully
+ * vectorized across quadrature points. Partially vectorized (scalar Number
+ * template argument) means the computed mapping data is provided in scalar
+ * format and only unit points are provided vectorized (for seamless
+ * interaction with FEPointEvaluation). Fully vectorized (vectorized Number
+ * template argument, e.g. VectorizedArray<double>) provides both mapping data
+ * and unit points in vectorized format. The Number template parameter of
+ * MappingInfo and FEPointEvaluation has to be identical.
*/
- template <int dim, int spacedim = dim>
+ template <int dim, int spacedim = dim, typename Number = double>
class MappingInfo : public Subscriptor
{
public:
- using MappingData =
- dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
- spacedim>;
+ /**
+ * The VectorizedArray type the unit points are stored in inside this class.
+ */
+ using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
+ Number>::vectorized_value_type;
/**
* Constructor.
const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
/**
- * Getter function for current unit points.
- *
- * @p cell_index and @p face_number are the indices
- * into the compressed, CRS like data storage of unit points.
- *
- * If you have initialized this object with reinit_cells() you can access
- * the stored unit points of the cell with the respective @p cell_index
- * (and the default argument for @p face_number).
- *
- * If you have initialized this object with reinit_faces() you can access
- * the stored unit points of the faces on the cell with the respective @p cell_index
- * and the respective local @p face_number.
- *
- * If you have initialized this object with reinit() you can access the
- * stored unit points of a single cell with the default arguments.
- *
- * The correct state of this object is checked in this call (in debug mode).
+ * Getter function for unit points. The offset can be obtained with
+ * compute_unit_point_index_offset().
+ */
+ const Point<dim, VectorizedArrayType> *
+ get_unit_point(const unsigned int offset) const;
+
+ /**
+ * Getter function for Jacobians. The offset can be obtained with
+ * compute_data_index_offset().
+ */
+ const DerivativeForm<1, dim, spacedim, Number> *
+ get_jacobian(const unsigned int offset) const;
+
+ /**
+ * Getter function for inverse Jacobians. The offset can be obtained with
+ * compute_data_index_offset().
*/
- const ArrayView<const Point<dim>>
- get_unit_points(
- const unsigned int cell_index = numbers::invalid_unsigned_int,
- const unsigned int face_number = numbers::invalid_unsigned_int) const;
+ const DerivativeForm<1, spacedim, dim, Number> *
+ get_inverse_jacobian(const unsigned int offset) const;
/**
- * Getter function for computed mapping data. This function accesses
- * internal data and is therefore not a stable interface.
+ * Getter function for normal vectors. The offset can be obtained with
+ * compute_data_index_offset().
*/
- const MappingData &
- get_mapping_data(
- const unsigned int cell_index = numbers::invalid_unsigned_int,
- const unsigned int face_number = numbers::invalid_unsigned_int) const;
+ const Tensor<1, spacedim, Number> *
+ get_normal_vector(const unsigned int offset) const;
+
+ /**
+ * Getter function for Jacobian times quadrature weight (JxW). The offset
+ * can be obtained with compute_data_index_offset().
+ */
+ const Number *
+ get_JxW(const unsigned int offset) const;
+
+ /**
+ * Getter function for real points. The offset can be obtained with
+ * compute_data_index_offset().
+ */
+ const Point<dim, Number> *
+ get_real_point(const unsigned int offset) const;
/**
* Getter function for underlying mapping.
UpdateFlags
get_update_flags() const;
+ /**
+ * Getter function for the mapping update flags.
+ */
+ UpdateFlags
+ get_update_flags_mapping() const;
+
/**
* Connects to is_reinitialized().
*/
boost::signals2::connection
connect_is_reinitialized(const std::function<void()> &set_is_reinitialized);
+ /**
+ * Compute the unit points index offset for the current cell/face.
+ */
+ unsigned int
+ compute_unit_point_index_offset(const unsigned int cell_index,
+ const unsigned int face_number) const;
+
+ /**
+ * Compute the data index offset for the current cell/face.
+ */
+ unsigned int
+ compute_data_index_offset(const unsigned int cell_index,
+ const unsigned int face_number) const;
+
+ /**
+ * Get number of unvectorized quadrature points.
+ */
+ unsigned int
+ get_n_q_points_unvectorized(const unsigned int cell_index,
+ const unsigned int face_number) const;
+
private:
+ using MappingData =
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
+ spacedim>;
+
/**
- * Enum class for reinitialized states.
+ * Compute number of quadrature point batches depending on NumberType.
*/
- enum class State
- {
- invalid,
- single_cell,
- cell_vector,
- faces_on_cells_in_vector
- };
+ template <typename NumberType>
+ unsigned int
+ compute_n_q_points(const unsigned int n_q_points_unvectorized);
/**
- * Enum class that stores the currently initialized state
- * upon the last call of reinit().
+ * Resize the unit_point data field.
*/
- State state;
+ void
+ resize_unit_points(const unsigned int n_unit_point_batches);
/**
- * Compute the mapping related data for the given @p mapping,
- * @p cell and @p unit_points that is required by the FEPointEvaluation
- * class.
+ * Resize the mapping data fields.
*/
void
- compute_mapping_data_for_quadrature(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- CellSimilarity::Similarity &cell_similarity,
- const Quadrature<dim> & quadrature,
- MappingData & mapping_data);
+ resize_data_fields(const unsigned int n_data_point_batches);
/**
- * Compute the mapping related data for the given @p mapping,
- * @p cell and @p quadrature that is required by the FEPointEvaluation
- * class.
+ * Store the unit points.
*/
void
- compute_mapping_data_for_immersed_surface_quadrature(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const ImmersedSurfaceQuadrature<dim> & quadrature,
- MappingData & mapping_data);
+ store_unit_points(const unsigned int unit_points_index_offset,
+ const unsigned int n_q_points,
+ const unsigned int n_q_points_unvectorized,
+ const std::vector<Point<dim>> &points);
/**
- * Compute the mapping related data for the given @p mapping, @p cell,
- * @p face_no and @p quadrature that is required by the FEPointEvaluation
- * class.
+ * Store the requested mapping data.
*/
void
- compute_mapping_data_for_face_quadrature(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
- MappingData & mapping_data);
+ store_mapping_data(const unsigned int unit_points_index_offset,
+ const unsigned int n_q_points,
+ const unsigned int n_q_points_unvectorized,
+ const MappingData &mapping_data);
+
+ /**
+ * Compute the compressed cell index.
+ */
+ unsigned int
+ compute_compressed_cell_index(const unsigned int cell_index) const;
+
+ /**
+ * Compute the geometry index offset of the current cell/face.
+ */
+ unsigned int
+ compute_geometry_index_offset(const unsigned int cell_index,
+ const unsigned int face_number) const;
+
+ /**
+ * Enum class for reinitialized states.
+ */
+ enum class State
+ {
+ invalid,
+ single_cell,
+ cell_vector,
+ faces_on_cells_in_vector
+ };
+
+ /**
+ * Enum class that stores the currently initialized state
+ * upon the last call of reinit().
+ */
+ State state;
/**
* The reference points specified at reinit().
+ *
+ * Indexed by @p unit_points_index.
*/
- std::vector<Point<dim>> unit_points;
+ AlignedVector<Point<dim, VectorizedArrayType>> unit_points;
/**
- * Offset to point to the first unit point of a cell/face
+ * Offset to point to the first unit point of a cell/face.
*/
- std::vector<unsigned int> unit_points_index;
+ AlignedVector<unsigned int> unit_points_index;
/**
* A pointer to the internal data of the underlying mapping.
*/
- std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+ std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
internal_mapping_data;
/**
UpdateFlags update_flags_mapping;
/**
- * The internal data container for mapping information. The implementation
- * is subject to future changes.
+ * Stores the index offset into the arrays @p JxW_values, @p jacobians,
+ * @p inverse_jacobians and @p normal_vectors.
*/
- std::vector<MappingData> mapping_data;
+ AlignedVector<unsigned int> data_index_offsets;
+
+ /**
+ * The storage of the Jacobian determinant times the quadrature weight on
+ * quadrature points.
+ *
+ * Indexed by @p data_index_offsets.
+ */
+ AlignedVector<Number> JxW_values;
+
+ /**
+ * Stores the normal vectors.
+ *
+ * Indexed by @p data_index_offsets.
+ */
+ AlignedVector<Tensor<1, spacedim, Number>> normal_vectors;
+
+ /**
+ * The storage of contravariant transformation on quadrature points, i.e.,
+ * the Jacobians of the transformation from the unit to the real cell.
+ *
+ * Indexed by @p data_index_offsets.
+ */
+ AlignedVector<DerivativeForm<1, dim, spacedim, Number>> jacobians;
+
+ /**
+ * The storage of covariant transformation on quadrature points, i.e.,
+ * the inverse Jacobians of the transformation from the
+ * unit to the real cell.
+ *
+ * Indexed by @p data_index_offsets.
+ */
+ AlignedVector<DerivativeForm<1, spacedim, dim, Number>> inverse_jacobians;
+
+ /**
+ * The mapped real points.
+ *
+ * Indexed by @p data_index_offsets.
+ */
+ AlignedVector<Point<spacedim, Number>> real_points;
+
+ /**
+ * Number of unvectorized unit points per geometric entity (cell/face).
+ */
+ std::vector<unsigned int> n_q_points_unvectorized;
/**
* Offset to point to the first element of a cell in internal data
// ----------------------- template functions ----------------------
- template <int dim, int spacedim>
- MappingInfo<dim, spacedim>::MappingInfo(const Mapping<dim> &mapping,
- const UpdateFlags update_flags)
+ template <int dim, int spacedim, typename Number>
+ MappingInfo<dim, spacedim, Number>::MappingInfo(
+ const Mapping<dim> &mapping,
+ const UpdateFlags update_flags)
: mapping(&mapping)
, update_flags(update_flags)
+ , update_flags_mapping(update_default)
{
- update_flags_mapping = update_default;
// translate update flags
- if (update_flags & update_jacobians || update_flags & update_JxW_values)
+ if (update_flags & update_jacobians)
update_flags_mapping |= update_jacobians;
if (update_flags & update_JxW_values)
update_flags_mapping |= update_JxW_values;
// always save quadrature points for now
update_flags_mapping |= update_quadrature_points;
+ update_flags_mapping =
+ internal::ComputeMappingDataHelper<dim, spacedim>::required_update_flags(
+ this->mapping, update_flags_mapping);
+
// construct internal_mapping_data for MappingQ to be able to reuse it in
// reinit() calls to avoid memory allocations
if (const MappingQ<dim, spacedim> *mapping_q =
- template <int dim, int spacedim>
+ template <int dim, int spacedim, typename Number>
void
- MappingInfo<dim, spacedim>::reinit(
+ MappingInfo<dim, spacedim, Number>::reinit(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const std::vector<Point<dim>> & unit_points_in)
{
- template <int dim, int spacedim>
+ template <int dim, int spacedim, typename Number>
void
- MappingInfo<dim, spacedim>::reinit(
+ MappingInfo<dim, spacedim, Number>::reinit(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const ArrayView<const Point<dim>> & unit_points_in)
{
- template <int dim, int spacedim>
+ template <int dim, int spacedim, typename Number>
void
- MappingInfo<dim, spacedim>::reinit(
+ MappingInfo<dim, spacedim, Number>::reinit(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const Quadrature<dim> & quadrature)
{
- unit_points = quadrature.get_points();
-
- mapping_data.resize(1);
- CellSimilarity::Similarity cell_similarity =
- CellSimilarity::Similarity::none;
- compute_mapping_data_for_quadrature(cell,
- cell_similarity,
- quadrature,
- mapping_data[0]);
+ n_q_points_unvectorized.resize(1);
+ n_q_points_unvectorized[0] = quadrature.size();
+
+ const unsigned int n_q_points =
+ compute_n_q_points<VectorizedArrayType>(n_q_points_unvectorized[0]);
+
+ const unsigned int n_q_points_data =
+ compute_n_q_points<Number>(n_q_points_unvectorized[0]);
+
+ // resize data vectors
+ resize_unit_points(n_q_points);
+ resize_data_fields(n_q_points_data);
+
+ // store unit points
+ store_unit_points(0,
+ n_q_points,
+ n_q_points_unvectorized[0],
+ quadrature.get_points());
+
+ // compute mapping data
+ MappingData mapping_data;
+ CellSimilarity::Similarity cell_similarity = CellSimilarity::none;
+ internal::ComputeMappingDataHelper<dim, spacedim>::
+ compute_mapping_data_for_quadrature(mapping,
+ update_flags_mapping,
+ cell,
+ cell_similarity,
+ quadrature,
+ internal_mapping_data,
+ mapping_data);
+
+ // store mapping data
+ store_mapping_data(0,
+ n_q_points_data,
+ n_q_points_unvectorized[0],
+ mapping_data);
state = State::single_cell;
is_reinitialized();
- template <int dim, int spacedim>
+ template <int dim, int spacedim, typename Number>
template <typename ContainerType>
void
- MappingInfo<dim, spacedim>::reinit_cells(
+ MappingInfo<dim, spacedim, Number>::reinit_cells(
const ContainerType & cell_iterator_range,
const std::vector<std::vector<Point<dim>>> &unit_points_vector,
const unsigned int n_unfiltered_cells)
- template <int dim, int spacedim>
+ template <int dim, int spacedim, typename Number>
template <typename ContainerType>
void
- MappingInfo<dim, spacedim>::reinit_cells(
+ MappingInfo<dim, spacedim, Number>::reinit_cells(
const ContainerType & cell_iterator_range,
const std::vector<Quadrature<dim>> &quadrature_vector,
const unsigned int n_unfiltered_cells)
std::distance(cell_iterator_range.begin(),
cell_iterator_range.end()));
+ n_q_points_unvectorized.reserve(n_cells);
+
// fill unit points index offset vector
unit_points_index.reserve(n_cells + 1);
unit_points_index.push_back(0);
+ data_index_offsets.reserve(n_cells + 1);
+ data_index_offsets.push_back(0);
for (const auto &quadrature : quadrature_vector)
- unit_points_index.push_back(unit_points_index.back() + quadrature.size());
+ {
+ const unsigned int n_points = quadrature.size();
+ n_q_points_unvectorized.push_back(n_points);
+
+ const unsigned int n_q_points =
+ compute_n_q_points<VectorizedArrayType>(n_points);
+ unit_points_index.push_back(unit_points_index.back() + n_q_points);
+
+ const unsigned int n_q_points_data =
+ compute_n_q_points<Number>(n_points);
+ data_index_offsets.push_back(data_index_offsets.back() +
+ n_q_points_data);
+ }
const unsigned int n_unit_points = unit_points_index.back();
+ const unsigned int n_data_points = data_index_offsets.back();
- unit_points.resize(n_unit_points);
- mapping_data.resize(n_cells);
+ // resize data vectors
+ resize_unit_points(n_unit_points);
+ resize_data_fields(n_data_points);
if (do_cell_index_compression)
cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
numbers::invalid_unsigned_int);
+
+ MappingData mapping_data;
CellSimilarity::Similarity cell_similarity =
CellSimilarity::Similarity::none;
unsigned int cell_index = 0;
for (const auto &cell : cell_iterator_range)
{
- auto it = unit_points.begin() + unit_points_index[cell_index];
- for (const auto &unit_point :
- quadrature_vector[cell_index].get_points())
- {
- *it = unit_point;
- ++it;
- }
-
- compute_mapping_data_for_quadrature(cell,
- cell_similarity,
- quadrature_vector[cell_index],
- mapping_data[cell_index]);
+ // store unit points
+ const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
+ n_q_points_unvectorized[cell_index]);
+ store_unit_points(unit_points_index[cell_index],
+ n_q_points,
+ n_q_points_unvectorized[cell_index],
+ quadrature_vector[cell_index].get_points());
+
+ // compute mapping data
+ internal::ComputeMappingDataHelper<dim, spacedim>::
+ compute_mapping_data_for_quadrature(mapping,
+ update_flags_mapping,
+ cell,
+ cell_similarity,
+ quadrature_vector[cell_index],
+ internal_mapping_data,
+ mapping_data);
+
+ // store mapping data
+ const unsigned int n_q_points_data =
+ compute_n_q_points<Number>(n_q_points_unvectorized[cell_index]);
+ store_mapping_data(data_index_offsets[cell_index],
+ n_q_points_data,
+ n_q_points_unvectorized[cell_index],
+ mapping_data);
if (do_cell_index_compression)
cell_index_to_compressed_cell_index[cell->active_cell_index()] =
- template <int dim, int spacedim>
+ template <int dim, int spacedim, typename Number>
template <typename Iterator>
void
- MappingInfo<dim, spacedim>::reinit_surface(
+ MappingInfo<dim, spacedim, Number>::reinit_surface(
const IteratorRange<Iterator> & cell_iterator_range,
const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
const unsigned int n_unfiltered_cells)
std::distance(cell_iterator_range.begin(),
cell_iterator_range.end()));
+ n_q_points_unvectorized.reserve(n_cells);
+
// fill unit points index offset vector
unit_points_index.reserve(n_cells + 1);
unit_points_index.push_back(0);
+ data_index_offsets.reserve(n_cells + 1);
+ data_index_offsets.push_back(0);
for (const auto &quadrature : quadrature_vector)
- unit_points_index.push_back(unit_points_index.back() +
- quadrature.get_points().size());
+ {
+ const unsigned int n_points = quadrature.size();
+ n_q_points_unvectorized.push_back(n_points);
+
+ const unsigned int n_q_points =
+ compute_n_q_points<VectorizedArrayType>(n_points);
+ unit_points_index.push_back(unit_points_index.back() + n_q_points);
+
+ const unsigned int n_q_points_data =
+ compute_n_q_points<Number>(n_points);
+ data_index_offsets.push_back(data_index_offsets.back() +
+ n_q_points_data);
+ }
const unsigned int n_unit_points = unit_points_index.back();
+ const unsigned int n_data_points = data_index_offsets.back();
- unit_points.resize(n_unit_points);
- mapping_data.resize(n_cells);
+ // resize data vectors
+ resize_unit_points(n_unit_points);
+ resize_data_fields(n_data_points);
if (do_cell_index_compression)
cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
numbers::invalid_unsigned_int);
+
+ MappingData mapping_data;
unsigned int cell_index = 0;
for (const auto &cell : cell_iterator_range)
{
const auto &quadrature = quadrature_vector[cell_index];
- auto it = unit_points.begin() + unit_points_index[cell_index];
- for (const auto &unit_point : quadrature.get_points())
- {
- *it = unit_point;
- ++it;
- }
-
- compute_mapping_data_for_immersed_surface_quadrature(
- cell, quadrature, mapping_data[cell_index]);
+ // store unit points
+ const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
+ n_q_points_unvectorized[cell_index]);
+ store_unit_points(unit_points_index[cell_index],
+ n_q_points,
+ n_q_points_unvectorized[cell_index],
+ quadrature_vector[cell_index].get_points());
+
+ // compute mapping data
+ internal::ComputeMappingDataHelper<dim, spacedim>::
+ compute_mapping_data_for_immersed_surface_quadrature(
+ mapping,
+ update_flags_mapping,
+ cell,
+ quadrature,
+ internal_mapping_data,
+ mapping_data);
+
+ // store mapping data
+ const unsigned int n_q_points_data =
+ compute_n_q_points<Number>(n_q_points_unvectorized[cell_index]);
+ store_mapping_data(data_index_offsets[cell_index],
+ n_q_points_data,
+ n_q_points_unvectorized[cell_index],
+ mapping_data);
if (do_cell_index_compression)
cell_index_to_compressed_cell_index[cell->active_cell_index()] =
- template <int dim, int spacedim>
+ template <int dim, int spacedim, typename Number>
template <typename Iterator>
void
- MappingInfo<dim, spacedim>::reinit_faces(
+ MappingInfo<dim, spacedim, Number>::reinit_faces(
const IteratorRange<Iterator> & cell_iterator_range,
const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector,
const unsigned int n_unfiltered_cells)
++cell_index;
}
+ n_q_points_unvectorized.reserve(n_faces);
+
// fill unit points index offset vector
unit_points_index.resize(n_faces + 1);
+ data_index_offsets.resize(n_faces + 1);
cell_index = 0;
unsigned int n_unit_points = 0;
+ unsigned int n_data_points = 0;
for (const auto &cell : cell_iterator_range)
{
for (const auto &f : cell->face_indices())
const unsigned int current_face_index =
cell_index_offset[cell_index] + f;
- unit_points_index[current_face_index] = n_unit_points;
- n_unit_points +=
- quadrature_vector[cell_index][f].get_points().size();
+ unit_points_index[current_face_index] = n_unit_points;
+ data_index_offsets[current_face_index] = n_data_points;
+
+ const unsigned int n_points =
+ quadrature_vector[cell_index][f].size();
+ n_q_points_unvectorized.push_back(n_points);
+
+ const unsigned int n_q_points =
+ compute_n_q_points<VectorizedArrayType>(n_points);
+ n_unit_points += n_q_points;
+
+ const unsigned int n_q_points_data =
+ compute_n_q_points<Number>(n_points);
+ n_data_points += n_q_points_data;
}
++cell_index;
}
- unit_points_index[n_faces] = n_unit_points;
+ unit_points_index[n_faces] = n_unit_points;
+ data_index_offsets[n_faces] = n_data_points;
// compress indices
if (do_cell_index_compression)
numbers::invalid_unsigned_int);
// fill unit points and mapping data for every face of all cells
- unit_points.resize(n_unit_points);
- mapping_data.resize(n_faces);
+ // resize data vectors
+ resize_unit_points(n_unit_points);
+ resize_data_fields(n_data_points);
+
+ MappingData mapping_data;
cell_index = 0;
QProjector<dim> q_projector;
for (const auto &cell : cell_iterator_range)
const unsigned int current_face_index =
cell_index_offset[cell_index] + f;
- auto it =
- unit_points.begin() + unit_points_index[current_face_index];
- for (const auto &unit_point : unit_points_on_cell)
- {
- *it = unit_point;
- ++it;
- }
-
- compute_mapping_data_for_face_quadrature(
- cell, f, quadrature_on_face, mapping_data[current_face_index]);
+ // store unit points
+ const unsigned int n_q_points =
+ compute_n_q_points<VectorizedArrayType>(
+ n_q_points_unvectorized[current_face_index]);
+ store_unit_points(unit_points_index[current_face_index],
+ n_q_points,
+ n_q_points_unvectorized[current_face_index],
+ quadrature_on_cell.get_points());
+
+ internal::ComputeMappingDataHelper<dim, spacedim>::
+ compute_mapping_data_for_face_quadrature(mapping,
+ update_flags_mapping,
+ cell,
+ f,
+ quadrature_on_face,
+ internal_mapping_data,
+ mapping_data);
+
+ const unsigned int n_q_points_data = compute_n_q_points<Number>(
+ n_q_points_unvectorized[current_face_index]);
+ store_mapping_data(data_index_offsets[current_face_index],
+ n_q_points_data,
+ n_q_points_unvectorized[current_face_index],
+ mapping_data);
}
if (do_cell_index_compression)
cell_index_to_compressed_cell_index[cell->active_cell_index()] =
- template <int dim, int spacedim>
- inline const ArrayView<const Point<dim>>
- MappingInfo<dim, spacedim>::get_unit_points(
+ template <int dim, int spacedim, typename Number>
+ unsigned int
+ MappingInfo<dim, spacedim, Number>::get_n_q_points_unvectorized(
const unsigned int cell_index,
const unsigned int face_number) const
{
Assert(state == State::single_cell,
ExcMessage(
"This mapping info is not reinitialized for a single cell!"));
- return unit_points;
+ return n_q_points_unvectorized[0];
+ }
+ else
+ {
+ return n_q_points_unvectorized[compute_geometry_index_offset(
+ cell_index, face_number)];
}
- else if (face_number == numbers::invalid_unsigned_int)
+ }
+
+
+
+ template <int dim, int spacedim, typename Number>
+ template <typename NumberType>
+ unsigned int
+ MappingInfo<dim, spacedim, Number>::compute_n_q_points(
+ const unsigned int n_q_points_unvectorized)
+ {
+ const unsigned int n_lanes =
+ dealii::internal::VectorizedArrayTrait<NumberType>::width;
+ const unsigned int n_filled_lanes_last_batch =
+ n_q_points_unvectorized % n_lanes;
+ unsigned int n_q_points = n_q_points_unvectorized / n_lanes;
+ if (n_filled_lanes_last_batch > 0)
+ ++n_q_points;
+ return n_q_points;
+ }
+
+
+
+ template <int dim, int spacedim, typename Number>
+ unsigned int
+ MappingInfo<dim, spacedim, Number>::compute_geometry_index_offset(
+ const unsigned int cell_index,
+ const unsigned int face_number) const
+ {
+ const unsigned int compressed_cell_index =
+ compute_compressed_cell_index(cell_index);
+ if (face_number == numbers::invalid_unsigned_int)
{
Assert(state == State::cell_vector,
ExcMessage(
"This mapping info is not reinitialized for a cell vector!"));
- if (do_cell_index_compression)
- {
- Assert(cell_index_to_compressed_cell_index[cell_index] !=
- numbers::invalid_unsigned_int,
- ExcMessage("Mapping info object was not initialized for this"
- " active cell index!"));
- const auto it_begin =
- unit_points.begin() +
- unit_points_index
- [cell_index_to_compressed_cell_index[cell_index]];
- const auto it_end =
- unit_points.begin() +
- unit_points_index
- [cell_index_to_compressed_cell_index[cell_index] + 1];
- return make_array_view(it_begin, it_end);
- }
- else
- {
- const auto it_begin =
- unit_points.begin() + unit_points_index[cell_index];
- const auto it_end =
- unit_points.begin() + unit_points_index[cell_index + 1];
- return make_array_view(it_begin, it_end);
- }
+ return compressed_cell_index;
}
- else if (cell_index != numbers::invalid_unsigned_int)
+ else
{
+ Assert(cell_index != numbers::invalid_unsigned_int,
+ ExcMessage(
+ "cell_index has to be set if face_number is specified!"));
Assert(state == State::faces_on_cells_in_vector,
ExcMessage("This mapping info is not reinitialized for faces"
" on cells in a vector!"));
- if (do_cell_index_compression)
- {
- Assert(
- cell_index_to_compressed_cell_index[cell_index] !=
- numbers::invalid_unsigned_int,
- ExcMessage(
- "Mapping info object was not initialized for this active cell index"
- " and corresponding face numbers!"));
- const unsigned int current_face_index =
- cell_index_offset
- [cell_index_to_compressed_cell_index[cell_index]] +
- face_number;
- const auto it_begin =
- unit_points.begin() + unit_points_index[current_face_index];
- const auto it_end =
- unit_points.begin() + unit_points_index[current_face_index + 1];
- return make_array_view(it_begin, it_end);
- }
- else
- {
- const unsigned int current_face_index =
- cell_index_offset[cell_index] + face_number;
- const auto it_begin =
- unit_points.begin() + unit_points_index[current_face_index];
- const auto it_end =
- unit_points.begin() + unit_points_index[current_face_index + 1];
- return make_array_view(it_begin, it_end);
- }
+ return cell_index_offset[compressed_cell_index] + face_number;
}
- else
- AssertThrow(
- false,
- ExcMessage(
- "cell_index has to be specified if face_number is specified!"));
}
- template <int dim, int spacedim>
- inline const typename MappingInfo<dim, spacedim>::MappingData &
- MappingInfo<dim, spacedim>::get_mapping_data(
- const unsigned int cell_index,
- const unsigned int face_number) const
+ template <int dim, int spacedim, typename Number>
+ unsigned int
+ MappingInfo<dim, spacedim, Number>::compute_compressed_cell_index(
+ const unsigned int cell_index) const
{
- if (cell_index == numbers::invalid_unsigned_int &&
- face_number == numbers::invalid_unsigned_int)
+ if (do_cell_index_compression)
{
- Assert(state == State::single_cell,
- ExcMessage(
- "This mapping info is not reinitialized for a single cell!"));
- return mapping_data[0];
+ Assert(cell_index_to_compressed_cell_index[cell_index] !=
+ numbers::invalid_unsigned_int,
+ ExcMessage("Mapping info object was not initialized for this"
+ " active cell index!"));
+ return cell_index_to_compressed_cell_index[cell_index];
}
- else if (face_number == numbers::invalid_unsigned_int)
+ else
+ return cell_index;
+ }
+
+
+ template <int dim, int spacedim, typename Number>
+ void
+ MappingInfo<dim, spacedim, Number>::store_unit_points(
+ const unsigned int unit_points_index_offset,
+ const unsigned int n_q_points,
+ const unsigned int n_q_points_unvectorized,
+ const std::vector<Point<dim>> &points)
+ {
+ const unsigned int n_lanes =
+ dealii::internal::VectorizedArrayTrait<VectorizedArrayType>::width;
+
+ for (unsigned int q = 0; q < n_q_points; ++q)
{
- Assert(state == State::cell_vector,
- ExcMessage(
- "This mapping info is not reinitialized for a cell vector!"));
- if (do_cell_index_compression)
- {
- Assert(cell_index_to_compressed_cell_index[cell_index] !=
- numbers::invalid_unsigned_int,
- ExcMessage("Mapping info object was not initialized for this"
- " active cell index!"));
- return mapping_data
- [cell_index_to_compressed_cell_index[cell_index]];
- }
- else
- return mapping_data[cell_index];
+ const unsigned int offset = unit_points_index_offset + q;
+ for (unsigned int v = 0;
+ v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
+ ++v)
+ for (unsigned int d = 0; d < dim; ++d)
+ dealii::internal::VectorizedArrayTrait<VectorizedArrayType>::get(
+ unit_points[offset][d], v) = points[q * n_lanes + v][d];
}
- else if (cell_index != numbers::invalid_unsigned_int)
+ }
+
+
+
+ template <int dim, int spacedim, typename Number>
+ void
+ MappingInfo<dim, spacedim, Number>::store_mapping_data(
+ const unsigned int unit_points_index_offset,
+ const unsigned int n_q_points,
+ const unsigned int n_q_points_unvectorized,
+ const MappingInfo::MappingData &mapping_data)
+ {
+ const unsigned int n_lanes =
+ dealii::internal::VectorizedArrayTrait<Number>::width;
+
+ for (unsigned int q = 0; q < n_q_points; ++q)
{
- Assert(state == State::faces_on_cells_in_vector,
- ExcMessage("This mapping info is not reinitialized for faces"
- " on cells in a vector!"));
- if (do_cell_index_compression)
+ const unsigned int offset = unit_points_index_offset + q;
+ for (unsigned int v = 0;
+ v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
+ ++v)
{
- Assert(
- cell_index_to_compressed_cell_index[cell_index] !=
- numbers::invalid_unsigned_int,
- ExcMessage(
- "Mapping info object was not initialized for this active cell index"
- " and corresponding face numbers!"));
- return mapping_data
- [cell_index_offset
- [cell_index_to_compressed_cell_index[cell_index]] +
- face_number];
+ if (update_flags_mapping & UpdateFlags::update_jacobians)
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int s = 0; s < spacedim; ++s)
+ dealii::internal::VectorizedArrayTrait<Number>::get(
+ jacobians[offset][d][s], v) =
+ mapping_data.jacobians[q * n_lanes + v][d][s];
+ if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int s = 0; s < spacedim; ++s)
+ dealii::internal::VectorizedArrayTrait<Number>::get(
+ inverse_jacobians[offset][s][d], v) =
+ mapping_data.inverse_jacobians[q * n_lanes + v][s][d];
+ if (update_flags_mapping & UpdateFlags::update_JxW_values)
+ dealii::internal::VectorizedArrayTrait<Number>::get(
+ JxW_values[offset], v) =
+ mapping_data.JxW_values[q * n_lanes + v];
+ if (update_flags_mapping & UpdateFlags::update_normal_vectors)
+ for (unsigned int s = 0; s < spacedim; ++s)
+ dealii::internal::VectorizedArrayTrait<Number>::get(
+ normal_vectors[offset][s], v) =
+ mapping_data.normal_vectors[q * n_lanes + v][s];
+ if (update_flags_mapping & UpdateFlags::update_quadrature_points)
+ for (unsigned int s = 0; s < spacedim; ++s)
+ dealii::internal::VectorizedArrayTrait<Number>::get(
+ real_points[offset][s], v) =
+ mapping_data.quadrature_points[q * n_lanes + v][s];
}
- else
- return mapping_data[cell_index_offset[cell_index] + face_number];
}
- else
- AssertThrow(
- false,
- ExcMessage(
- "cell_index has to be specified if face_number is specified!"));
}
- template <int dim, int spacedim>
- const Mapping<dim, spacedim> &
- MappingInfo<dim, spacedim>::get_mapping() const
+ template <int dim, int spacedim, typename Number>
+ void
+ MappingInfo<dim, spacedim, Number>::resize_unit_points(
+ const unsigned int n_unit_point_batches)
{
- return *mapping;
+ unit_points.resize(n_unit_point_batches);
}
- template <int dim, int spacedim>
- UpdateFlags
- MappingInfo<dim, spacedim>::get_update_flags() const
+ template <int dim, int spacedim, typename Number>
+ void
+ MappingInfo<dim, spacedim, Number>::resize_data_fields(
+ const unsigned int n_data_point_batches)
{
- return update_flags;
+ if (update_flags_mapping & UpdateFlags::update_jacobians)
+ jacobians.resize(n_data_point_batches);
+ if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
+ inverse_jacobians.resize(n_data_point_batches);
+ if (update_flags_mapping & UpdateFlags::update_JxW_values)
+ JxW_values.resize(n_data_point_batches);
+ if (update_flags_mapping & UpdateFlags::update_normal_vectors)
+ normal_vectors.resize(n_data_point_batches);
+ if (update_flags_mapping & UpdateFlags::update_quadrature_points)
+ real_points.resize(n_data_point_batches);
}
- template <int dim, int spacedim>
- boost::signals2::connection
- MappingInfo<dim, spacedim>::connect_is_reinitialized(
- const std::function<void()> &set_is_reinitialized)
+ template <int dim, int spacedim, typename Number>
+ inline const Point<
+ dim,
+ typename MappingInfo<dim, spacedim, Number>::VectorizedArrayType> *
+ MappingInfo<dim, spacedim, Number>::get_unit_point(
+ const unsigned int offset) const
{
- return is_reinitialized.connect(set_is_reinitialized);
+ return &unit_points[offset];
}
- template <int dim, int spacedim>
- void
- MappingInfo<dim, spacedim>::compute_mapping_data_for_quadrature(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- CellSimilarity::Similarity & cell_similarity,
- const Quadrature<dim> & quadrature,
- MappingData & mapping_data)
+ template <int dim, int spacedim, typename Number>
+ inline const Point<dim, Number> *
+ MappingInfo<dim, spacedim, Number>::get_real_point(
+ const unsigned int offset) const
{
- update_flags_mapping |=
- mapping->requires_update_flags(update_flags_mapping);
+ return &real_points[offset];
+ }
- mapping_data.initialize(quadrature.size(), update_flags_mapping);
- // reuse internal_mapping_data for MappingQ to avoid memory allocations
- if (const MappingQ<dim, spacedim> *mapping_q =
- dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+
+ template <int dim, int spacedim, typename Number>
+ unsigned int
+ MappingInfo<dim, spacedim, Number>::compute_unit_point_index_offset(
+ const unsigned int cell_index,
+ const unsigned int face_number) const
+ {
+ if (cell_index == numbers::invalid_unsigned_int &&
+ face_number == numbers::invalid_unsigned_int)
{
- (void)mapping_q;
- auto &data =
- dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
- *internal_mapping_data);
- data.initialize(update_flags_mapping, quadrature, quadrature.size());
+ Assert(state == State::single_cell,
+ ExcMessage(
+ "This mapping info is not reinitialized for a single cell!"));
+ return 0;
}
else
{
- internal_mapping_data =
- mapping->get_data(update_flags_mapping, quadrature);
+ const unsigned int offset =
+ compute_geometry_index_offset(cell_index, face_number);
+ return unit_points_index[offset];
}
-
- cell_similarity = mapping->fill_fe_values(
- cell, cell_similarity, quadrature, *internal_mapping_data, mapping_data);
}
- template <int dim, int spacedim>
- void
- MappingInfo<dim, spacedim>::
- compute_mapping_data_for_immersed_surface_quadrature(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const ImmersedSurfaceQuadrature<dim> & quadrature,
- MappingData & mapping_data)
+ template <int dim, int spacedim, typename Number>
+ unsigned int
+ MappingInfo<dim, spacedim, Number>::compute_data_index_offset(
+ const unsigned int cell_index,
+ const unsigned int face_number) const
{
- update_flags_mapping |=
- mapping->requires_update_flags(update_flags_mapping);
-
- mapping_data.initialize(quadrature.size(), update_flags_mapping);
-
- // reuse internal_mapping_data for MappingQ to avoid memory allocations
- if (const MappingQ<dim, spacedim> *mapping_q =
- dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+ if (cell_index == numbers::invalid_unsigned_int &&
+ face_number == numbers::invalid_unsigned_int)
{
- (void)mapping_q;
- auto &data =
- dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
- *internal_mapping_data);
- data.initialize(update_flags_mapping, quadrature, quadrature.size());
+ Assert(state == State::single_cell,
+ ExcMessage(
+ "This mapping info is not reinitialized for a single cell!"));
+ return 0;
}
else
{
- internal_mapping_data =
- mapping->get_data(update_flags_mapping, quadrature);
+ const unsigned int offset =
+ compute_geometry_index_offset(cell_index, face_number);
+ return data_index_offsets[offset];
}
+ }
+
- mapping->fill_fe_immersed_surface_values(cell,
- quadrature,
- *internal_mapping_data,
- mapping_data);
+
+ template <int dim, int spacedim, typename Number>
+ inline const DerivativeForm<1, dim, spacedim, Number> *
+ MappingInfo<dim, spacedim, Number>::get_jacobian(
+ const unsigned int offset) const
+ {
+ return &jacobians[offset];
}
- template <int dim, int spacedim>
- void
- MappingInfo<dim, spacedim>::compute_mapping_data_for_face_quadrature(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
- MappingData & mapping_data)
+ template <int dim, int spacedim, typename Number>
+ inline const DerivativeForm<1, spacedim, dim, Number> *
+ MappingInfo<dim, spacedim, Number>::get_inverse_jacobian(
+ const unsigned int offset) const
{
- update_flags_mapping |=
- mapping->requires_update_flags(update_flags_mapping);
+ return &inverse_jacobians[offset];
+ }
- mapping_data.initialize(quadrature.size(), update_flags_mapping);
- // reuse internal_mapping_data for MappingQ to avoid memory allocations
- if (const MappingQ<dim, spacedim> *mapping_q =
- dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
- {
- auto &data =
- dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
- *internal_mapping_data);
- data.initialize_face(update_flags_mapping,
- QProjector<dim>::project_to_oriented_face(
- ReferenceCells::get_hypercube<dim>(),
- quadrature,
- face_no,
- cell->face_orientation(face_no),
- cell->face_flip(face_no),
- cell->face_rotation(face_no)),
- quadrature.size());
-
- mapping_q->fill_mapping_data_for_face_quadrature(
- cell, face_no, quadrature, *internal_mapping_data, mapping_data);
- }
- else
- {
- auto internal_mapping_data =
- mapping->get_face_data(update_flags_mapping,
- hp::QCollection<dim - 1>(quadrature));
-
- mapping->fill_fe_face_values(cell,
- face_no,
- hp::QCollection<dim - 1>(quadrature),
- *internal_mapping_data,
- mapping_data);
- }
+ template <int dim, int spacedim, typename Number>
+ inline const Tensor<1, spacedim, Number> *
+ MappingInfo<dim, spacedim, Number>::get_normal_vector(
+ const unsigned int offset) const
+ {
+ return &normal_vectors[offset];
+ }
+
+
+
+ template <int dim, int spacedim, typename Number>
+ inline const Number *
+ MappingInfo<dim, spacedim, Number>::get_JxW(const unsigned int offset) const
+ {
+ return &JxW_values[offset];
+ }
+
+
+
+ template <int dim, int spacedim, typename Number>
+ const Mapping<dim, spacedim> &
+ MappingInfo<dim, spacedim, Number>::get_mapping() const
+ {
+ return *mapping;
+ }
+
+
+
+ template <int dim, int spacedim, typename Number>
+ UpdateFlags
+ MappingInfo<dim, spacedim, Number>::get_update_flags() const
+ {
+ return update_flags;
+ }
+
+
+
+ template <int dim, int spacedim, typename Number>
+ UpdateFlags
+ MappingInfo<dim, spacedim, Number>::get_update_flags_mapping() const
+ {
+ return update_flags_mapping;
+ }
+
+
+
+ template <int dim, int spacedim, typename Number>
+ boost::signals2::connection
+ MappingInfo<dim, spacedim, Number>::connect_is_reinitialized(
+ const std::function<void()> &set_is_reinitialized)
+ {
+ return is_reinitialized.connect(set_is_reinitialized);
}
} // namespace NonMatching
if (this->update_each &
(update_boundary_forms | update_normal_vectors | update_JxW_values))
{
- aux.resize(dim - 1,
- AlignedVector<Tensor<1, spacedim>>(n_original_q_points));
+ aux.resize(dim - 1);
+ aux[0].resize(n_original_q_points);
+ if (dim > 2)
+ aux[1].resize(n_original_q_points);
// Compute tangentials to the unit cell.
for (const unsigned int i : GeometryInfo<dim>::face_indices())
using namespace dealii;
+template <int dim>
void
test(const bool filtered_compression)
{
- constexpr unsigned int dim = 2;
constexpr unsigned int degree = 1;
FE_Q<dim> fe_q(degree);
initlog();
- test(true);
+ test<2>(true);
deallog << std::endl;
- test(false);
+ test<2>(false);
+ deallog << std::endl;
+ test<3>(true);
+ deallog << std::endl;
+ test<3>(false);
}
DEAL::check difference l2 norm cell: 0.00000
DEAL::check difference l2 norm surface: 0.00000
DEAL::check difference l2 norm faces: 0.00000
+DEAL::
+DEAL::check difference l2 norm cell: 0.00000
+DEAL::check difference l2 norm surface: 0.00000
+DEAL::check difference l2 norm faces: 0.00000
+DEAL::
+DEAL::check difference l2 norm cell: 0.00000
+DEAL::check difference l2 norm surface: 0.00000
+DEAL::check difference l2 norm faces: 0.00000
projects / dealii.git / commitdiff