const bool evaluate_gradients,
const unsigned int active_fe_index,
const unsigned int first_selected_component,
- const unsigned int cell,
- const std::array<unsigned int, n_face_orientations> face_no,
+ const std::array<unsigned int, n_face_orientations> cells,
+ const std::array<unsigned int, n_face_orientations> face_nos,
const unsigned int subface_index,
const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
- const std::array<unsigned int, n_face_orientations> face_orientation,
+ const std::array<unsigned int, n_face_orientations> face_orientations,
const Table<2, unsigned int> & orientation_map)
{
return process_and_io( //
evaluate_gradients,
active_fe_index,
first_selected_component,
- cell,
- face_no,
+ cells,
+ face_nos,
subface_index,
dof_access_index,
- face_orientation,
+ face_orientations,
orientation_map,
[](auto & temp_1,
auto & temp_2,
const bool integrate_gradients,
const unsigned int active_fe_index,
const unsigned int first_selected_component,
- const unsigned int cell,
- const std::array<unsigned int, n_face_orientations> face_no,
+ const std::array<unsigned int, n_face_orientations> cells,
+ const std::array<unsigned int, n_face_orientations> face_nos,
const unsigned int subface_index,
const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
- const std::array<unsigned int, n_face_orientations> face_orientation,
+ const std::array<unsigned int, n_face_orientations> face_orientations,
const Table<2, unsigned int> & orientation_map)
{
return process_and_io( //
integrate_gradients,
active_fe_index,
first_selected_component,
- cell,
- face_no,
+ cells,
+ face_nos,
subface_index,
dof_access_index,
- face_orientation,
+ face_orientations,
orientation_map,
[](const auto &temp_1,
const auto &temp_2,
// case 5: default vector access, must be handled separately, just do
// the face-normal interpolation
- AssertDimension(face_no.size(), 1);
+ AssertDimension(face_nos.size(), 1);
FEFaceNormalEvaluationImpl<dim,
fe_degree,
n_components,
VectorizedArrayType>::
template interpolate<false, false>(
- data, temp1, values_array, integrate_gradients, face_no[0]);
+ data, temp1, values_array, integrate_gradients, face_nos[0]);
},
[&](auto &temp1, const auto &dofs_per_face) {
if (fe_degree > -1 &&
const bool do_gradients,
const unsigned int active_fe_index,
const unsigned int first_selected_component,
- const unsigned int cell,
- const std::array<unsigned int, n_face_orientations> face_no,
+ const std::array<unsigned int, n_face_orientations> cells,
+ const std::array<unsigned int, n_face_orientations> face_nos,
const unsigned int subface_index,
const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
- const std::array<unsigned int, n_face_orientations> face_orientation,
+ const std::array<unsigned int, n_face_orientations> face_orientations,
const Table<2, unsigned int> & orientation_map,
const Function1a & function_1a,
const Function1b & function_1b,
const Function5 & function_5,
const Function0 & function_0)
{
- (void)subface_index;
+ const unsigned int cell = cells[0];
if (integrate &&
- (face_orientation[0] > 0 &&
+ (face_orientations[0] > 0 &&
subface_index < GeometryInfo<dim>::max_children_per_cell))
{
- AssertDimension(face_orientation.size(), 1);
- adjust_for_face_orientation(face_orientation[0],
+ AssertDimension(face_orientations.size(), 1);
+ adjust_for_face_orientation(face_orientations[0],
orientation_map,
true,
do_values,
data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
data.data.front()
.shape_data_on_face[0][fe_degree + (integrate ?
- (2 - (face_no[0] % 2)) :
- (1 + (face_no[0] % 2)))] :
+ (2 - (face_nos[0] % 2)) :
+ (1 + (face_nos[0] % 2)))] :
VectorizedArrayType(0.0 /*dummy*/);
constexpr unsigned int static_dofs_per_component =
// re-orientation
std::array<const unsigned int *, n_face_orientations> orientation;
- orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ?
- &data.face_orientations[face_orientation[0]][0] :
- &dummy;
+ if (n_face_orientations == 1)
+ orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ?
+ &data.face_orientations[face_orientations[0]][0] :
+ &dummy;
+ else
+ {
+ for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+ {
+ // the loop breaks once an invalid_unsigned_int is hit for
+ // all cases except the exterior faces in the ECL loop (where
+ // some faces might be at the boundaries but others not)
+ if (cells[v] == numbers::invalid_unsigned_int)
+ continue;
+
+ orientation[v] =
+ (data.data.front().nodal_at_cell_boundaries == true) ?
+ &data.face_orientations[face_orientations[v]][0] :
+ &dummy;
+ }
+ }
// face_to_cell_index_hermite
std::array<const unsigned int *, n_face_orientations> index_array_hermite;
- index_array_hermite[0] =
- (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
- data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
- &data.face_to_cell_index_hermite(face_no[0], 0) :
- &dummy;
+ if (n_face_orientations == 1)
+ index_array_hermite[0] =
+ (data.data.front().nodal_at_cell_boundaries == true &&
+ fe_degree > 1 &&
+ data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
+ &data.face_to_cell_index_hermite(face_nos[0], 0) :
+ &dummy;
if (n_face_orientations > 1 &&
data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite)
{
- const unsigned int n_filled_lanes =
- dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
-
- for (unsigned int v = 1; v < n_filled_lanes; ++v)
+ for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
{
+ if (cells[v] == numbers::invalid_unsigned_int)
+ continue;
+
grad_weight[v] =
data.data.front().shape_data_on_face
- [0][fe_degree + (integrate ? (2 - (face_no[v] % 2)) :
- (1 + (face_no[v] % 2)))][v];
+ [0][fe_degree + (integrate ? (2 - (face_nos[v] % 2)) :
+ (1 + (face_nos[v] % 2)))][v];
index_array_hermite[v] =
- &data.face_to_cell_index_hermite(face_no[v], 0);
+ &data.face_to_cell_index_hermite(face_nos[v], 0);
}
}
// face_to_cell_index_nodal
std::array<const unsigned int *, n_face_orientations> index_array_nodal;
- index_array_nodal[0] =
- (data.data.front().nodal_at_cell_boundaries == true) ?
- &data.face_to_cell_index_nodal(face_no[0], 0) :
- &dummy;
+ if (n_face_orientations == 1)
+ index_array_nodal[0] =
+ (data.data.front().nodal_at_cell_boundaries == true) ?
+ &data.face_to_cell_index_nodal(face_nos[0], 0) :
+ &dummy;
if (n_face_orientations > 1 &&
(data.data.front().nodal_at_cell_boundaries == true))
{
- const unsigned int n_filled_lanes =
- dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
+ for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+ {
+ if (cells[v] == numbers::invalid_unsigned_int)
+ continue;
- for (unsigned int v = 1; v < n_filled_lanes; ++v)
- index_array_nodal[v] =
- &data.face_to_cell_index_nodal(face_no[v], 0);
+ index_array_nodal[v] =
+ &data.face_to_cell_index_nodal(face_nos[v], 0);
+ }
}
const auto reorientate = [&](const unsigned int v, const unsigned int i) {
- return (dim < 3 || face_orientation[0] == 0 ||
+ return (dim < 3 ||
+ face_orientations[n_face_orientations == 1 ? 0 : v] == 0 ||
subface_index < GeometryInfo<dim>::max_children_per_cell) ?
i :
orientation[v][i];
};
- // case 1: contiguous and interleaved indices
- if (((do_gradients == false &&
- data.data.front().nodal_at_cell_boundaries == true) ||
- (data.element_type ==
- MatrixFreeFunctions::tensor_symmetric_hermite &&
- fe_degree > 1)) &&
- dof_info.index_storage_variants[dof_access_index][cell] ==
- MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
- interleaved_contiguous)
+ if ((do_gradients == false &&
+ data.data.front().nodal_at_cell_boundaries == true) ||
+ (data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite &&
+ fe_degree > 1))
{
- AssertDimension(
- dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
- VectorizedArrayType::size());
- const unsigned int dof_index =
- dof_info
- .dof_indices_contiguous[dof_access_index]
- [cell * VectorizedArrayType::size()] +
- dof_info.component_dof_indices_offset[active_fe_index]
+ // case 1: contiguous and interleaved indices
+ if (n_face_orientations == 1 &&
+ dof_info.index_storage_variants[dof_access_index][cell] ==
+ MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous)
+ {
+ AssertDimension(n_face_orientations, 1);
+
+ AssertDimension(
+ dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+ VectorizedArrayType::size());
+ const unsigned int dof_index =
+ dof_info
+ .dof_indices_contiguous[dof_access_index]
+ [cell * VectorizedArrayType::size()] +
+ dof_info
+ .component_dof_indices_offset[active_fe_index]
[first_selected_component] *
- VectorizedArrayType::size();
+ VectorizedArrayType::size();
- if (fe_degree > 1 && do_gradients == true)
- {
- for (unsigned int i = 0; i < dofs_per_face; ++i)
+ if (fe_degree > 1 && do_gradients == true)
{
- if (n_face_orientations == 1)
- {
- const unsigned int ind1 = index_array_hermite[0][2 * i];
- const unsigned int ind2 =
- index_array_hermite[0][2 * i + 1];
- AssertIndexRange(ind1, data.dofs_per_component_on_cell);
- AssertIndexRange(ind2, data.dofs_per_component_on_cell);
- const unsigned int i_ = reorientate(0, i);
- for (unsigned int comp = 0; comp < n_components; ++comp)
- function_1a(
- temp1[i_ + 2 * comp * dofs_per_face],
- temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face],
- global_vector_ptr + dof_index +
- (ind1 + comp * static_dofs_per_component) *
- VectorizedArrayType::size(),
- global_vector_ptr + dof_index +
- (ind2 + comp * static_dofs_per_component) *
- VectorizedArrayType::size(),
- grad_weight);
- }
- else
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
{
- Assert(false, ExcNotImplemented());
-
- const unsigned int n_filled_lanes =
- dof_info
- .n_vectorization_lanes_filled[dof_access_index][cell];
-
- for (unsigned int comp = 0; comp < n_components; ++comp)
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
- function_3a(
- temp1[reorientate(v, i) + 2 * comp * dofs_per_face]
- [v],
- temp1[reorientate(v, i) + dofs_per_face +
- 2 * comp * dofs_per_face][v],
- global_vector_ptr[dof_index +
- (index_array_hermite[v][2 * i] +
- comp *
- static_dofs_per_component) *
- VectorizedArrayType::size() +
- v],
- global_vector_ptr
- [dof_index +
- (index_array_hermite[v][2 * i + 1] +
- comp * static_dofs_per_component) *
- VectorizedArrayType::size() +
- v],
- grad_weight[v]);
+ if (n_face_orientations == 1)
+ {
+ const unsigned int ind1 =
+ index_array_hermite[0][2 * i];
+ const unsigned int ind2 =
+ index_array_hermite[0][2 * i + 1];
+ AssertIndexRange(ind1,
+ data.dofs_per_component_on_cell);
+ AssertIndexRange(ind2,
+ data.dofs_per_component_on_cell);
+ const unsigned int i_ = reorientate(0, i);
+ for (unsigned int comp = 0; comp < n_components;
+ ++comp)
+ function_1a(
+ temp1[i_ + 2 * comp * dofs_per_face],
+ temp1[i_ + dofs_per_face +
+ 2 * comp * dofs_per_face],
+ global_vector_ptr + dof_index +
+ (ind1 + comp * static_dofs_per_component) *
+ VectorizedArrayType::size(),
+ global_vector_ptr + dof_index +
+ (ind2 + comp * static_dofs_per_component) *
+ VectorizedArrayType::size(),
+ grad_weight);
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
}
}
- }
- else
- {
- for (unsigned int i = 0; i < dofs_per_face; ++i)
+ else
{
- if (n_face_orientations == 1)
- {
- const unsigned int i_ = reorientate(0, i);
- const unsigned int ind = index_array_nodal[0][i];
- for (unsigned int comp = 0; comp < n_components; ++comp)
- function_1b(temp1[i_ + 2 * comp * dofs_per_face],
- global_vector_ptr + dof_index +
- (ind + comp * static_dofs_per_component) *
- VectorizedArrayType::size());
- }
- else
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
{
- Assert(false, ExcNotImplemented());
-
- const unsigned int n_filled_lanes =
- dof_info
- .n_vectorization_lanes_filled[dof_access_index][cell];
-
- for (unsigned int comp = 0; comp < n_components; ++comp)
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
- function_3b(
- temp1[reorientate(v, i) + 2 * comp * dofs_per_face]
- [v],
- global_vector_ptr[dof_index +
- (index_array_nodal[v][i] +
- comp *
- static_dofs_per_component) *
- VectorizedArrayType::size() +
- v]);
+ if (n_face_orientations == 1)
+ {
+ const unsigned int i_ = reorientate(0, i);
+ const unsigned int ind = index_array_nodal[0][i];
+ for (unsigned int comp = 0; comp < n_components;
+ ++comp)
+ function_1b(temp1[i_ + 2 * comp * dofs_per_face],
+ global_vector_ptr + dof_index +
+ (ind +
+ comp * static_dofs_per_component) *
+ VectorizedArrayType::size());
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
}
}
}
- }
- // case 2: contiguous and interleaved indices with fixed stride
- else if (((do_gradients == false &&
- data.data.front().nodal_at_cell_boundaries == true) ||
- (data.element_type ==
- MatrixFreeFunctions::tensor_symmetric_hermite &&
- fe_degree > 1)) &&
- dof_info.index_storage_variants[dof_access_index][cell] ==
- MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
- interleaved_contiguous_strided)
- {
- AssertDimension(
- dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
- VectorizedArrayType::size());
- const unsigned int *indices =
- &dof_info
- .dof_indices_contiguous[dof_access_index]
- [cell * VectorizedArrayType::size()];
- if (fe_degree > 1 && do_gradients == true)
+ // case 2: contiguous and interleaved indices with fixed stride
+ else if (n_face_orientations == 1 &&
+ dof_info.index_storage_variants[dof_access_index][cell] ==
+ MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous_strided)
{
- for (unsigned int i = 0; i < dofs_per_face; ++i)
+ AssertDimension(n_face_orientations, 1);
+
+ AssertDimension(
+ dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+ VectorizedArrayType::size());
+ const unsigned int *indices =
+ &dof_info
+ .dof_indices_contiguous[dof_access_index]
+ [cell * VectorizedArrayType::size()];
+ if (fe_degree > 1 && do_gradients == true)
{
- if (n_face_orientations == 1)
- {
- const unsigned int i_ = reorientate(0, i);
-
- const unsigned int ind1 = index_array_hermite[0][2 * i] *
- VectorizedArrayType::size();
- const unsigned int ind2 =
- index_array_hermite[0][2 * i + 1] *
- VectorizedArrayType::size();
- for (unsigned int comp = 0; comp < n_components; ++comp)
- function_2a(
- temp1[i_ + 2 * comp * dofs_per_face],
- temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face],
- global_vector_ptr + ind1 +
- comp * static_dofs_per_component *
- VectorizedArrayType::size() +
- dof_info.component_dof_indices_offset
- [active_fe_index][first_selected_component] *
- VectorizedArrayType::size(),
- global_vector_ptr + ind2 +
- comp * static_dofs_per_component *
- VectorizedArrayType::size() +
- dof_info.component_dof_indices_offset
- [active_fe_index][first_selected_component] *
- VectorizedArrayType::size(),
- grad_weight,
- indices,
- indices);
- }
- else
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
{
- Assert(false, ExcNotImplemented());
-
- const unsigned int n_filled_lanes =
- dof_info
- .n_vectorization_lanes_filled[dof_access_index][cell];
+ if (n_face_orientations == 1)
+ {
+ const unsigned int i_ = reorientate(0, i);
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
- for (unsigned int comp = 0; comp < n_components; ++comp)
- {
- const unsigned int i_ = reorientate(0, i);
- function_3a(temp1[i_ + 2 * comp * dofs_per_face][v],
+ const unsigned int ind1 =
+ index_array_hermite[0][2 * i] *
+ VectorizedArrayType::size();
+ const unsigned int ind2 =
+ index_array_hermite[0][2 * i + 1] *
+ VectorizedArrayType::size();
+ for (unsigned int comp = 0; comp < n_components;
+ ++comp)
+ function_2a(temp1[i_ + 2 * comp * dofs_per_face],
temp1[i_ + dofs_per_face +
- 2 * comp * dofs_per_face][v],
- global_vector_ptr
- [index_array_hermite[v][2 * i] *
- VectorizedArrayType::size() +
- comp * static_dofs_per_component *
- VectorizedArrayType::size() +
- dof_info.component_dof_indices_offset
- [active_fe_index]
- [first_selected_component] *
- VectorizedArrayType::size() +
- indices[v] + v],
- global_vector_ptr
- [index_array_hermite[v][2 * i + 1] *
- VectorizedArrayType::size() +
- comp * static_dofs_per_component *
- VectorizedArrayType::size() +
- dof_info.component_dof_indices_offset
- [active_fe_index]
- [first_selected_component] *
- VectorizedArrayType::size() +
- indices[v] + v],
- grad_weight[v]);
- }
+ 2 * comp * dofs_per_face],
+ global_vector_ptr + ind1 +
+ comp * static_dofs_per_component *
+ VectorizedArrayType::size() +
+ dof_info.component_dof_indices_offset
+ [active_fe_index]
+ [first_selected_component] *
+ VectorizedArrayType::size(),
+ global_vector_ptr + ind2 +
+ comp * static_dofs_per_component *
+ VectorizedArrayType::size() +
+ dof_info.component_dof_indices_offset
+ [active_fe_index]
+ [first_selected_component] *
+ VectorizedArrayType::size(),
+ grad_weight,
+ indices,
+ indices);
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
}
}
- }
- else
- {
- for (unsigned int i = 0; i < dofs_per_face; ++i)
+ else
{
- if (n_face_orientations == 1)
- {
- const unsigned int i_ = reorientate(0, i);
-
- const unsigned int ind =
- index_array_nodal[0][i] * VectorizedArrayType::size();
- for (unsigned int comp = 0; comp < n_components; ++comp)
- function_2b(
- temp1[i_ + 2 * comp * dofs_per_face],
- global_vector_ptr + ind +
- comp * static_dofs_per_component *
- VectorizedArrayType::size() +
- dof_info.component_dof_indices_offset
- [active_fe_index][first_selected_component] *
- VectorizedArrayType::size(),
- indices);
- }
- else
- {
- Assert(false, ExcNotImplemented());
-
- const unsigned int n_filled_lanes =
- dof_info
- .n_vectorization_lanes_filled[dof_access_index][cell];
-
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
- for (unsigned int comp = 0; comp < n_components; ++comp)
- function_3b(
- temp1[reorientate(v, i) + 2 * comp * dofs_per_face]
- [v],
- global_vector_ptr
- [index_array_nodal[v][i] +
- comp * static_dofs_per_component *
- VectorizedArrayType::size() +
- dof_info.component_dof_indices_offset
- [active_fe_index][first_selected_component] *
- VectorizedArrayType::size() +
- indices[v] + v]);
- }
- }
- }
- }
-
- // case 3: contiguous and interleaved indices with mixed stride
- else if (((do_gradients == false &&
- data.data.front().nodal_at_cell_boundaries == true) ||
- (data.element_type ==
- MatrixFreeFunctions::tensor_symmetric_hermite &&
- fe_degree > 1)) &&
- dof_info.index_storage_variants[dof_access_index][cell] ==
- MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
- interleaved_contiguous_mixed_strides)
- {
- const unsigned int *strides =
- &dof_info.dof_indices_interleave_strides
- [dof_access_index][cell * VectorizedArrayType::size()];
- unsigned int indices[VectorizedArrayType::size()];
- for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
- indices[v] =
- dof_info.dof_indices_contiguous
- [dof_access_index][cell * VectorizedArrayType::size() + v] +
- dof_info.component_dof_indices_offset[active_fe_index]
- [first_selected_component] *
- strides[v];
- const unsigned int n_filled_lanes =
- dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
-
- if (fe_degree > 1 && do_gradients == true)
- {
- if (n_filled_lanes == VectorizedArrayType::size())
- for (unsigned int comp = 0; comp < n_components; ++comp)
for (unsigned int i = 0; i < dofs_per_face; ++i)
{
if (n_face_orientations == 1)
{
const unsigned int i_ = reorientate(0, i);
- unsigned int ind1[VectorizedArrayType::size()];
- DEAL_II_OPENMP_SIMD_PRAGMA
- for (unsigned int v = 0;
- v < VectorizedArrayType::size();
- ++v)
- ind1[v] = indices[v] +
- (comp * static_dofs_per_component +
- index_array_hermite[0 /*TODO*/][2 * i]) *
- strides[v];
- unsigned int ind2[VectorizedArrayType::size()];
- DEAL_II_OPENMP_SIMD_PRAGMA
- for (unsigned int v = 0;
- v < VectorizedArrayType::size();
- ++v)
- ind2[v] =
- indices[v] +
- (comp * static_dofs_per_component +
- index_array_hermite[0 /*TODO*/][2 * i + 1]) *
- strides[v];
- function_2a(temp1[i_ + 2 * comp * dofs_per_face],
- temp1[i_ + dofs_per_face +
- 2 * comp * dofs_per_face],
- global_vector_ptr,
- global_vector_ptr,
- grad_weight,
- ind1,
- ind2);
+
+ const unsigned int ind = index_array_nodal[0][i] *
+ VectorizedArrayType::size();
+ for (unsigned int comp = 0; comp < n_components;
+ ++comp)
+ function_2b(temp1[i_ + 2 * comp * dofs_per_face],
+ global_vector_ptr + ind +
+ comp * static_dofs_per_component *
+ VectorizedArrayType::size() +
+ dof_info.component_dof_indices_offset
+ [active_fe_index]
+ [first_selected_component] *
+ VectorizedArrayType::size(),
+ indices);
}
else
{
Assert(false, ExcNotImplemented());
+ }
+ }
+ }
+ }
- const unsigned int n_filled_lanes =
- dof_info
- .n_vectorization_lanes_filled[dof_access_index]
- [cell];
+ // case 3: contiguous and interleaved indices with mixed stride
+ else if (n_face_orientations == 1 &&
+ dof_info.index_storage_variants[dof_access_index][cell] ==
+ MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous_mixed_strides)
+ {
+ AssertDimension(n_face_orientations, 1);
+
+ const unsigned int *strides =
+ &dof_info.dof_indices_interleave_strides
+ [dof_access_index][cell * VectorizedArrayType::size()];
+ unsigned int indices[VectorizedArrayType::size()];
+ for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+ indices[v] =
+ dof_info.dof_indices_contiguous
+ [dof_access_index][cell * VectorizedArrayType::size() + v] +
+ dof_info
+ .component_dof_indices_offset[active_fe_index]
+ [first_selected_component] *
+ strides[v];
+ const unsigned int n_filled_lanes =
+ dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ if (fe_degree > 1 && do_gradients == true)
+ {
+ if (n_filled_lanes == VectorizedArrayType::size())
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ if (n_face_orientations == 1)
+ {
+ const unsigned int i_ = reorientate(0, i);
+ unsigned int ind1[VectorizedArrayType::size()];
+ DEAL_II_OPENMP_SIMD_PRAGMA
+ for (unsigned int v = 0;
+ v < VectorizedArrayType::size();
+ ++v)
+ ind1[v] =
+ indices[v] +
+ (comp * static_dofs_per_component +
+ index_array_hermite[0 /*TODO*/][2 * i]) *
+ strides[v];
+ unsigned int ind2[VectorizedArrayType::size()];
+ DEAL_II_OPENMP_SIMD_PRAGMA
+ for (unsigned int v = 0;
+ v < VectorizedArrayType::size();
+ ++v)
+ ind2[v] =
+ indices[v] +
+ (comp * static_dofs_per_component +
+ index_array_hermite[0 /*TODO*/][2 * i + 1]) *
+ strides[v];
+ function_2a(temp1[i_ + 2 * comp * dofs_per_face],
+ temp1[i_ + dofs_per_face +
+ 2 * comp * dofs_per_face],
+ global_vector_ptr,
+ global_vector_ptr,
+ grad_weight,
+ ind1,
+ ind2);
+ }
+ else
{
- const unsigned int i_ = reorientate(v, i);
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ else
+ {
+ if (integrate == false)
+ for (unsigned int i = 0;
+ i < n_components * 2 * dofs_per_face;
+ ++i)
+ temp1[i] = VectorizedArrayType();
+
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int i_ =
+ reorientate(n_face_orientations == 1 ? 0 : v,
+ i);
function_3a(
temp1[i_ + 2 * comp * dofs_per_face][v],
temp1[i_ + dofs_per_face +
global_vector_ptr
[indices[v] +
(comp * static_dofs_per_component +
- index_array_hermite[v][2 * i]) *
+ index_array_hermite
+ [n_face_orientations == 1 ? 0 : v]
+ [2 * i]) *
strides[v]],
global_vector_ptr
[indices[v] +
(comp * static_dofs_per_component +
- index_array_hermite[v][2 * i + 1]) *
+ index_array_hermite
+ [n_face_orientations == 1 ? 0 : v]
+ [2 * i + 1]) *
strides[v]],
- grad_weight[v]);
+ grad_weight[n_face_orientations == 1 ? 0 : v]);
}
- }
}
+ }
else
{
- if (integrate == false)
- for (unsigned int i = 0;
- i < n_components * 2 * dofs_per_face;
- ++i)
- temp1[i] = VectorizedArrayType();
-
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ if (n_filled_lanes == VectorizedArrayType::size())
for (unsigned int comp = 0; comp < n_components; ++comp)
for (unsigned int i = 0; i < dofs_per_face; ++i)
{
- const unsigned int i_ =
- reorientate(n_face_orientations == 1 ? 0 : v, i);
- function_3a(
- temp1[i_ + 2 * comp * dofs_per_face][v],
- temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face]
- [v],
- global_vector_ptr
- [indices[v] +
- (comp * static_dofs_per_component +
- index_array_hermite
- [n_face_orientations == 1 ? 0 : v][2 * i]) *
- strides[v]],
- global_vector_ptr
- [indices[v] +
- (comp * static_dofs_per_component +
- index_array_hermite[n_face_orientations == 1 ?
- 0 :
- v][2 * i + 1]) *
- strides[v]],
- grad_weight[n_face_orientations == 1 ? 0 : v]);
+ if (n_face_orientations == 1)
+ {
+ unsigned int ind[VectorizedArrayType::size()];
+ DEAL_II_OPENMP_SIMD_PRAGMA
+ for (unsigned int v = 0;
+ v < VectorizedArrayType::size();
+ ++v)
+ ind[v] = indices[v] +
+ (comp * static_dofs_per_component +
+ index_array_nodal[0][i]) *
+ strides[v];
+ const unsigned int i_ = reorientate(0, i);
+ function_2b(temp1[i_ + 2 * comp * dofs_per_face],
+ global_vector_ptr,
+ ind);
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
}
+ else
+ {
+ if (integrate == false)
+ for (unsigned int i = 0;
+ i < n_components * dofs_per_face;
+ ++i)
+ temp1[i] = VectorizedArrayType();
+
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ function_3b(
+ temp1[reorientate(
+ n_face_orientations == 1 ? 0 : v, i) +
+ 2 * comp * dofs_per_face][v],
+ global_vector_ptr
+ [indices[v] +
+ (comp * static_dofs_per_component +
+ index_array_nodal
+ [n_face_orientations == 1 ? 0 : v][i]) *
+ strides[v]]);
+ }
}
}
- else
+
+ // case 4: contiguous indices without interleaving
+ else if (n_face_orientations > 1 ||
+ dof_info.index_storage_variants[dof_access_index][cell] ==
+ MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ contiguous)
{
- if (n_filled_lanes == VectorizedArrayType::size())
- for (unsigned int comp = 0; comp < n_components; ++comp)
+ const unsigned int *indices =
+ &dof_info
+ .dof_indices_contiguous[dof_access_index]
+ [cell * VectorizedArrayType::size()];
+
+ if (do_gradients == true &&
+ data.element_type ==
+ MatrixFreeFunctions::tensor_symmetric_hermite)
+ {
for (unsigned int i = 0; i < dofs_per_face; ++i)
{
- if (n_face_orientations == 1)
+ if (n_face_orientations == 1 &&
+ dof_info
+ .n_vectorization_lanes_filled[dof_access_index]
+ [cell] ==
+ VectorizedArrayType::size())
{
- unsigned int ind[VectorizedArrayType::size()];
- DEAL_II_OPENMP_SIMD_PRAGMA
- for (unsigned int v = 0;
- v < VectorizedArrayType::size();
- ++v)
- ind[v] =
- indices[v] + (comp * static_dofs_per_component +
- index_array_nodal[0 /*TODO*/][i]) *
- strides[v];
+ const unsigned int ind1 =
+ index_array_hermite[0][2 * i];
+ const unsigned int ind2 =
+ index_array_hermite[0][2 * i + 1];
const unsigned int i_ = reorientate(0, i);
- function_2b(temp1[i_ + 2 * comp * dofs_per_face],
- global_vector_ptr,
- ind);
+
+ for (unsigned int comp = 0; comp < n_components;
+ ++comp)
+ function_2a(
+ temp1[i_ + 2 * comp * dofs_per_face],
+ temp1[i_ + dofs_per_face +
+ 2 * comp * dofs_per_face],
+ global_vector_ptr +
+ comp * static_dofs_per_component + ind1 +
+ dof_info.component_dof_indices_offset
+ [active_fe_index][first_selected_component],
+ global_vector_ptr +
+ comp * static_dofs_per_component + ind2 +
+ dof_info.component_dof_indices_offset
+ [active_fe_index][first_selected_component],
+ grad_weight,
+ indices,
+ indices);
+ }
+ else if (n_face_orientations == 1)
+ {
+ const unsigned int ind1 =
+ index_array_hermite[0][2 * i];
+ const unsigned int ind2 =
+ index_array_hermite[0][2 * i + 1];
+ const unsigned int i_ = reorientate(0, i);
+
+ const unsigned int n_filled_lanes =
+ dof_info
+ .n_vectorization_lanes_filled[dof_access_index]
+ [cell];
+
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ for (unsigned int comp = 0; comp < n_components;
+ ++comp)
+ function_3a(
+ temp1[i_ + 2 * comp * dofs_per_face][v],
+ temp1[i_ + dofs_per_face +
+ 2 * comp * dofs_per_face][v],
+ global_vector_ptr
+ [comp * static_dofs_per_component + ind1 +
+ dof_info.component_dof_indices_offset
+ [active_fe_index]
+ [first_selected_component] +
+ indices[v]],
+ global_vector_ptr
+ [comp * static_dofs_per_component + ind2 +
+ dof_info.component_dof_indices_offset
+ [active_fe_index]
+ [first_selected_component] +
+ indices[v]],
+ grad_weight[v]);
+
+ if (integrate == false)
+ for (unsigned int v = n_filled_lanes;
+ v < VectorizedArrayType::size();
+ ++v)
+ for (unsigned int comp = 0; comp < n_components;
+ ++comp)
+ {
+ temp1[i_ + 2 * comp * dofs_per_face][v] = 0.0;
+ temp1[i_ + dofs_per_face +
+ 2 * comp * dofs_per_face][v] = 0.0;
+ }
}
else
{
Assert(false, ExcNotImplemented());
-
const unsigned int n_filled_lanes =
dof_info
.n_vectorization_lanes_filled[dof_access_index]
[cell];
for (unsigned int v = 0; v < n_filled_lanes; ++v)
- function_3b(
- temp1[reorientate(v, i) +
- 2 * comp * dofs_per_face][v],
- global_vector_ptr[indices[v] +
- (comp *
- static_dofs_per_component +
- index_array_nodal[v][i]) *
- strides[v]]);
+ for (unsigned int comp = 0; comp < n_components;
+ ++comp)
+ function_3a(
+ temp1[reorientate(v, i) +
+ 2 * comp * dofs_per_face][v],
+ temp1[reorientate(v, i) + dofs_per_face +
+ 2 * comp * dofs_per_face][v],
+ global_vector_ptr
+ [comp * static_dofs_per_component +
+ index_array_hermite[v][2 * i] +
+ dof_info.component_dof_indices_offset
+ [active_fe_index]
+ [first_selected_component] +
+ indices[v]],
+ global_vector_ptr
+ [comp * static_dofs_per_component +
+ index_array_hermite[v][2 * i + 1] +
+ dof_info.component_dof_indices_offset
+ [active_fe_index]
+ [first_selected_component] +
+ indices[v]],
+ grad_weight[v]);
}
}
+ }
else
{
- if (integrate == false)
- for (unsigned int i = 0; i < n_components * dofs_per_face;
- ++i)
- temp1[i] = VectorizedArrayType();
-
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
for (unsigned int comp = 0; comp < n_components; ++comp)
- for (unsigned int i = 0; i < dofs_per_face; ++i)
- function_3b(
- temp1[reorientate(n_face_orientations == 1 ? 0 : v,
- i) +
- 2 * comp * dofs_per_face][v],
- global_vector_ptr
- [indices[v] +
- (comp * static_dofs_per_component +
- index_array_nodal
- [n_face_orientations == 1 ? 0 : v][i]) *
- strides[v]]);
- }
- }
- }
-
- // case 4: contiguous indices without interleaving
- else if (((do_gradients == false &&
- data.data.front().nodal_at_cell_boundaries == true) ||
- (data.element_type ==
- MatrixFreeFunctions::tensor_symmetric_hermite &&
- fe_degree > 1)) &&
- dof_info.index_storage_variants[dof_access_index][cell] ==
- MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
- contiguous &&
- dof_info.n_vectorization_lanes_filled[dof_access_index][cell] ==
- VectorizedArrayType::size())
- {
- const unsigned int *indices =
- &dof_info
- .dof_indices_contiguous[dof_access_index]
- [cell * VectorizedArrayType::size()];
-
- if (do_gradients == true &&
- data.element_type ==
- MatrixFreeFunctions::tensor_symmetric_hermite)
- {
- for (unsigned int i = 0; i < dofs_per_face; ++i)
- {
- if (n_face_orientations == 1)
- {
- const unsigned int ind1 = index_array_hermite[0][2 * i];
- const unsigned int ind2 =
- index_array_hermite[0][2 * i + 1];
- const unsigned int i_ = reorientate(0, i);
-
- for (unsigned int comp = 0; comp < n_components; ++comp)
- function_2a(
- temp1[i_ + 2 * comp * dofs_per_face],
- temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face],
- global_vector_ptr + comp * static_dofs_per_component +
- ind1 +
- dof_info.component_dof_indices_offset
- [active_fe_index][first_selected_component],
- global_vector_ptr + comp * static_dofs_per_component +
- ind2 +
- dof_info.component_dof_indices_offset
- [active_fe_index][first_selected_component],
- grad_weight,
- indices,
- indices);
- }
- else
- {
- Assert(false, ExcNotImplemented());
+ {
+ if (n_face_orientations == 1 &&
+ dof_info
+ .n_vectorization_lanes_filled[dof_access_index]
+ [cell] ==
+ VectorizedArrayType::size())
+ {
+ const unsigned int ind = index_array_nodal[0][i];
+ const unsigned int i_ = reorientate(0, i);
+
+ function_2b(
+ temp1[i_ + 2 * comp * dofs_per_face],
+ global_vector_ptr +
+ comp * static_dofs_per_component + ind +
+ dof_info.component_dof_indices_offset
+ [active_fe_index][first_selected_component],
+ indices);
+ }
+ else if (n_face_orientations == 1)
+ {
+ const unsigned int ind = index_array_nodal[0][i];
+ const unsigned int i_ = reorientate(0, i);
- const unsigned int n_filled_lanes =
- dof_info
- .n_vectorization_lanes_filled[dof_access_index][cell];
+ const unsigned int n_filled_lanes =
+ dof_info
+ .n_vectorization_lanes_filled[dof_access_index]
+ [cell];
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
- for (unsigned int comp = 0; comp < n_components; ++comp)
- function_3a(
- temp1[reorientate(v, i) + 2 * comp * dofs_per_face]
- [v],
- temp1[reorientate(v, i) + dofs_per_face +
- 2 * comp * dofs_per_face][v],
- global_vector_ptr
- [comp * static_dofs_per_component +
- index_array_hermite[v][2 * i] +
- dof_info.component_dof_indices_offset
- [active_fe_index][first_selected_component] +
- indices[v]],
- global_vector_ptr
- [comp * static_dofs_per_component +
- index_array_hermite[v][2 * i + 1] +
- dof_info.component_dof_indices_offset
- [active_fe_index][first_selected_component] +
- indices[v]],
- grad_weight[v]);
- }
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ function_3b(
+ temp1[i_ + 2 * comp * dofs_per_face][v],
+ global_vector_ptr
+ [comp * static_dofs_per_component + ind +
+ dof_info.component_dof_indices_offset
+ [active_fe_index]
+ [first_selected_component] +
+ indices[v]]);
+
+ if (integrate == false)
+ for (unsigned int v = n_filled_lanes;
+ v < VectorizedArrayType::size();
+ ++v)
+ temp1[i_ + 2 * comp * dofs_per_face][v] = 0.0;
+ }
+ else
+ {
+ for (unsigned int v = 0;
+ v < VectorizedArrayType::size();
+ ++v)
+ if (cells[v] != numbers::invalid_unsigned_int)
+ function_3b(
+ temp1[reorientate(v, i) +
+ 2 * comp * dofs_per_face][v],
+ global_vector_ptr
+ [index_array_nodal[v][i] +
+ comp * static_dofs_per_component +
+ dof_info.component_dof_indices_offset
+ [active_fe_index]
+ [first_selected_component] +
+ dof_info.dof_indices_contiguous
+ [dof_access_index][cells[v]]]);
+ }
+ }
}
}
else
{
- for (unsigned int i = 0; i < dofs_per_face; ++i)
- for (unsigned int comp = 0; comp < n_components; ++comp)
- {
- if (n_face_orientations == 1)
- {
- const unsigned int ind = index_array_nodal[0][i];
- const unsigned int i_ = reorientate(0, i);
-
- function_2b(
- temp1[i_ + 2 * comp * dofs_per_face],
- global_vector_ptr + comp * static_dofs_per_component +
- ind +
- dof_info.component_dof_indices_offset
- [active_fe_index][first_selected_component],
- indices);
- }
- else
- {
- Assert(false, ExcNotImplemented());
-
- const unsigned int n_filled_lanes =
- dof_info
- .n_vectorization_lanes_filled[dof_access_index]
- [cell];
-
- for (unsigned int v = 0; v < n_filled_lanes; ++v)
- function_3b(
- temp1[reorientate(v, i) + 2 * comp * dofs_per_face]
- [v],
- global_vector_ptr
- [comp * static_dofs_per_component +
- index_array_nodal[v][i] +
- dof_info.component_dof_indices_offset
- [active_fe_index][first_selected_component] +
- indices[v]]);
- }
- }
+ // case 5: default vector access
+ function_5(temp1);
+ return false;
}
}
-
- // case 5: default vector access
else
{
+ // case 5: default vector access
function_5(temp1);
return false;
}
function_0(temp1, dofs_per_face);
if (!integrate &&
- (face_orientation[0] > 0 &&
+ (face_orientations[0] > 0 &&
subface_index < GeometryInfo<dim>::max_children_per_cell))
{
- AssertDimension(face_orientation.size(), 1);
- adjust_for_face_orientation(face_orientation[0],
+ AssertDimension(face_orientations.size(), 1);
+ adjust_for_face_orientation(face_orientations[0],
orientation_map,
false,
do_values,
std::array<unsigned int, VectorizedArrayType::size()>
get_cell_ids() const;
+
+ /**
+ * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is
+ * associated with.
+ */
+ std::array<unsigned int, VectorizedArrayType::size()>
+ get_cell_or_face_ids() const;
+
//@}
/**
* points is inaccurate and this value must be used instead.
*/
const unsigned int n_q_points;
+
+
+private:
+ /**
+ * Return face number of each face of the current face batch.
+ */
+ std::array<unsigned int, VectorizedArrayType::size()>
+ compute_face_no_data();
+
+ /**
+ * Determine the orientation of each face of the current face batch.
+ */
+ std::array<unsigned int, VectorizedArrayType::size()>
+ compute_face_orientations();
};
+template <int dim,
+ int n_components_,
+ typename Number,
+ bool is_face,
+ typename VectorizedArrayType>
+std::array<unsigned int, VectorizedArrayType::size()>
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+ get_cell_or_face_ids() const
+{
+ const unsigned int v_len = VectorizedArrayType::size();
+ std::array<unsigned int, VectorizedArrayType::size()> cells;
+
+ // initialize array
+ for (unsigned int i = 0; i < v_len; ++i)
+ cells[i] = numbers::invalid_unsigned_int;
+
+ if (is_face &&
+ this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ this->is_interior_face == false)
+ {
+ // cell-based face-loop: plus face
+ for (unsigned int i = 0; i < v_len; i++)
+ {
+ // compute actual (non vectorized) cell ID
+ const unsigned int cell_this = this->cell * v_len + i;
+ // compute face ID
+ unsigned int fn =
+ this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
+ this->face_no,
+ i);
+
+ if (fn == numbers::invalid_unsigned_int)
+ continue; // invalid face ID: no neighbor on boundary
+
+ // get cell ID on both sides of face
+ auto cell_m = this->matrix_info->get_face_info(fn / v_len)
+ .cells_interior[fn % v_len];
+ auto cell_p = this->matrix_info->get_face_info(fn / v_len)
+ .cells_exterior[fn % v_len];
+
+ // compare the IDs with the given cell ID
+ if (cell_m == cell_this)
+ cells[i] = cell_p; // neighbor has the other ID
+ else if (cell_p == cell_this)
+ cells[i] = cell_m;
+ }
+ }
+ else
+ {
+ for (unsigned int i = 0; i < v_len; ++i)
+ cells[i] = cell * v_len + i;
+ }
+
+ return cells;
+}
+
+
+
template <int dim,
int n_components_,
typename Number,
ExcMessage(
"Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
- if (!internal::FEFaceEvaluationSelector<dim,
- fe_degree,
- n_q_points_1d,
- n_components,
- Number,
- VectorizedArrayType>::
- gather_evaluate(input_vector.begin(),
- *this->data,
- *this->dof_info,
- this->begin_values(),
- this->begin_gradients(),
- this->scratch_data,
- evaluation_flag & EvaluationFlags::values,
- evaluation_flag & EvaluationFlags::gradients,
- this->active_fe_index,
- this->first_selected_component,
- this->cell,
- std::array<unsigned int, 1>{{this->face_no}},
- this->subface_index,
- this->dof_access_index,
- std::array<unsigned int, 1>{{this->face_orientation}},
- this->mapping_data->descriptor[this->active_fe_index]
- .face_orientations))
+ const auto fu = [&]() {
+ if (this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ this->is_interior_face == false)
+ {
+ const auto cells = this->get_cell_or_face_ids();
+ const auto face_nos = this->compute_face_no_data();
+ const auto face_orientations = this->compute_face_orientations();
+
+ return internal::FEFaceEvaluationSelector<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType>::
+ gather_evaluate(input_vector.begin(),
+ *this->data,
+ *this->dof_info,
+ this->begin_values(),
+ this->begin_gradients(),
+ this->scratch_data,
+ evaluation_flag & EvaluationFlags::values,
+ evaluation_flag & EvaluationFlags::gradients,
+ this->active_fe_index,
+ this->first_selected_component,
+ cells,
+ face_nos,
+ this->subface_index,
+ this->dof_access_index,
+ face_orientations,
+ this->mapping_data->descriptor[this->active_fe_index]
+ .face_orientations);
+ }
+ else
+ {
+ return internal::FEFaceEvaluationSelector<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType>::
+ gather_evaluate(input_vector.begin(),
+ *this->data,
+ *this->dof_info,
+ this->begin_values(),
+ this->begin_gradients(),
+ this->scratch_data,
+ evaluation_flag & EvaluationFlags::values,
+ evaluation_flag & EvaluationFlags::gradients,
+ this->active_fe_index,
+ this->first_selected_component,
+ std::array<unsigned int, 1>{{this->cell}},
+ std::array<unsigned int, 1>{{this->face_no}},
+ this->subface_index,
+ this->dof_access_index,
+ std::array<unsigned int, 1>{{this->face_orientation}},
+ this->mapping_data->descriptor[this->active_fe_index]
+ .face_orientations);
+ }
+ };
+
+ if (!fu())
{
this->read_dof_values(input_vector);
this->evaluate(evaluation_flag);
"Use integrate() followed by distribute_local_to_global() "
"instead.");
+ Assert((this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ this->is_interior_face == false) == false,
+ ExcNotImplemented());
+
if (!internal::FEFaceEvaluationSelector<dim,
fe_degree,
n_q_points_1d,
evaluation_flag & EvaluationFlags::gradients,
this->active_fe_index,
this->first_selected_component,
- this->cell,
+ std::array<unsigned int, 1>{{this->cell}},
std::array<unsigned int, 1>{{this->face_no}},
this->subface_index,
this->dof_access_index,
+template <int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components_,
+ typename Number,
+ typename VectorizedArrayType>
+std::array<unsigned int, VectorizedArrayType::size()>
+FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components_,
+ Number,
+ VectorizedArrayType>::compute_face_no_data()
+{
+ std::array<unsigned int, VectorizedArrayType::size()> face_no_data;
+
+ if (this->dof_access_index !=
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ||
+ this->is_interior_face == true)
+ {
+ std::fill(face_no_data.begin(),
+ face_no_data.begin() +
+ this->dof_info->n_vectorization_lanes_filled
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [this->cell],
+ this->face_no);
+ }
+ else
+ {
+ std::fill(face_no_data.begin(),
+ face_no_data.end(),
+ numbers::invalid_unsigned_int);
+
+ for (unsigned int i = 0;
+ i < this->dof_info->n_vectorization_lanes_filled
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [this->cell];
+ i++)
+ {
+ // compute actual (non vectorized) cell ID
+ const unsigned int cell_this =
+ this->cell * VectorizedArrayType::size() + i;
+ // compute face ID
+ const unsigned int face_index =
+ this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
+ this->face_no,
+ i);
+
+ Assert(face_index != numbers::invalid_unsigned_int,
+ ExcNotInitialized());
+
+ // get cell ID on both sides of face
+ auto cell_m =
+ this->matrix_info
+ ->get_face_info(face_index / VectorizedArrayType::size())
+ .cells_interior[face_index % VectorizedArrayType::size()];
+
+ // compare the IDs with the given cell ID
+ face_no_data[i] =
+ (cell_m == cell_this) ?
+ this->matrix_info
+ ->get_face_info(face_index / VectorizedArrayType::size())
+ .exterior_face_no :
+ this->matrix_info
+ ->get_face_info(face_index / VectorizedArrayType::size())
+ .interior_face_no;
+ }
+ }
+
+ return face_no_data;
+}
+
+
+
+template <int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components_,
+ typename Number,
+ typename VectorizedArrayType>
+std::array<unsigned int, VectorizedArrayType::size()>
+FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components_,
+ Number,
+ VectorizedArrayType>::compute_face_orientations()
+{
+ std::array<unsigned int, VectorizedArrayType::size()> face_no_data;
+
+ if (this->dof_access_index !=
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ||
+ this->is_interior_face == true)
+ {
+ Assert(false, ExcNotImplemented());
+ }
+ else
+ {
+ std::fill(face_no_data.begin(), face_no_data.end(), 0);
+
+ if (dim == 3)
+ {
+ for (unsigned int i = 0;
+ i < this->dof_info->n_vectorization_lanes_filled
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [this->cell];
+ i++)
+ {
+ // compute actual (non vectorized) cell ID
+ const unsigned int cell_this =
+ this->cell * VectorizedArrayType::size() + i;
+ // compute face ID
+ const unsigned int face_index =
+ this->matrix_info->get_cell_and_face_to_plain_faces()(
+ this->cell, this->face_no, i);
+
+ Assert(face_index != numbers::invalid_unsigned_int,
+ ExcNotInitialized());
+
+ const unsigned int macro =
+ face_index / VectorizedArrayType::size();
+ const unsigned int lane =
+ face_index % VectorizedArrayType::size();
+
+ const auto &faces = this->matrix_info->get_face_info(macro);
+
+ // get cell ID on both sides of face
+ auto cell_m = faces.cells_interior[lane];
+
+ const bool is_interior_face = cell_m != cell_this;
+ const bool fo_interior_face = faces.face_orientation >= 8;
+
+ unsigned int face_orientation = faces.face_orientation % 8;
+
+ if (is_interior_face != fo_interior_face)
+ {
+ // invert (see also:
+ // Triangulation::update_periodic_face_map())
+ static const std::array<unsigned int, 8> table{
+ {0, 1, 0, 3, 6, 5, 4, 7}};
+
+ face_orientation = table[face_orientation];
+ }
+
+ // compare the IDs with the given cell ID
+ face_no_data[i] = face_orientation;
+ }
+ }
+ }
+
+ return face_no_data;
+}
+
+
+
/*------------------------- end FEFaceEvaluation ------------------------- */
projects / dealii.git / commitdiff