dim,
n_points_1d,
n_points_1d,
- VectorizedArrayType>
- eval(AlignedVector<VectorizedArrayType>(),
+ VectorizedArrayType,
+ Number>
+ eval({},
data.get_shape_info().data[0].shape_gradients_collocation_eo,
- AlignedVector<VectorizedArrayType>());
+ {});
AlignedVector<VectorizedArrayType> buffer(phi.static_n_q_points *
phi.n_components);
* we need to store the evaluations of the 1d polynomials at
* the 1d quadrature points. That is what this variable is for.
*/
- internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> shape_info;
+ internal::MatrixFreeFunctions::ShapeInfo<double> shape_info;
/**
* In case the quadrature rule given represents a tensor product
std::vector<typename Number::value_type> constraint_pool_data;
std::vector<unsigned int> constraint_pool_row_index;
- std::vector<ShapeInfo<Number>> shape_infos;
+ std::vector<ShapeInfo<typename Number::value_type>> shape_infos;
std::vector<compressed_constraint_kind> hanging_node_constraint_masks;
std::vector<unsigned int> active_fe_indices;
// might have non-symmetric quadrature formula, so use the more
// conservative 'evaluate_general' scheme rather than 'even_odd' as the
// Hessians are not used very often
- const MatrixFreeFunctions::UnivariateShapeData<Number> &data =
+ const MatrixFreeFunctions::UnivariateShapeData<Number2> &data =
fe_eval.get_shape_info().data[0];
AssertDimension(data.shape_gradients_collocation.size(),
data.n_q_points_1d * data.n_q_points_1d);
using Number2 =
typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
- const MatrixFreeFunctions::UnivariateShapeData<Number> &data =
+ const MatrixFreeFunctions::UnivariateShapeData<Number2> &data =
fe_eval.get_shape_info().data[0];
AssertDimension(data.shape_gradients_collocation.size(),
data.n_q_points_1d * data.n_q_points_1d);
const auto face_no = fe_eval.get_face_no(v);
grad_weight[v] =
- shape_data.shape_data_on_face[0][fe_degree +
- (integrate ?
- (2 - (face_no % 2)) :
- (1 + (face_no % 2)))][0];
+ shape_data
+ .shape_data_on_face[0][fe_degree + (integrate ?
+ (2 - (face_no % 2)) :
+ (1 + (face_no % 2)))];
index_array_hermite[v] =
&fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
Number>
{
private:
- template <int structdim, unsigned int direction, bool transpose>
+ template <int structdim,
+ unsigned int direction,
+ bool transpose,
+ typename Number2>
static void
- interpolate(const unsigned int offset,
- const unsigned int outer_stride,
- const unsigned int given_degree,
- const Number mask_weight,
- const Number mask_write,
- const Number *DEAL_II_RESTRICT weights,
- Number *DEAL_II_RESTRICT values)
+ interpolate(const unsigned int offset,
+ const unsigned int outer_stride,
+ const unsigned int given_degree,
+ const Number mask_weight,
+ const Number mask_write,
+ const Number2 *DEAL_II_RESTRICT weights,
+ Number *DEAL_II_RESTRICT values)
{
static constexpr unsigned int max_n_points_1D = 40;
}
public:
- template <bool transpose>
+ template <bool transpose, typename Number2>
static void
run_internal(
- const unsigned int n_components,
- const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const unsigned int n_components,
+ const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
const std::array<MatrixFreeFunctions::compressed_constraint_kind,
- Number::size()> &constraint_mask,
- Number *values)
+ Number::size()> &constraint_mask,
+ Number *values)
{
const unsigned int given_degree =
fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
- const Number *DEAL_II_RESTRICT weights =
+ const Number2 *DEAL_II_RESTRICT weights =
shape_info.data.front().subface_interpolation_matrices[0].data();
const unsigned int points = given_degree + 1;
}
public:
- template <bool transpose>
+ template <bool transpose, typename Number2>
static void
run_internal(
- const unsigned int n_desired_components,
- const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const unsigned int n_desired_components,
+ const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
const std::array<MatrixFreeFunctions::compressed_constraint_kind,
- Number::size()> &constraint_mask,
- Number *values)
+ Number::size()> &constraint_mask,
+ Number *values)
{
const unsigned int given_degree =
fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
struct FEEvaluationImplHangingNodes
{
public:
- template <int fe_degree>
+ template <int fe_degree, typename Number2>
static bool
- run(const unsigned int n_desired_components,
- const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
- const bool transpose,
+ run(const unsigned int n_desired_components,
+ const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
+ const bool transpose,
const std::array<MatrixFreeFunctions::compressed_constraint_kind,
- Number::size()> &c_mask,
- Number *values)
+ Number::size()> &c_mask,
+ Number *values)
{
using RunnerType =
FEEvaluationImplHangingNodesRunner<used_runner_type<fe_degree>(),
struct FEEvaluationHangingNodesFactory
{
static void
- apply(const unsigned int n_components,
- const unsigned int fe_degree,
- const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &shape_info,
- const bool transpose,
+ apply(const unsigned int n_components,
+ const unsigned int fe_degree,
+ const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const bool transpose,
const std::array<MatrixFreeFunctions::compressed_constraint_kind,
- VectorizedArrayType::size()> &c_mask,
- VectorizedArrayType *values);
+ VectorizedArrayType::size()> &c_mask,
+ VectorizedArrayType *values);
};
} // end of namespace internal
template <int dim, typename Number, typename VectorizedArrayType>
void
FEEvaluationHangingNodesFactory<dim, Number, VectorizedArrayType>::apply(
- const unsigned int n_components,
- const unsigned int fe_degree,
- const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &shape_info,
- const bool transpose,
+ const unsigned int n_components,
+ const unsigned int fe_degree,
+ const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const bool transpose,
const std::array<MatrixFreeFunctions::compressed_constraint_kind,
- VectorizedArrayType::size()> &c_mask,
- VectorizedArrayType *values)
+ VectorizedArrayType::size()> &c_mask,
+ VectorizedArrayType *values)
{
instantiation_helper_degree_run<
1,
(void)base_element_number;
Assert(this->data == nullptr, ExcInternalError());
- this->data =
- new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
- Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
- fe,
- fe.component_to_base_index(first_selected_component).first);
+ this->data = new internal::MatrixFreeFunctions::ShapeInfo<Number>(
+ Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
+ fe,
+ fe.component_to_base_index(first_selected_component).first);
this->set_data_pointers(scratch_data_array, n_components_);
}
{
Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
this->data =
- new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
- *other.data);
+ new internal::MatrixFreeFunctions::ShapeInfo<Number>(*other.data);
// Create deep copy of mapped geometry for use in parallel
this->mapped_geometry =
{
Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
this->data =
- new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
- *other.data);
+ new internal::MatrixFreeFunctions::ShapeInfo<Number>(*other.data);
scratch_data_array = new AlignedVector<VectorizedArrayType>();
// Create deep copy of mapped geometry for use in parallel
template <int dim, typename Number, bool is_face>
class FEEvaluationData
{
- using ShapeInfoType = internal::MatrixFreeFunctions::ShapeInfo<Number>;
+ using ShapeInfoType = internal::MatrixFreeFunctions::ShapeInfo<
+ typename internal::VectorizedArrayTrait<Number>::value_type>;
using MappingInfoStorageType = internal::MatrixFreeFunctions::
MappingInfoStorage<(is_face ? dim - 1 : dim), dim, Number>;
using DoFInfo = internal::MatrixFreeFunctions::DoFInfo;
public:
static constexpr unsigned int dimension = dim;
- using NumberType = Number;
- using shape_info_number_type = Number;
+ using NumberType = Number;
using ScalarNumber =
typename internal::VectorizedArrayTrait<Number>::value_type;
+ using shape_info_number_type = ScalarNumber;
static constexpr unsigned int n_lanes =
internal::VectorizedArrayTrait<Number>::width();
typename VectorizedDouble>
void
mapping_q_compute_range(
- const unsigned int begin_cell,
- const unsigned int end_cell,
- const std::vector<GeometryType> &cell_type,
- const std::vector<bool> &process_cell,
- const UpdateFlags update_flags_cells,
- const AlignedVector<double> &plain_quadrature_points,
- const ShapeInfo<VectorizedDouble> &shape_info,
+ const unsigned int begin_cell,
+ const unsigned int end_cell,
+ const std::vector<GeometryType> &cell_type,
+ const std::vector<bool> &process_cell,
+ const UpdateFlags update_flags_cells,
+ const AlignedVector<double> &plain_quadrature_points,
+ const ShapeInfo<double> &shape_info,
MappingInfoStorage<dim, dim, VectorizedArrayType> &my_data)
{
constexpr unsigned int n_lanes = VectorizedArrayType::size();
const unsigned int begin_face,
const unsigned int end_face,
const std::vector<FaceToCellTopology<VectorizedArrayType::size()>>
- &faces,
- const std::vector<GeometryType> &face_type,
- const std::vector<bool> &process_face,
- const UpdateFlags update_flags_faces,
- const AlignedVector<double> &plain_quadrature_points,
- const ShapeInfo<VectorizedDouble> &shape_info,
+ &faces,
+ const std::vector<GeometryType> &face_type,
+ const std::vector<bool> &process_face,
+ const UpdateFlags update_flags_faces,
+ const AlignedVector<double> &plain_quadrature_points,
+ const ShapeInfo<double> &shape_info,
MappingInfoStorage<dim - 1, dim, VectorizedArrayType> &my_data)
{
constexpr unsigned int n_lanes = VectorizedArrayType::size();
// functions or the quadrature points; shape info is merely a vehicle to
// return us the right interpolation matrices from the cell support
// points to the cell and face quadrature points.
- std::vector<ShapeInfo<VectorizedDouble>> shape_infos(cell_data.size());
+ std::vector<ShapeInfo<double>> shape_infos(cell_data.size());
{
FE_DGQ<dim> fe_geometry(mapping_degree);
for (unsigned int my_q = 0; my_q < cell_data.size(); ++my_q)
/**
* Return the unit cell information for given hp-index.
*/
- const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> &
get_shape_info(const unsigned int dof_handler_index_component = 0,
const unsigned int quad_index = 0,
const unsigned int fe_base_element = 0,
/**
* Contains shape value information on the unit cell.
*/
- Table<4, internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>>
- shape_info;
+ Table<4, internal::MatrixFreeFunctions::ShapeInfo<Number>> shape_info;
/**
* Describes how the cells are gone through. With the cell level (first
template <int dim, typename Number, typename VectorizedArrayType>
-inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
+inline const internal::MatrixFreeFunctions::ShapeInfo<Number> &
MatrixFree<dim, Number, VectorizedArrayType>::get_shape_info(
const unsigned int dof_handler_index,
const unsigned int index_quad,
zero_vector_region(const unsigned int range_index, VectorType &vec) const
{
if (range_index == numbers::invalid_unsigned_int || range_index == 0)
- vec = typename VectorType::value_type();
+ {
+ if constexpr (std::is_same_v<
+ ArrayView<typename VectorType::value_type>,
+ VectorType>)
+ {
+ for (unsigned int i = 0; i < vec.size(); ++i)
+ vec[i] = typename VectorType::value_type();
+ }
+ else
+ vec = typename VectorType::value_type();
+ }
}
};
- const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
- *shape_info_base;
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> *shape_info_base;
void
reinit(unsigned int cell_batch_index,
- const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
- *shape_info = nullptr);
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> *shape_info =
+ nullptr);
void
reinit(unsigned int cell_batch_index,
unsigned int face_number,
- const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
- *shape_info = nullptr);
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> *shape_info =
+ nullptr);
};
template <int dim,
Number,
VectorizedArrayType>::
reinit(const unsigned int cell_batch_index,
- const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
- *shape_info)
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> *shape_info)
{
Assert(this->mapped_geometry == nullptr,
ExcMessage(
VectorizedArrayType>::
reinit(const unsigned int cell_batch_index,
const unsigned int face_number,
- const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
- *shape_info)
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> *shape_info)
{
Assert(
this->quad_no <
matrix_free.initialize_dof_vector(dst);
- std::array<internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>,
- 2 * dim>
+ std::array<internal::MatrixFreeFunctions::ShapeInfo<Number>, 2 * dim>
shape_info_shift;
const QGauss<1> quadrature_1D(n_points);
const Quadrature<1> quadrature_1D_shift_plus =
shift_1d_quadrature(quadrature_1D, +1.);
- dealii::internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
- shape_info_base(quadrature_1D, fe);
- dealii::internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
+ dealii::internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info_base(
+ quadrature_1D, fe);
+ dealii::internal::MatrixFreeFunctions::ShapeInfo<Number>
shape_info_shift_plus(quadrature_1D_shift_plus, fe);
- dealii::internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
+ dealii::internal::MatrixFreeFunctions::ShapeInfo<Number>
shape_info_shift_minus(quadrature_1D_shift_minus, fe);
for (unsigned int f = 0; f < 2 * dim; ++f)
const auto w =
transpose ?
constraint_weights[i * (fe_degree + 1) +
- interp_idx][v] :
+ interp_idx] :
constraint_weights[interp_idx *
(fe_degree + 1) +
- i][v];
+ i];
t += w * values_temp[real_idx][v];
}
}
transpose ?
constraint_weights[(fe_degree - i) *
(fe_degree + 1) +
- fe_degree - interp_idx]
- [v] :
+ fe_degree - interp_idx] :
constraint_weights[(fe_degree - interp_idx) *
(fe_degree + 1) +
- fe_degree - i][v];
+ fe_degree - i];
t += w * values_temp[real_idx][v];
}
}
const auto w =
transpose ?
constraint_weights[i * (fe_degree + 1) +
- interp_idx][v] :
+ interp_idx] :
constraint_weights[interp_idx *
(fe_degree + 1) +
- i][v];
+ i];
t += w * values_temp[real_idx][v];
}
}
transpose ?
constraint_weights[(fe_degree - i) *
(fe_degree + 1) +
- fe_degree - interp_idx]
- [v] :
+ fe_degree - interp_idx] :
constraint_weights[(fe_degree -
interp_idx) *
(fe_degree + 1) +
- fe_degree - i][v];
+ fe_degree - i];
t += w * values_temp[real_idx][v];
}
}
dim,
n_points_1d,
n_points_1d,
- VectorizedArrayType>
- eval(AlignedVector<VectorizedArrayType>(),
+ VectorizedArrayType,
+ Number>
+ eval({},
data.get_shape_info().data[0].shape_gradients_collocation_eo,
- AlignedVector<VectorizedArrayType>());
+ {});
internal::EvaluatorTensorProduct<
internal::EvaluatorVariant::evaluate_evenodd,
dim - 1,
n_q_points_1d,
n_q_points_1d,
- VectorizedArrayType>
+ VectorizedArrayType,
+ Number>
eval_face({},
data.get_shape_info().data[0].shape_gradients_collocation_eo,
{});
phi_m.reinit(cell, face);
- const AlignedVector<VectorizedArrayType> &shape_data =
+ const AlignedVector<Number> &shape_data =
data.get_shape_info().data.front().quadrature_data_on_face[face %
2];
const std::array<int, 2> n_blocks{