return;
}
- internal::MatrixFreeFunctions::
- MappingInfoStorage<dim, dim, VectorizedArrayType>
- temp_data;
- temp_data.descriptor.resize(1);
- temp_data.descriptor[0].n_q_points = n_q_points;
- FEEvaluationData<dim, VectorizedArrayType, false> eval(
- std::make_tuple(&data.shape_info, nullptr, &temp_data, 0, 0));
+ FEEvaluationData<dim, VectorizedArrayType, false> eval(data.shape_info);
// prepare arrays
if (evaluation_flag != EvaluationFlags::nothing)
dim,
n_components,
evaluation_flag,
- &eval.get_orientation_map()(face_orientation, 0),
+ &eval.get_shape_info().face_orientations_quad(face_orientation, 0),
false,
data.n_q_points_face,
temp,
dim,
n_components,
integration_flag,
- &eval.get_orientation_map()(face_orientation, 0),
+ &eval.get_shape_info().face_orientations_quad(face_orientation, 0),
true,
data.n_q_points_face,
temp,
if (shape_data.nodal_at_cell_boundaries &&
eval.get_all_face_orientations()[v] != 0)
- orientation[v] =
- &eval.get_shape_info()
- .face_orientations[eval.get_all_face_orientations()[v]][0];
+ orientation[v] = &eval.get_shape_info().face_orientations_dofs(
+ eval.get_all_face_orientations()[v], 0);
}
else if (eval.get_face_orientation() != 0)
- orientation[0] =
- &eval.get_orientation_map()(eval.get_face_orientation(), 0);
+ orientation[0] = &eval.get_shape_info().face_orientations_dofs(
+ eval.get_face_orientation(), 0);
// face_to_cell_index_hermite
std::array<const unsigned int *, n_face_orientations> index_array_hermite =
const bool vectorization_possible =
all_faces_are_same && (sm_ptr == nullptr);
- std::array<Number2_ *, n_lanes> vector_ptrs;
+ std::array<Number2_ *, n_lanes> vector_ptrs;
+ std::array<unsigned int, n_lanes> reordered_indices;
if (vectorization_possible == false)
{
Assert(false, ExcNotImplemented());
}
}
+ else if (n_face_orientations == n_lanes)
+ {
+ for (unsigned int v = 0; v < n_lanes; ++v)
+ reordered_indices[v] =
+ dof_info.dof_indices_contiguous[dof_access_index]
+ [eval.get_cell_ids()[v]];
+ dof_indices = reordered_indices.data();
+ }
if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
{
n_components,
v,
evaluation_flag,
- &eval.get_orientation_map()
- [eval.get_all_face_orientations()[v]][0],
+ &eval.get_shape_info().face_orientations_quad(
+ eval.get_all_face_orientations()[v], 0),
false,
Utilities::pow(n_q_points_1d, dim - 1),
&temp[0][0],
dim,
n_components,
evaluation_flag,
- &eval.get_orientation_map()[eval.get_face_orientation()][0],
+ &eval.get_shape_info().face_orientations_quad(
+ eval.get_face_orientation(), 0),
false,
Utilities::pow(n_q_points_1d, dim - 1),
temp,
n_components,
v,
integration_flag,
- &eval.get_orientation_map()
- [eval.get_all_face_orientations()[v]][0],
+ &eval.get_shape_info().face_orientations_quad(
+ eval.get_all_face_orientations()[v], 0),
true,
Utilities::pow(n_q_points_1d, dim - 1),
&temp[0][0],
dim,
n_components,
integration_flag,
- &eval.get_orientation_map()[eval.get_face_orientation()][0],
+ &eval.get_shape_info().face_orientations_quad(
+ eval.get_face_orientation(), 0),
true,
Utilities::pow(n_q_points_1d, dim - 1),
temp,
int dim,
typename Number,
typename VectorizedArrayType>
- inline std::tuple<
- const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *,
- const internal::MatrixFreeFunctions::DoFInfo *,
- const internal::MatrixFreeFunctions::
- MappingInfoStorage<dim - is_face, dim, VectorizedArrayType> *,
- unsigned int,
- unsigned int>
- get_shape_info_and_indices(
- const MatrixFree<dim, Number, VectorizedArrayType> &data,
- const unsigned int dof_no,
- const unsigned int first_selected_component,
- const unsigned int quad_no,
- const unsigned int fe_degree,
- const unsigned int n_q_points,
- const unsigned int active_fe_index_given,
- const unsigned int active_quad_index_given)
+ inline typename FEEvaluationData<dim, VectorizedArrayType, is_face>::
+ InitializationData
+ get_shape_info_and_indices(
+ const MatrixFree<dim, Number, VectorizedArrayType> &data,
+ const unsigned int dof_no,
+ const unsigned int first_selected_component,
+ const unsigned int quad_no,
+ const unsigned int fe_degree,
+ const unsigned int n_q_points,
+ const unsigned int active_fe_index_given,
+ const unsigned int active_quad_index_given,
+ const unsigned int face_type)
{
- const auto &dof_info = data.get_dof_info(dof_no);
- const auto &mapping_storage = internal::MatrixFreeFunctions::
- MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
- data.get_mapping_info(), quad_no);
- const unsigned int active_fe_index =
+ typename FEEvaluationData<dim, VectorizedArrayType, is_face>::
+ InitializationData init_data;
+
+ init_data.dof_info = &data.get_dof_info(dof_no);
+ init_data.mapping_data =
+ &internal::MatrixFreeFunctions::
+ MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
+ data.get_mapping_info(), quad_no);
+
+ init_data.active_fe_index =
fe_degree != numbers::invalid_unsigned_int ?
- dof_info.fe_index_from_degree(first_selected_component, fe_degree) :
+ init_data.dof_info->fe_index_from_degree(first_selected_component,
+ fe_degree) :
(active_fe_index_given != numbers::invalid_unsigned_int ?
active_fe_index_given :
0);
- const unsigned int active_quad_index =
+ init_data.active_quad_index =
fe_degree == numbers::invalid_unsigned_int ?
(active_quad_index_given != numbers::invalid_unsigned_int ?
active_quad_index_given :
- std::min<unsigned int>(active_fe_index,
- mapping_storage.descriptor.size() - 1)) :
- mapping_storage.quad_index_from_n_q_points(n_q_points);
- return std::make_tuple(
- &data.get_shape_info(
- dof_no,
- quad_no,
- dof_info.component_to_base_index[first_selected_component],
- active_fe_index,
- active_quad_index),
- &dof_info,
- &mapping_storage,
- active_fe_index,
- active_quad_index);
+ std::min<unsigned int>(init_data.active_fe_index,
+ init_data.mapping_data->descriptor.size() -
+ 1)) :
+ init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
+
+ init_data.shape_info = &data.get_shape_info(
+ dof_no,
+ quad_no,
+ init_data.dof_info->component_to_base_index[first_selected_component],
+ init_data.active_fe_index,
+ init_data.active_quad_index);
+ init_data.descriptor =
+ &init_data.mapping_data->descriptor
+ [is_face ?
+ (init_data.active_quad_index * std::max<unsigned int>(1, dim - 1) +
+ (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
+ init_data.active_quad_index];
+
+ return init_data;
}
} // namespace internal
fe_degree,
n_q_points,
active_fe_index,
- active_quad_index),
+ active_quad_index,
+ face_type),
is_interior_face,
quad_no,
- face_type,
first_selected_component)
, scratch_data_array(data_in.acquire_scratch_data())
, matrix_info(&data_in)
static constexpr unsigned int n_lanes = sizeof(Number) / sizeof(ScalarNumber);
/**
- * Constructor, taking data stored in MatrixFree. If applied to problems
- * with more than one quadrature formula selected during construction of
- * `matrix_free`, `quad_no` allows to select the appropriate formula.
+ * Constructor, taking a single ShapeInfo object to inject the capability
+ * for evaluation and set up some data containers, compatible with the
+ * internal evaluation kernels of FEEvaluationImpl and friends. For actual
+ * use, select one of the derived classes FEEvaluation or FEFaceEvaluation
+ * with their respective arguments.
*/
- FEEvaluationData(const std::tuple<const ShapeInfoType *,
- const DoFInfo *,
- const MappingInfoStorageType *,
- unsigned int,
- unsigned int> &info,
- const bool is_interior_face = true,
- const unsigned int quad_no = 0,
- const unsigned int face_type = 0,
- const unsigned int first_selected_component = 0);
+ FEEvaluationData(const ShapeInfoType &info,
+ const bool is_interior_face = true);
/**
* Copy constructor.
/**
* If is_face is true, this function returns the orientation index within an
- * array of orientations as returned by get_orientation_map().
+ * array of orientations as stored in ShapeInfo for unknowns and quadrature
+ * points.
*
* @note This function depends on the internal representation of data, which
* is not stable between releases of deal.II, and is hence mostly for
unsigned int
get_face_orientation() const;
- /**
- * If is_face is true, this function returns the map re-ordering from face
- * quadrature points in standard orientation to a specific orientation.
- *
- * @note This function depends on the internal representation of data, which
- * is not stable between releases of deal.II, and is hence mostly for
- * internal use.
- */
- const Table<2, unsigned int> &
- get_orientation_map() const;
-
/**
* Return the current index in the access to compressed indices.
*
//@}
+ /**
+ * This data structure is used for the initialization by the derived
+ * FEEvaluation classes.
+ */
+ struct InitializationData
+ {
+ const ShapeInfoType * shape_info;
+ const DoFInfo * dof_info;
+ const MappingInfoStorageType *mapping_data;
+ unsigned int active_fe_index;
+ unsigned int active_quad_index;
+ const typename MappingInfoStorageType::QuadratureDescriptor *descriptor;
+ };
+
protected:
+ /**
+ * Constructor filling information available from a MatrixFree class,
+ * collected in an internal data structure to reduce overhead of
+ * initialization. Used in derived classes when this class is initialized
+ * from a MatrixFree object.
+ */
+ FEEvaluationData(const InitializationData &info,
+ const bool is_interior_face,
+ const unsigned int quad_no,
+ const unsigned int first_selected_component);
+
/**
* Constructor that comes with reduced functionality and works similarly as
* FEValues.
template <int dim, typename Number, bool is_face>
inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
- const std::tuple<const ShapeInfoType *,
- const DoFInfo *,
- const MappingInfoStorageType *,
- unsigned int,
- unsigned int> &shape_dof_info,
- const bool is_interior_face,
- const unsigned int quad_no,
- const unsigned int face_type,
- const unsigned int first_selected_component)
- : data(std::get<0>(shape_dof_info))
- , dof_info(std::get<1>(shape_dof_info))
- , mapping_data(std::get<2>(shape_dof_info))
+ const ShapeInfoType &shape_info,
+ const bool is_interior_face)
+ : FEEvaluationData(
+ InitializationData{&shape_info, nullptr, nullptr, 0, 0, nullptr},
+ is_interior_face,
+ 0,
+ 0)
+{}
+
+
+template <int dim, typename Number, bool is_face>
+inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
+ const InitializationData &initialization_data,
+ const bool is_interior_face,
+ const unsigned int quad_no,
+ const unsigned int first_selected_component)
+ : data(initialization_data.shape_info)
+ , dof_info(initialization_data.dof_info)
+ , mapping_data(initialization_data.mapping_data)
, quad_no(quad_no)
, n_fe_components(dof_info != nullptr ? dof_info->start_components.back() : 0)
, first_selected_component(first_selected_component)
- , active_fe_index(std::get<3>(shape_dof_info))
- , active_quad_index(std::get<4>(shape_dof_info))
- , descriptor(
- &mapping_data->descriptor
- [is_face ?
- (active_quad_index * std::max<unsigned int>(1, dim - 1) +
- (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
- active_quad_index])
- , n_quadrature_points(descriptor->n_q_points)
+ , active_fe_index(initialization_data.active_fe_index)
+ , active_quad_index(initialization_data.active_quad_index)
+ , descriptor(initialization_data.descriptor)
+ , n_quadrature_points(is_face ? data->n_q_points_face : data->n_q_points)
, jacobian(nullptr)
, J_value(nullptr)
, normal_vectors(nullptr)
, normal_x_jacobian(nullptr)
- , quadrature_weights(descriptor->quadrature_weights.begin())
+ , quadrature_weights(
+ descriptor != nullptr ? descriptor->quadrature_weights.begin() : nullptr)
+# ifdef DEBUG
, dof_values_initialized(false)
, values_quad_initialized(false)
, gradients_quad_initialized(false)
, hessians_quad_initialized(false)
, values_quad_submitted(false)
, gradients_quad_submitted(false)
+# endif
, cell(numbers::invalid_unsigned_int)
, is_interior_face(is_interior_face)
, dof_access_index(
, face_orientation(0)
, subface_index(0)
, cell_type(internal::MatrixFreeFunctions::general)
-{}
+{
+ // TODO: This currently fails for simplex/matrix_free_face_integral_01
+ // if (descriptor != nullptr)
+ // AssertDimension(n_quadrature_points, descriptor->n_q_points);
+}
normal_x_jacobian = nullptr;
quadrature_weights = other.quadrature_weights;
+# ifdef DEBUG
dof_values_initialized = false;
values_quad_initialized = false;
gradients_quad_initialized = false;
hessians_quad_initialized = false;
values_quad_submitted = false;
gradients_quad_submitted = false;
+# endif
cell = numbers::invalid_unsigned_int;
is_interior_face = other.is_interior_face;
"update_values|update_gradients"));
if (cell_type <= internal::MatrixFreeFunctions::affine)
{
- Assert(quadrature_weights != nullptr, ExcInternalError());
+ Assert(quadrature_weights != nullptr, ExcNotInitialized());
return J_value[0] * quadrature_weights[q_point];
}
else
-template <int dim, typename Number, bool is_face>
-inline const Table<2, unsigned int> &
-FEEvaluationData<dim, Number, is_face>::get_orientation_map() const
-{
- Assert(descriptor != nullptr, ExcNotInitialized());
- return descriptor->face_orientations;
-}
-
-
-
template <int dim, typename Number, bool is_face>
inline internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex
FEEvaluationData<dim, Number, is_face>::get_dof_access_index() const
if (flag == false)
{
- cell_data[my_q].descriptor[hpq].initialize(quad[my_q][hpq],
- update_default);
+ cell_data[my_q].descriptor[hpq].initialize(quad[my_q][hpq]);
const auto quad_face =
get_unique_face_quadratures(quad[my_q][hpq]);
if (quad_face.first.size() > 0) // line, quad
{
face_data[my_q].descriptor[hpq * scale].initialize(
- quad_face.first, update_default);
+ quad_face.first);
face_data_by_cells[my_q]
.descriptor[hpq * scale]
- .initialize(quad_face.first, update_default);
+ .initialize(quad_face.first);
}
if (quad_face.second.size() > 0) // triangle
AssertDimension(dim, 3);
face_data[my_q]
.descriptor[hpq * scale + (is_mixed_mesh ? 1 : 0)]
- .initialize(quad_face.second, update_default);
+ .initialize(quad_face.second);
face_data_by_cells[my_q]
.descriptor[hpq * scale + (is_mixed_mesh ? 1 : 0)]
- .initialize(quad_face.second, update_default);
+ .initialize(quad_face.second);
}
const auto face_quadrature_collection =
else
{
cell_data[my_q].descriptor[hpq].initialize(
- quad[my_q][hpq].get_tensor_basis()[0], update_default);
+ quad[my_q][hpq].get_tensor_basis()[0]);
const auto quad_face = quad[my_q][hpq].get_tensor_basis()[0];
- face_data[my_q].descriptor[hpq * scale].initialize(
- quad_face, this->update_flags_boundary_faces);
+ face_data[my_q].descriptor[hpq * scale].initialize(quad_face);
face_data[my_q].q_collection[hpq] =
dealii::hp::QCollection<dim - 1>(quad_face);
face_data_by_cells[my_q].descriptor[hpq * scale].initialize(
- quad_face, update_default);
+ quad_face);
reference_cell_types[my_q][hpq] =
dealii::ReferenceCells::get_hypercube<dim>();
}
shape_info.dofs_per_component_on_cell;
constexpr unsigned int hess_dim = dim * (dim + 1) / 2;
- MappingInfoStorage<dim, dim, VectorizedDouble> temp_data;
- temp_data.copy_descriptor(my_data);
- FEEvaluationData<dim, VectorizedDouble, false> eval(
- std::make_tuple(&shape_info, nullptr, &temp_data, 0, 0));
+ FEEvaluationData<dim, VectorizedDouble, false> eval(shape_info);
AlignedVector<VectorizedDouble> evaluation_data;
eval.set_data_pointers(&evaluation_data, dim);
shape_info.dofs_per_component_on_cell;
constexpr unsigned int hess_dim = dim * (dim + 1) / 2;
- MappingInfoStorage<dim - 1, dim, VectorizedDouble> temp_data;
- temp_data.copy_descriptor(my_data);
- FEEvaluationData<dim, VectorizedDouble, true> eval_int(
- std::make_tuple(&shape_info, nullptr, &temp_data, 0, 0), true);
- FEEvaluationData<dim, VectorizedDouble, true> eval_ext(
- std::make_tuple(&shape_info, nullptr, &temp_data, 0, 0), false);
+ FEEvaluationData<dim, VectorizedDouble, true> eval_int(shape_info,
+ true);
+ FEEvaluationData<dim, VectorizedDouble, true> eval_ext(shape_info,
+ false);
- AlignedVector<VectorizedDouble> evaluation_data, evaluation_data_ext;
+ // Let both evaluators use the same array as their use will not
+ // overlap
+ AlignedVector<VectorizedDouble> evaluation_data;
eval_int.set_data_pointers(&evaluation_data, dim);
- eval_ext.set_data_pointers(&evaluation_data_ext, dim);
+ eval_ext.set_data_pointers(&evaluation_data, dim);
for (unsigned int face = begin_face; face < end_face; ++face)
for (unsigned vv = 0; vv < n_lanes; vv += n_lanes_d)
{
internal::MatrixFreeFunctions::FaceToCellTopology<
VectorizedDouble::size()>
- face_double = {};
+ face_double = {};
face_double.interior_face_no = faces[face].interior_face_no;
face_double.exterior_face_no = faces[face].exterior_face_no;
face_double.face_orientation = faces[face].face_orientation;
*/
template <int dim_q>
void
- initialize(const Quadrature<dim_q> &quadrature,
- const UpdateFlags update_flags_inner_faces = update_default);
+ initialize(const Quadrature<dim_q> &quadrature);
/**
* Set up the lengths in the various members of this struct.
*/
void
- initialize(const Quadrature<1> &quadrature_1d,
- const UpdateFlags update_flags_inner_faces = update_default);
+ initialize(const Quadrature<1> &quadrature_1d);
/**
* Returns the memory consumption in bytes.
* when it is used in a vectorized context).
*/
AlignedVector<ScalarNumber> quadrature_weights;
-
- /**
- * For quadrature on faces, the evaluation of basis functions is not
- * in the correct order if a face is not in the standard orientation
- * to a given element. This data structure is used to re-order the
- * data evaluated on quadrature points to represent the correct order.
- */
- dealii::Table<2, unsigned int> face_orientations;
};
/**
unsigned int
quad_index_from_n_q_points(const unsigned int n_q_points) const;
- /**
- * Get the descriptor from another instance of QuadratureDescriptor,
- * clearing all other data fields.
- */
- template <typename Number2>
- void
- copy_descriptor(
- const MappingInfoStorage<structdim, spacedim, Number2> &other);
-
/**
* Helper function to determine which update flags must be set in the
* internal functions to initialize all data as requested by the user.
template <int dim_q>
void
MappingInfoStorage<structdim, spacedim, Number>::QuadratureDescriptor::
- initialize(const Quadrature<dim_q> &quadrature,
- const UpdateFlags update_flags_inner_faces)
+ initialize(const Quadrature<dim_q> &quadrature)
{
- Assert(structdim + 1 <= spacedim ||
- update_flags_inner_faces == update_default,
- ExcMessage("Volume cells do not allow for setting inner faces"));
this->quadrature = quadrature;
n_q_points = quadrature.size();
quadrature_weights.resize(n_q_points);
quadrature_weights[i] = quadrature.weight(i);
// note: quadrature_1d and tensor_quadrature_weights are not set up
-
- // TODO: set up face_orientations
- (void)update_flags_inner_faces;
}
template <int structdim, int spacedim, typename Number>
void
MappingInfoStorage<structdim, spacedim, Number>::QuadratureDescriptor::
- initialize(const Quadrature<1> &quadrature_1d,
- const UpdateFlags update_flags_inner_faces)
+ initialize(const Quadrature<1> &quadrature_1d)
{
- Assert(structdim + 1 <= spacedim ||
- update_flags_inner_faces == update_default,
- ExcMessage("Volume cells do not allow for setting inner faces"));
this->quadrature_1d = quadrature_1d;
quadrature = Quadrature<structdim>(quadrature_1d);
n_q_points = quadrature.size();
for (unsigned int i = 0; i < quadrature_1d.size(); ++i)
tensor_quadrature_weights[d][i] = quadrature_1d.weight(i);
}
-
- // face orientation for faces in 3D
- if (structdim == spacedim - 1 && spacedim == 3 &&
- update_flags_inner_faces != update_default)
- {
- const unsigned int n = quadrature_1d.size();
- face_orientations.reinit(8, n * n);
- for (unsigned int j = 0, i = 0; j < n; ++j)
- for (unsigned int k = 0; k < n; ++k, ++i)
- {
- // face_orientation=true, face_flip=false, face_rotation=false
- face_orientations[0][i] = i;
- // face_orientation=false, face_flip=false, face_rotation=false
- face_orientations[1][i] = j + k * n;
- // face_orientation=true, face_flip=true, face_rotation=false
- face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n;
- // face_orientation=false, face_flip=true, face_rotation=false
- face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n;
- // face_orientation=true, face_flip=false, face_rotation=true
- face_orientations[4][i] = j + (n - 1 - k) * n;
- // face_orientation=false, face_flip=false, face_rotation=true
- face_orientations[5][i] = k + (n - 1 - j) * n;
- // face_orientation=true, face_flip=true, face_rotation=true
- face_orientations[6][i] = (n - 1 - j) + k * n;
- // face_orientation=false, face_flip=true, face_rotation=true
- face_orientations[7][i] = (n - 1 - k) + j * n;
- }
- }
}
memory_consumption() const
{
std::size_t memory = sizeof(this) + quadrature.memory_consumption() +
- quadrature_weights.memory_consumption() +
- face_orientations.memory_consumption();
+ quadrature_weights.memory_consumption();
for (int d = 0; d < structdim; ++d)
memory += tensor_quadrature_weights[d].memory_consumption();
return memory;
- template <int structdim, int spacedim, typename Number>
- template <typename Number2>
- void
- MappingInfoStorage<structdim, spacedim, Number>::copy_descriptor(
- const MappingInfoStorage<structdim, spacedim, Number2> &other)
- {
- clear_data_fields();
- descriptor.clear();
- descriptor.resize(other.descriptor.size());
- for (unsigned int i = 0; i < descriptor.size(); ++i)
- {
- descriptor[i].n_q_points = other.descriptor[i].n_q_points;
- descriptor[i].quadrature_1d = other.descriptor[i].quadrature_1d;
- descriptor[i].quadrature = other.descriptor[i].quadrature;
- for (unsigned int d = 0; d < structdim; ++d)
- {
- descriptor[i].tensor_quadrature_weights[d].resize(
- other.descriptor[i].tensor_quadrature_weights[d].size());
- std::copy(
- other.descriptor[i].tensor_quadrature_weights[d].begin(),
- other.descriptor[i].tensor_quadrature_weights[d].end(),
- descriptor[i].tensor_quadrature_weights[d].begin());
- }
- descriptor[i].quadrature_weights.resize(
- other.descriptor[i].quadrature_weights.size());
- std::copy(other.descriptor[i].quadrature_weights.begin(),
- other.descriptor[i].quadrature_weights.end(),
- descriptor[i].quadrature_weights.begin());
- descriptor[i].face_orientations =
- other.descriptor[i].face_orientations;
- }
- }
-
-
-
template <int structdim, int spacedim, typename Number>
UpdateFlags
MappingInfoStorage<structdim, spacedim, Number>::compute_update_flags(
dealii::Table<2, unsigned int> face_to_cell_index_hermite;
/**
- * For degrees on faces, the basis functions are not
+ * For unknowns located on faces, the basis functions are not
* in the correct order if a face is not in the standard orientation
* to a given element. This data structure is used to re-order the
* basis functions to represent the correct order.
*/
- dealii::Table<2, unsigned int> face_orientations;
+ dealii::Table<2, unsigned int> face_orientations_dofs;
+
+ /**
+ * For interpretation of values at quadrature points, the order of
+ * points is not correct if a face is not in the standard orientation to
+ * a given element. This data structure is used to re-order the
+ * quadrature points to represent the correct order.
+ */
+ dealii::Table<2, unsigned int> face_orientations_quad;
private:
/**
// (similar to MappingInfoStorage::QuadratureDescriptor::initialize)
if (dim == 3)
{
- const unsigned int n = fe_degree + 1;
- face_orientations.reinit(8, n * n);
- for (unsigned int j = 0, i = 0; j < n; ++j)
- for (unsigned int k = 0; k < n; ++k, ++i)
- {
- // face_orientation=true, face_flip=false,
- // face_rotation=false
- face_orientations[0][i] = i;
- // face_orientation=false, face_flip=false,
- // face_rotation=false
- face_orientations[1][i] = j + k * n;
- // face_orientation=true, face_flip=true,
- // face_rotation=false
- face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n;
- // face_orientation=false, face_flip=true,
- // face_rotation=false
- face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n;
- // face_orientation=true, face_flip=false,
- // face_rotation=true
- face_orientations[4][i] = j + (n - 1 - k) * n;
- // face_orientation=false, face_flip=false,
- // face_rotation=true
- face_orientations[5][i] = k + (n - 1 - j) * n;
- // face_orientation=true, face_flip=true,
- // face_rotation=true
- face_orientations[6][i] = (n - 1 - j) + k * n;
- // face_orientation=false, face_flip=true,
- // face_rotation=true
- face_orientations[7][i] = (n - 1 - k) + j * n;
- }
+ auto compute_orientations =
+ [](const unsigned int n,
+ Table<2, unsigned int> &face_orientations) {
+ face_orientations.reinit(8, n * n);
+ for (unsigned int j = 0, i = 0; j < n; ++j)
+ for (unsigned int k = 0; k < n; ++k, ++i)
+ {
+ // face_orientation=true, face_flip=false,
+ // face_rotation=false
+ face_orientations[0][i] = i;
+ // face_orientation=false, face_flip=false,
+ // face_rotation=false
+ face_orientations[1][i] = j + k * n;
+ // face_orientation=true, face_flip=true,
+ // face_rotation=false
+ face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n;
+ // face_orientation=false, face_flip=true,
+ // face_rotation=false
+ face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n;
+ // face_orientation=true, face_flip=false,
+ // face_rotation=true
+ face_orientations[4][i] = j + (n - 1 - k) * n;
+ // face_orientation=false, face_flip=false,
+ // face_rotation=true
+ face_orientations[5][i] = k + (n - 1 - j) * n;
+ // face_orientation=true, face_flip=true,
+ // face_rotation=true
+ face_orientations[6][i] = (n - 1 - j) + k * n;
+ // face_orientation=false, face_flip=true,
+ // face_rotation=true
+ face_orientations[7][i] = (n - 1 - k) + j * n;
+ }
+ };
+ compute_orientations(fe_degree + 1, face_orientations_dofs);
+ compute_orientations(n_q_points_1d, face_orientations_quad);
}
else
{
- face_orientations.reinit(1, 1);
+ face_orientations_dofs.reinit(1, 1);
+ face_orientations_quad.reinit(1, 1);
}
}
projects / dealii.git / commitdiff