unsigned int padding_length;
const unsigned int mf_object_id;
- const unsigned int constraint_mask;
+ const dealii::internal::MatrixFreeFunctions::ConstraintTypes
+ constraint_mask;
const bool use_coloring;
# include <deal.II/lac/cuda_vector.h>
# include <deal.II/lac/la_parallel_vector.h>
+# include <deal.II/matrix_free/hanging_nodes_internal.h>
+
DEAL_II_NAMESPACE_OPEN
/**
* Mask deciding where constraints are set on a given cell.
*/
- unsigned int *constraint_mask;
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes *constraint_mask;
/**
* If true, use graph coloring has been used and we can simply add into
/**
* Mask deciding where constraints are set on a given cell.
*/
- std::vector<unsigned int *> constraint_mask;
+ std::vector<dealii::internal::MatrixFreeFunctions::ConstraintTypes *>
+ constraint_mask;
/**
* Grid dimensions associated to the different colors. The grid dimensions
/**
* Mask deciding where constraints are set on a given cell.
*/
- std::vector<unsigned int> constraint_mask;
+ std::vector<dealii::internal::MatrixFreeFunctions::ConstraintTypes>
+ constraint_mask;
/**
* If true, use graph coloring has been used and we can simply add into
# include <deal.II/fe/fe_dgq.h>
# include <deal.II/fe/fe_values.h>
-# include <deal.II/matrix_free/hanging_nodes_internal.h>
# include <deal.II/matrix_free/shape_info.h>
# include <cuda_runtime_api.h>
std::vector<Point<dim, Number>> q_points_host;
std::vector<Number> JxW_host;
std::vector<Number> inv_jacobian_host;
- std::vector<unsigned int> constraint_mask_host;
+ std::vector<dealii::internal::MatrixFreeFunctions::ConstraintTypes>
+ constraint_mask_host;
// Local buffer
std::vector<types::global_dof_index> local_dof_indices;
FEValues<dim> fe_values;
for (unsigned int i = 0; i < dofs_per_cell; ++i)
lexicographic_dof_indices[i] = local_dof_indices[lexicographic_inv[i]];
- const ArrayView<unsigned int> cell_id_view(constraint_mask_host[cell_id]);
+ const ArrayView<dealii::internal::MatrixFreeFunctions::ConstraintTypes>
+ cell_id_view(constraint_mask_host[cell_id]);
hanging_nodes.setup_constraints(cell,
partitioner,
n_cells * dim * dim * padding_length);
}
- alloc_and_copy(&data->constraint_mask[color],
- ArrayView<const unsigned int>(constraint_mask_host.data(),
- constraint_mask_host.size()),
- n_cells);
+ alloc_and_copy(
+ &data->constraint_mask[color],
+ ArrayView<const dealii::internal::MatrixFreeFunctions::ConstraintTypes>(
+ constraint_mask_host.data(), constraint_mask_host.size()),
+ n_cells);
}
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/matrix_free/face_info.h>
+#include <deal.II/matrix_free/hanging_nodes_internal.h>
#include <deal.II/matrix_free/mapping_info.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/task_info.h>
* Masks indicating for each cell and component if the optimized
* hanging-node constraint is applicable and if yes which type.
*/
- std::vector<unsigned int> hanging_node_constraint_masks;
+ std::vector<ConstraintTypes> hanging_node_constraint_masks;
/**
* This variable describes the position of constraints in terms of the
const TriaIterator<DoFCellAccessor<dim, dim, false>> &cell,
std::vector<types::global_dof_index> & dof_indices)
{
- const ArrayView<unsigned int> mask_view(
+ const ArrayView<ConstraintTypes> mask_view(
hanging_node_constraint_masks.data() +
cell_number * cell->get_fe().n_components(),
cell->get_fe().n_components());
{
template <int fe_degree, int n_q_points_1d>
static bool
- run(const unsigned int n_desired_components,
+ run(const unsigned int n_desired_components,
const FEEvaluationBaseData<dim,
typename Number::value_type,
is_face,
- Number> & fe_eval,
- const bool transpose,
- const std::array<unsigned int, Number::size()> &c_mask,
- Number * values)
+ Number> &fe_eval,
+ const bool transpose,
+ const std::array<MatrixFreeFunctions::ConstraintTypes, Number::size()>
+ & c_mask,
+ Number *values)
{
if (transpose)
run_internal<fe_degree, true>(n_desired_components,
template <int fe_degree, bool transpose>
static void
- run_internal(
- const unsigned int n_desired_components,
- const FEEvaluationBaseData<dim,
- typename Number::value_type,
- is_face,
- Number> & fe_eval,
- const std::array<unsigned int, Number::size()> &constraint_mask,
- Number * values)
+ run_internal(const unsigned int n_desired_components,
+ const FEEvaluationBaseData<dim,
+ typename Number::value_type,
+ is_face,
+ Number> &fe_eval,
+ const std::array<MatrixFreeFunctions::ConstraintTypes,
+ Number::size()> & constraint_mask,
+ Number * values)
{
const Number *weights = fe_eval.get_shape_info()
.data.front()
{
for (unsigned int v = 0; v < Number::size(); ++v)
{
- if (constraint_mask[v] == 0)
- continue;
-
const auto mask = constraint_mask[v];
+ if (mask == MatrixFreeFunctions::ConstraintTypes::unconstrained)
+ continue;
+
if (dim == 2) // 2D: only faces
{
// direction 0:
#include <deal.II/matrix_free/dof_info.h>
#include <deal.II/matrix_free/evaluation_flags.h>
+#include <deal.II/matrix_free/hanging_nodes_internal.h>
#include <deal.II/matrix_free/shape_info.h>
apply(const unsigned int n_components,
const unsigned int fe_degree,
const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType>
- & fe_eval,
- const bool transpose,
- const std::array<unsigned int, VectorizedArrayType::size()> &c_mask,
- VectorizedArrayType * values);
+ & fe_eval,
+ const bool transpose,
+ const std::array<MatrixFreeFunctions::ConstraintTypes,
+ VectorizedArrayType::size()> &c_mask,
+ VectorizedArrayType * values);
/**
* For faces.
apply(const unsigned int n_components,
const unsigned int fe_degree,
const FEEvaluationBaseData<dim, Number, true, VectorizedArrayType>
- & fe_eval,
- const bool transpose,
- const std::array<unsigned int, VectorizedArrayType::size()> &c_mask,
- VectorizedArrayType * values);
+ & fe_eval,
+ const bool transpose,
+ const std::array<MatrixFreeFunctions::ConstraintTypes,
+ VectorizedArrayType::size()> &c_mask,
+ VectorizedArrayType * values);
};
} // end of namespace internal
const unsigned int n_components,
const unsigned int fe_degree,
const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType>
- & fe_eval,
- const bool transpose,
- const std::array<unsigned int, VectorizedArrayType::size()> &c_mask,
- VectorizedArrayType * values)
+ & fe_eval,
+ const bool transpose,
+ const std::array<MatrixFreeFunctions::ConstraintTypes,
+ VectorizedArrayType::size()> &c_mask,
+ VectorizedArrayType * values)
{
instantiation_helper_run<
1,
const unsigned int n_components,
const unsigned int fe_degree,
const FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &fe_eval,
- const bool transpose,
- const std::array<unsigned int, VectorizedArrayType::size()> &c_mask,
- VectorizedArrayType * values)
+ const bool transpose,
+ const std::array<MatrixFreeFunctions::ConstraintTypes,
+ VectorizedArrayType::size()> &c_mask,
+ VectorizedArrayType * values)
{
instantiation_helper_run<
1,
this->dof_info
->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
- constexpr unsigned int n_lanes = VectorizedArrayType::size();
- std::array<unsigned int, n_lanes> constraint_mask;
+ constexpr unsigned int n_lanes = VectorizedArrayType::size();
+ std::array<internal::MatrixFreeFunctions::ConstraintTypes, n_lanes>
+ constraint_mask;
bool hn_available = false;
this->dof_info->hanging_node_constraint_masks[cell_dof_index];
constraint_mask[v] = mask;
- hn_available |= (mask != 0);
+ hn_available |=
+ (mask != internal::MatrixFreeFunctions::ConstraintTypes::unconstrained);
}
if (hn_available == false)
return; // no hanging node on cell batch -> nothing to do
for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
- constraint_mask[v] = 0;
+ constraint_mask[v] =
+ internal::MatrixFreeFunctions::ConstraintTypes::unconstrained;
internal::FEEvaluationHangingNodesFactory<dim, Number, VectorizedArrayType>::
apply(n_components,
{
namespace MatrixFreeFunctions
{
+ /**
+ * Here is the system for how we store constraint types in a binary mask.
+ * This is not a complete contradiction-free system, i.e., there are
+ * invalid states that we just assume that we never get.
+ *
+ * If the mask is zero, there are no constraints. Then, there are three
+ * different fields with one bit per dimension. The first field determines
+ * the type, or the position of an element along each direction. The
+ * second field determines if there is a constrained face with that
+ * direction as normal. The last field determines if there is a
+ * constrained edge of a given pair of coordinate planes, but where
+ * neither of the corresponding faces are constrained (only valid in 3D).
+ *
+ * The element is placed in the 'first position' along *-axis. These also
+ * determine which face is constrained. For example, in 2D, if
+ * face_x and type are set, then x = 0 is constrained.
+ */
+ enum ConstraintTypes : unsigned short
+ {
+ unconstrained = 0,
+
+ type_x = 1 << 0,
+ type_y = 1 << 1,
+ type_z = 1 << 2,
+
+ // Element has as a constraint at * = 0 or * = fe_degree face
+ face_x = 1 << 3,
+ face_y = 1 << 4,
+ face_z = 1 << 5,
+
+ // Element has as a constraint at * = 0 or * = fe_degree edge
+ edge_xy = 1 << 6,
+ edge_yz = 1 << 7,
+ edge_zx = 1 << 8
+ };
+
+
+
+ /**
+ * Global operator which returns an object in which all bits are set which
+ * are either set in the first or the second argument. This operator exists
+ * since if it did not then the result of the bit-or <tt>operator |</tt>
+ * would be an integer which would in turn trigger a compiler warning when
+ * we tried to assign it to an object of type UpdateFlags.
+ */
+ inline ConstraintTypes
+ operator|(const ConstraintTypes f1, const ConstraintTypes f2)
+ {
+ return static_cast<ConstraintTypes>(static_cast<unsigned short>(f1) |
+ static_cast<unsigned short>(f2));
+ }
+
+
+
+ /**
+ * Global operator which sets the bits from the second argument also in the
+ * first one.
+ */
+ inline ConstraintTypes &
+ operator|=(ConstraintTypes &f1, const ConstraintTypes f2)
+ {
+ f1 = f1 | f2;
+ return f1;
+ }
+
+
+
+#ifdef DEAL_II_COMPILER_CUDA_AWARE
+ __device__ inline ConstraintTypes
+ operator|(const ConstraintTypes f1, const ConstraintTypes f2)
+ {
+ return static_cast<ConstraintTypes>(static_cast<unsigned short>(f1) |
+ static_cast<unsigned short>(f2));
+ }
+
+
+
+ __device__ inline ConstraintTypes &
+ operator|=(ConstraintTypes &f1, const ConstraintTypes f2)
+ {
+ f1 = f1 | f2;
+ return f1;
+ }
+#endif
+
+
+
/**
* This class creates the mask used in the treatment of hanging nodes in
* CUDAWrappers::MatrixFree.
const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
const std::vector<unsigned int> & lexicographic_mapping,
std::vector<types::global_dof_index> &dof_indices,
- const ArrayView<unsigned int> & mask) const;
+ const ArrayView<ConstraintTypes> & mask) const;
private:
/**
line_to_cells;
};
- /**
- * Here is the system for how we store constraint types in a binary mask.
- * This is not a complete contradiction-free system, i.e., there are
- * invalid states that we just assume that we never get.
- *
- * If the mask is zero, there are no constraints. Then, there are three
- * different fields with one bit per dimension. The first field determines
- * the type, or the position of an element along each direction. The
- * second field determines if there is a constrained face with that
- * direction as normal. The last field determines if there is a
- * constrained edge of a given pair of coordinate planes, but where
- * neither of the corresponding faces are constrained (only valid in 3D).
- *
- * The element is placed in the 'first position' along *-axis. These also
- * determine which face is constrained. For example, in 2D, if
- * face_x and type are set, then x = 0 is constrained.
- */
- enum ConstraintTypes : unsigned int
- {
- unconstrained = 0,
-
- type_x = 1 << 0,
- type_y = 1 << 1,
- type_z = 1 << 2,
-
- // Element has as a constraint at * = 0 or * = fe_degree face
- face_x = 1 << 3,
- face_y = 1 << 4,
- face_z = 1 << 5,
-
- // Element has as a constraint at * = 0 or * = fe_degree edge
- edge_xy = 1 << 6,
- edge_yz = 1 << 7,
- edge_zx = 1 << 8
- };
template <int dim>
const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
const std::vector<unsigned int> & lexicographic_mapping,
std::vector<types::global_dof_index> &dof_indices,
- const ArrayView<unsigned int> & masks) const
+ const ArrayView<ConstraintTypes> & masks) const
{
bool cell_has_hanging_node_constraints = false;
++c, ++comp)
{
auto &mask = masks[comp];
- mask = 0;
+ mask = ConstraintTypes::unconstrained;
const auto &fe = cell->get_fe().base_element(base_element_index);
// For each line on cell, which faces does it belong to, what is
// the edge mask, what is the types of the faces it belong to,
// and what is the type along the edge.
- const unsigned int line_to_edge[12][4] = {
+ const ConstraintTypes line_to_edge[12][4] = {
{ConstraintTypes::face_x | ConstraintTypes::face_z,
ConstraintTypes::edge_zx,
ConstraintTypes::type_x | ConstraintTypes::type_z,
ConstraintTypes::type_y},
{ConstraintTypes::face_x | ConstraintTypes::face_z,
ConstraintTypes::edge_zx,
- 0,
+ ConstraintTypes::unconstrained,
ConstraintTypes::type_y},
{ConstraintTypes::face_y | ConstraintTypes::face_z,
ConstraintTypes::edge_yz,
ConstraintTypes::type_x},
{ConstraintTypes::face_y | ConstraintTypes::face_z,
ConstraintTypes::edge_yz,
- 0,
+ ConstraintTypes::unconstrained,
ConstraintTypes::type_x},
{ConstraintTypes::face_x | ConstraintTypes::face_y,
ConstraintTypes::edge_xy,
ConstraintTypes::type_z},
{ConstraintTypes::face_x | ConstraintTypes::face_y,
ConstraintTypes::edge_xy,
- 0,
+ ConstraintTypes::unconstrained,
ConstraintTypes::type_z}};
for (unsigned int local_line = 0;
bool transpose,
typename Number>
__device__ inline void
- interpolate_boundary_2d(const unsigned int constraint_mask, Number *values)
+ interpolate_boundary_2d(
+ const dealii::internal::MatrixFreeFunctions::ConstraintTypes
+ constraint_mask,
+ Number *values)
{
const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
const unsigned int y_idx = threadIdx.y;
(constraint_mask &
(((direction == 0) ?
dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y :
- 0) |
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::
+ unconstrained) |
((direction == 1) ?
dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x :
- 0)));
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::
+ unconstrained)));
// Flag is true if for the given direction, the dof is constrained with
// the right type and is on the correct side (left (= 0) or right (=
bool transpose,
typename Number>
__device__ inline void
- interpolate_boundary_3d(const unsigned int constraint_mask, Number *values)
+ interpolate_boundary_3d(
+ const dealii::internal::MatrixFreeFunctions::ConstraintTypes
+ constraint_mask,
+ Number *values)
{
const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
const unsigned int y_idx = threadIdx.y;
*/
template <int dim, int fe_degree, bool transpose, typename Number>
__device__ void
- resolve_hanging_nodes(const unsigned int constraint_mask, Number *values)
+ resolve_hanging_nodes(
+ const dealii::internal::MatrixFreeFunctions::ConstraintTypes
+ constraint_mask,
+ Number *values)
{
if (dim == 2)
{
AlignedVector<VectorizedArrayType> values_dofs(
dofs_per_component);
- std::array<unsigned int, VectorizedArrayType::size()>
+ std::array<dealii::internal::MatrixFreeFunctions::
+ ConstraintTypes,
+ VectorizedArrayType::size()>
constraint_mask;
constraint_mask[0] = mask;
new_dof_indices.reserve(dof_indices.size());
new_constraint_indicator.reserve(constraint_indicator.size());
- std::vector<unsigned int> new_hanging_node_constraint_masks;
+ std::vector<ConstraintTypes> new_hanging_node_constraint_masks;
new_hanging_node_constraint_masks.reserve(
new_hanging_node_constraint_masks.size());
.second = new_constraint_indicator.size();
if (hanging_node_constraint_masks.size() > 0)
- new_hanging_node_constraint_masks.push_back(0);
+ new_hanging_node_constraint_masks.push_back(
+ ConstraintTypes::unconstrained);
}
position_cell += n_vect;
}
{
template <int fe_degree, int n_q_points_1d>
static bool
- run(const FEEvaluationBaseData<dim,
- typename Number::value_type,
- is_face,
- Number> & fe_eval,
- const bool transpose,
- const std::array<unsigned int, Number::size()> &c_mask,
- Number * values)
+ run(
+ const FEEvaluationBaseData<dim,
+ typename Number::value_type,
+ is_face,
+ Number> &fe_eval,
+ const bool transpose,
+ const std::array<dealii::internal::MatrixFreeFunctions::ConstraintTypes,
+ Number::size()> & c_mask,
+ Number * values)
{
Assert(is_face == false, ExcInternalError());
template <int fe_degree_, unsigned int direction, bool transpose>
static void
- run_2D(const FEEvaluationBaseData<dim,
- typename Number::value_type,
- is_face,
- Number> & fe_eval,
- const std::array<unsigned int, Number::size()> &constraint_mask,
- Number * values)
+ run_2D(
+ const FEEvaluationBaseData<dim,
+ typename Number::value_type,
+ is_face,
+ Number> &fe_eval,
+ const std::array<dealii::internal::MatrixFreeFunctions::ConstraintTypes,
+ Number::size()> & constraint_mask,
+ Number * values)
{
const auto &constraint_weights =
fe_eval.get_shape_info().data.front().subface_interpolation_matrix;
template <int fe_degree_, unsigned int direction, bool transpose>
static void
- run_3D(const FEEvaluationBaseData<dim,
- typename Number::value_type,
- is_face,
- Number> & fe_eval,
- const std::array<unsigned int, Number::size()> &constraint_mask,
- Number * values)
+ run_3D(
+ const FEEvaluationBaseData<dim,
+ typename Number::value_type,
+ is_face,
+ Number> &fe_eval,
+ const std::array<dealii::internal::MatrixFreeFunctions::ConstraintTypes,
+ Number::size()> & constraint_mask,
+ Number * values)
{
const auto &constraint_weights =
fe_eval.get_shape_info().data.front().subface_interpolation_matrix;
template <int dim>
void
-test(const unsigned int degree, const unsigned int mask_value)
+test(const unsigned int degree,
+ const dealii::internal::MatrixFreeFunctions::ConstraintTypes mask_value)
{
Triangulation<dim> tria;
GridGenerator::subdivided_hyper_cube(tria, 2);
FEEvaluation<dim, -1, 0, 1, double> eval(matrix_free);
eval.reinit(0);
- std::array<unsigned int, VectorizedArray<double>::size()> cmask;
- std::fill(cmask.begin(), cmask.end(), 0);
+ std::array<dealii::internal::MatrixFreeFunctions::ConstraintTypes,
+ VectorizedArray<double>::size()>
+ cmask;
+ std::fill(
+ cmask.begin(),
+ cmask.end(),
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::unconstrained);
cmask[0] = mask_value;
for (unsigned int b = 0; b < 2; ++b)
for (unsigned int degree = 1; degree <= 3; ++degree)
{
- test<2>(degree, 0);
+ test<2>(
+ degree,
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::unconstrained);
deallog << std::endl;
test<2>(degree,
- 0 |
- dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x |
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x |
dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_x |
dealii::internal::MatrixFreeFunctions::ConstraintTypes::
type_y); // face 0/0
test<2>(degree,
- 0 |
- dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x |
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x |
dealii::internal::MatrixFreeFunctions::ConstraintTypes::
type_x); // face 0/1
test<2>(degree,
- 0 |
- dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x |
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x |
dealii::internal::MatrixFreeFunctions::ConstraintTypes::
type_y); // face 1/0
- test<2>(degree,
- 0 | dealii::internal::MatrixFreeFunctions::ConstraintTypes::
- face_x); // face 1/1
+ test<2>(
+ degree,
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x); // face
+ // 1/1
deallog << std::endl;
test<2>(degree,
- 0 |
- dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y |
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y |
dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_y |
dealii::internal::MatrixFreeFunctions::ConstraintTypes::
type_x); // face 2/0
test<2>(degree,
- 0 |
- dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y |
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y |
dealii::internal::MatrixFreeFunctions::ConstraintTypes::
type_y); // face 2/1
test<2>(degree,
- 0 |
- dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y |
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y |
dealii::internal::MatrixFreeFunctions::ConstraintTypes::
type_x); // face 3/0
- test<2>(degree,
- 0 | dealii::internal::MatrixFreeFunctions::ConstraintTypes::
- face_y); // face 3/1
+ test<2>(
+ degree,
+ dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y); // face
+ // 3/1
deallog << std::endl;
}
projects / dealii.git / commitdiff