start_components.clear();
row_starts_plain_indices.clear();
plain_dof_indices.clear();
+ dof_indices_interleaved.clear();
+ for (unsigned int i = 0; i < 3; ++i)
+ {
+ index_storage_variants[i].clear();
+ dof_indices_contiguous[i].clear();
+ dof_indices_interleave_strides[i].clear();
+ n_vectorization_lanes_filled[i].clear();
+ }
store_plain_indices = false;
cell_active_fe_index.clear();
max_fe_index = 0;
numbers::invalid_unsigned_int);
dof_indices_interleaved.resize(dof_indices.size(),
numbers::invalid_unsigned_int);
+ dof_indices_interleave_strides[dof_access_cell].resize(
+ irregular_cells.size() * vectorization_length,
+ numbers::invalid_unsigned_int);
std::vector<unsigned int> index_kinds(
- static_cast<unsigned int>(IndexStorageVariants::contiguous) + 1);
+ static_cast<unsigned int>(
+ IndexStorageVariants::interleaved_contiguous_mixed_strides) +
+ 1);
std::vector<unsigned int> offsets(vectorization_length);
for (unsigned int i = 0; i < irregular_cells.size(); ++i)
{
break;
}
}
+
bool indices_are_interleaved_and_contiguous =
(ndofs > 1 && n_comp == vectorization_length);
+
{
const unsigned int *dof_indices =
this->dof_indices.data() +
break;
}
}
+
if (indices_are_contiguous ||
indices_are_interleaved_and_contiguous)
{
.first];
}
- if (indices_are_contiguous)
+ if (indices_are_interleaved_and_contiguous)
+ {
+ Assert(n_comp == vectorization_length, ExcInternalError());
+ index_storage_variants[dof_access_cell][i] =
+ IndexStorageVariants::interleaved_contiguous;
+ for (unsigned int j = 0; j < n_comp; ++j)
+ dof_indices_interleave_strides[2][i * vectorization_length +
+ j] = n_comp;
+ }
+ else if (indices_are_contiguous)
{
index_storage_variants[dof_access_cell][i] =
IndexStorageVariants::contiguous;
+ for (unsigned int j = 0; j < n_comp; ++j)
+ dof_indices_interleave_strides[2][i * vectorization_length +
+ j] = 1;
}
else
{
+ int indices_are_interleaved_and_mixed = 2;
const unsigned int *dof_indices =
- this->dof_indices.data() +
- row_starts[i * vectorization_length * n_components].first;
- if (n_comp == vectorization_length)
- index_storage_variants[dof_access_cell][i] =
- IndexStorageVariants::interleaved;
- else
- index_storage_variants[dof_access_cell][i] =
- IndexStorageVariants::full;
-
- // do not use interleaved storage if two vectorized
- // components point to the same field (scatter not possible)
+ &this->dof_indices[row_starts[i * vectorization_length *
+ n_components]
+ .first];
+ for (unsigned int j = 0; j < n_comp; ++j)
+ offsets[j] =
+ dof_indices[j * ndofs + 1] - dof_indices[j * ndofs];
for (unsigned int k = 0; k < ndofs; ++k)
- for (unsigned int l = 0; l < n_comp; ++l)
- for (unsigned int j = l + 1; j < n_comp; ++j)
- if (dof_indices[j * ndofs + k] ==
- dof_indices[l * ndofs + k])
+ for (unsigned int j = 0; j < n_comp; ++j)
+ if (dof_indices[j * ndofs + k] !=
+ dof_indices[j * ndofs] + k * offsets[j])
+ {
+ indices_are_interleaved_and_mixed = 0;
+ break;
+ }
+ if (indices_are_interleaved_and_mixed == 2)
+ {
+ for (unsigned int j = 0; j < n_comp; ++j)
+ dof_indices_interleave_strides
+ [dof_access_cell][i * vectorization_length + j] =
+ offsets[j];
+ for (unsigned int j = 0; j < n_comp; ++j)
+ dof_indices_contiguous[dof_access_cell]
+ [i * vectorization_length + j] =
+ dof_indices[j * ndofs];
+ for (unsigned int j = 0; j < n_comp; ++j)
+ if (offsets[j] != vectorization_length)
{
- index_storage_variants[dof_access_cell][i] =
- IndexStorageVariants::full;
+ indices_are_interleaved_and_mixed = 1;
break;
}
- if (index_storage_variants[dof_access_cell][i] !=
- IndexStorageVariants::full)
+ if (indices_are_interleaved_and_mixed == 1 ||
+ n_comp != vectorization_length)
+ index_storage_variants[dof_access_cell][i] =
+ IndexStorageVariants::
+ interleaved_contiguous_mixed_strides;
+ else
+ index_storage_variants[dof_access_cell][i] =
+ IndexStorageVariants::interleaved_contiguous_strided;
+ }
+ else
{
- unsigned int *interleaved_dof_indices =
- this->dof_indices_interleaved.data() +
+ const unsigned int *dof_indices =
+ this->dof_indices.data() +
row_starts[i * vectorization_length * n_components]
.first;
+ if (n_comp == vectorization_length)
+ index_storage_variants[dof_access_cell][i] =
+ IndexStorageVariants::interleaved;
+ else
+ index_storage_variants[dof_access_cell][i] =
+ IndexStorageVariants::full;
+
+ // do not use interleaved storage if two vectorized
+ // components point to the same field (scatter not
+ // possible)
for (unsigned int k = 0; k < ndofs; ++k)
- for (unsigned int j = 0; j < n_comp; ++j)
- interleaved_dof_indices[k * n_comp + j] =
- dof_indices[j * ndofs + k];
+ for (unsigned int l = 0; l < n_comp; ++l)
+ for (unsigned int j = l + 1; j < n_comp; ++j)
+ if (dof_indices[j * ndofs + k] ==
+ dof_indices[l * ndofs + k])
+ {
+ index_storage_variants[dof_access_cell][i] =
+ IndexStorageVariants::full;
+ break;
+ }
}
}
}
index_kinds[static_cast<unsigned int>(
index_storage_variants[dof_access_cell][i])]++;
}
+
+ // Cleanup phase: we want to avoid single cells with different properties
+ // than the bulk of the domain in order to avoid extra checks in the face
+ // identification.
+
+ // Step 1: check whether the interleaved indices were only assigned to
+ // the single cell within a vectorized array.
+ auto fix_single_interleaved_indices =
+ [&](const IndexStorageVariants variant) {
+ if (index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::interleaved_contiguous_mixed_strides)] >
+ 0 &&
+ index_kinds[static_cast<unsigned int>(variant)] > 0)
+ for (unsigned int i = 0; i < irregular_cells.size(); ++i)
+ {
+ if (index_storage_variants[dof_access_cell][i] ==
+ IndexStorageVariants::
+ interleaved_contiguous_mixed_strides &&
+ n_vectorization_lanes_filled[dof_access_cell][i] == 1 &&
+ (variant != IndexStorageVariants::contiguous ||
+ dof_indices_interleave_strides[dof_access_cell]
+ [i * vectorization_length] ==
+ 1))
+ {
+ index_storage_variants[dof_access_cell][i] = variant;
+ index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::
+ interleaved_contiguous_mixed_strides)]--;
+ index_kinds[static_cast<unsigned int>(variant)]++;
+ }
+ }
+ };
+
+ fix_single_interleaved_indices(IndexStorageVariants::full);
+ fix_single_interleaved_indices(IndexStorageVariants::contiguous);
+ fix_single_interleaved_indices(IndexStorageVariants::interleaved);
+
+ unsigned int n_interleaved =
+ index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::interleaved_contiguous)] +
+ index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::interleaved_contiguous_strided)] +
+ index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::interleaved_contiguous_mixed_strides)];
+
+ // Step 2: fix single contiguous cell among others with interleaved
+ // storage
+ if (n_interleaved > 0 && index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::contiguous)])
+ for (unsigned int i = 0; i < irregular_cells.size(); ++i)
+ if (index_storage_variants[dof_access_cell][i] ==
+ IndexStorageVariants::contiguous)
+ {
+ index_storage_variants[dof_access_cell][i] =
+ IndexStorageVariants::interleaved_contiguous_mixed_strides;
+ index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::contiguous)]--;
+ index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::interleaved_contiguous_mixed_strides)]++;
+ }
+
+ n_interleaved =
+ index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::interleaved_contiguous)] +
+ index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::interleaved_contiguous_strided)] +
+ index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::interleaved_contiguous_mixed_strides)];
+
+ // Step 3: Interleaved cells are left but also some non-contiguous ones
+ // -> revert all to full storage
+ if (n_interleaved > 0 &&
+ index_kinds[static_cast<unsigned int>(IndexStorageVariants::full)] +
+ index_kinds[static_cast<unsigned int>(
+ IndexStorageVariants::interleaved)] >
+ 0)
+ for (unsigned int i = 0; i < irregular_cells.size(); ++i)
+ if (index_storage_variants[dof_access_cell][i] >
+ IndexStorageVariants::contiguous)
+ {
+ index_kinds[static_cast<unsigned int>(
+ index_storage_variants[2][i])]--;
+ if (n_vectorization_lanes_filled[dof_access_cell][i] ==
+ vectorization_length)
+ index_storage_variants[dof_access_cell][i] =
+ IndexStorageVariants::interleaved;
+ else
+ index_storage_variants[dof_access_cell][i] =
+ IndexStorageVariants::full;
+ index_kinds[static_cast<unsigned int>(
+ index_storage_variants[dof_access_cell][i])]++;
+ }
+
+ // Step 4: Copy the interleaved indices into their own data structure
+ for (unsigned int i = 0; i < irregular_cells.size(); ++i)
+ if (index_storage_variants[dof_access_cell][i] ==
+ IndexStorageVariants::interleaved)
+ {
+ const unsigned int ndofs =
+ dofs_per_cell[have_hp ? cell_active_fe_index[i] : 0];
+ const unsigned int *dof_indices =
+ &this->dof_indices
+ [row_starts[i * vectorization_length * n_components].first];
+ unsigned int *interleaved_dof_indices =
+ &this->dof_indices_interleaved
+ [row_starts[i * vectorization_length * n_components].first];
+ for (unsigned int k = 0; k < ndofs; ++k)
+ for (unsigned int j = 0; j < vectorization_length; ++j)
+ interleaved_dof_indices[k * vectorization_length + j] =
+ dof_indices[j * ndofs + k];
+ }
}
faces.size(), IndexStorageVariants::full);
dof_indices_contiguous[dof_access_face_interior].resize(
faces.size() * length, numbers::invalid_unsigned_int);
+ dof_indices_interleave_strides[dof_access_face_interior].resize(
+ faces.size() * length, numbers::invalid_unsigned_int);
n_vectorization_lanes_filled[dof_access_face_interior].resize(
faces.size());
n_exterior_faces, IndexStorageVariants::full);
dof_indices_contiguous[dof_access_face_exterior].resize(
n_exterior_faces * length, numbers::invalid_unsigned_int);
+ dof_indices_interleave_strides[dof_access_face_exterior].resize(
+ faces.size() * length, numbers::invalid_unsigned_int);
n_vectorization_lanes_filled[dof_access_face_exterior].resize(
n_exterior_faces);
auto face_computation = [&](const DoFAccessIndex face_index,
const unsigned int * cell_indices_face) {
bool is_contiguous = false;
+ bool is_interleaved = false;
bool needs_full_storage = false;
for (unsigned int v = 0;
v < length &&
++v)
{
n_vectorization_lanes_filled[face_index][face]++;
+ if (index_storage_variants[dof_access_cell]
+ [cell_indices_face[v] / length] >=
+ IndexStorageVariants::interleaved_contiguous)
+ is_interleaved = true;
if (index_storage_variants[dof_access_cell]
[cell_indices_face[v] / length] ==
IndexStorageVariants::contiguous)
is_contiguous = true;
+ if (index_storage_variants[dof_access_cell]
+ [cell_indices_face[v] / length] >=
+ IndexStorageVariants::contiguous)
+ dof_indices_interleave_strides[face_index][face * length +
+ v] =
+ dof_indices_interleave_strides[dof_access_cell]
+ [cell_indices_face[v]];
if (index_storage_variants[dof_access_cell]
[cell_indices_face[v] / length] <
IndexStorageVariants::contiguous)
needs_full_storage = true;
}
- if (is_contiguous)
+ Assert(!(is_interleaved && is_contiguous),
+ ExcMessage("Unsupported index compression found"));
+
+ if (is_interleaved || is_contiguous)
for (unsigned int v = 0;
v < n_vectorization_lanes_filled[face_index][face];
++v)
dof_indices_contiguous[face_index][face * length + v] =
dof_indices_contiguous[dof_access_cell][cell_indices_face[v]];
- if (is_contiguous && !needs_full_storage)
+ if (is_interleaved)
+ {
+ bool is_also_contiguous =
+ n_vectorization_lanes_filled[face_index][face] == length;
+ for (unsigned int v = 0;
+ v < n_vectorization_lanes_filled[face_index][face];
+ ++v)
+ if (dof_indices_contiguous[face_index][face * length + v] !=
+ dof_indices_contiguous[face_index][face * length] + v ||
+ dof_indices_interleave_strides[dof_access_cell]
+ [cell_indices_face[v]] !=
+ length)
+ is_also_contiguous = false;
+
+ if (is_also_contiguous)
+ {
+ index_storage_variants[face_index][face] =
+ IndexStorageVariants::interleaved_contiguous;
+ }
+ else
+ {
+ bool all_indices_same_offset =
+ n_vectorization_lanes_filled[face_index][face] == length;
+ for (unsigned int v = 0;
+ v < n_vectorization_lanes_filled[face_index][face];
+ ++v)
+ if (dof_indices_interleave_strides
+ [dof_access_cell][cell_indices_face[v]] != length)
+ all_indices_same_offset = false;
+ if (all_indices_same_offset)
+ index_storage_variants[face_index][face] =
+ IndexStorageVariants::interleaved_contiguous_strided;
+ else
+ index_storage_variants[face_index][face] =
+ IndexStorageVariants::
+ interleaved_contiguous_mixed_strides;
+ }
+ }
+ else if (is_contiguous && !needs_full_storage)
index_storage_variants[face_index][face] =
IndexStorageVariants::contiguous;
else
res = vector_access(vec, index);
}
+ template <typename VectorType>
+ void
+ process_dofs_vectorized(const unsigned int dofs_per_cell,
+ const unsigned int dof_index,
+ VectorType & vec,
+ VectorizedArray<Number> *dof_values,
+ std::integral_constant<bool, true>) const
+ {
+ const Number *vec_ptr = vec.begin() + dof_index;
+ for (unsigned int i = 0; i < dofs_per_cell;
+ ++i, vec_ptr += VectorizedArray<Number>::n_array_elements)
+ dof_values[i].load(vec_ptr);
+ }
+
+
+ template <typename VectorType>
+ void
+ process_dofs_vectorized(const unsigned int dofs_per_cell,
+ const unsigned int dof_index,
+ VectorType & vec,
+ VectorizedArray<Number> *dof_values,
+ std::integral_constant<bool, false>) const
+ {
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ dof_values[i][v] = vector_access(
+ vec, dof_index + v + i * VectorizedArray<Number>::n_array_elements);
+ }
+
template <typename VectorType>
void
process_dofs_vectorized_transpose(const unsigned int dofs_per_cell,
vector_access(vec, index) += res;
}
+ template <typename VectorType>
+ void
+ process_dofs_vectorized(const unsigned int dofs_per_cell,
+ const unsigned int dof_index,
+ VectorType & vec,
+ VectorizedArray<Number> *dof_values,
+ std::integral_constant<bool, true>) const
+ {
+ Number *vec_ptr = vec.begin() + dof_index;
+ for (unsigned int i = 0; i < dofs_per_cell;
+ ++i, vec_ptr += VectorizedArray<Number>::n_array_elements)
+ {
+ VectorizedArray<Number> tmp;
+ tmp.load(vec_ptr);
+ tmp += dof_values[i];
+ tmp.store(vec_ptr);
+ }
+ }
+
+ template <typename VectorType>
+ void
+ process_dofs_vectorized(const unsigned int dofs_per_cell,
+ const unsigned int dof_index,
+ VectorType & vec,
+ VectorizedArray<Number> *dof_values,
+ std::integral_constant<bool, false>) const
+ {
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ vector_access(vec,
+ dof_index + v +
+ i * VectorizedArray<Number>::n_array_elements) +=
+ dof_values[i][v];
+ }
+
template <typename VectorType>
void
process_dofs_vectorized_transpose(const unsigned int dofs_per_cell,
vector_access(vec, index) = res;
}
+ template <typename VectorType>
+ void
+ process_dofs_vectorized(const unsigned int dofs_per_cell,
+ const unsigned int dof_index,
+ VectorType & vec,
+ VectorizedArray<Number> *dof_values,
+ std::integral_constant<bool, true>) const
+ {
+ Number *vec_ptr = vec.begin() + dof_index;
+ for (unsigned int i = 0; i < dofs_per_cell;
+ ++i, vec_ptr += VectorizedArray<Number>::n_array_elements)
+ dof_values[i].store(vec_ptr);
+ }
+
+ template <typename VectorType>
+ void
+ process_dofs_vectorized(const unsigned int dofs_per_cell,
+ const unsigned int dof_index,
+ VectorType & vec,
+ VectorizedArray<Number> *dof_values,
+ std::integral_constant<bool, false>) const
+ {
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ vector_access(vec,
+ dof_index + v +
+ i * VectorizedArray<Number>::n_array_elements) =
+ dof_values[i][v];
+ }
+
template <typename VectorType>
void
process_dofs_vectorized_transpose(const unsigned int dofs_per_cell,
const std::vector<unsigned int> &dof_indices_cont =
dof_info->dof_indices_contiguous[ind];
+
+ // Simple case: We have contiguous storage, so we can simply copy out the
+ // data
+ if (dof_info->index_storage_variants[ind][cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous)
+ {
+ const unsigned int dof_index =
+ dof_indices_cont[cell * VectorizedArray<Number>::n_array_elements] +
+ dof_info->component_dof_indices_offset[active_fe_index]
+ [first_selected_component] *
+ VectorizedArray<Number>::n_array_elements;
+ if (n_components == 1 || n_fe_components == 1)
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ operation.process_dofs_vectorized(data->dofs_per_component_on_cell,
+ dof_index,
+ *src[comp],
+ values_dofs[comp],
+ vector_selector);
+ else
+ operation.process_dofs_vectorized(data->dofs_per_component_on_cell *
+ n_components,
+ dof_index,
+ *src[0],
+ values_dofs[0],
+ vector_selector);
+ return;
+ }
+
+ // More general case: Must go through the components one by one and apply
+ // some transformations
const unsigned int vectorization_populated =
dof_info->n_vectorization_lanes_filled[ind][this->cell];
unsigned int dof_indices[VectorizedArray<Number>::n_array_elements];
dof_indices[v] =
dof_indices_cont[cell * VectorizedArray<Number>::n_array_elements + v] +
dof_info->component_dof_indices_offset[active_fe_index]
- [first_selected_component];
+ [first_selected_component] *
+ dof_info->dof_indices_interleave_strides
+ [ind][cell * VectorizedArray<Number>::n_array_elements + v];
for (unsigned int v = vectorization_populated;
v < VectorizedArray<Number>::n_array_elements;
++v)
// constraints and that the indices within each element are contiguous
if (vectorization_populated == VectorizedArray<Number>::n_array_elements)
{
- if (n_components == 1 || n_fe_components == 1)
- for (unsigned int comp = 0; comp < n_components; ++comp)
- operation.process_dofs_vectorized_transpose(
- data->dofs_per_component_on_cell,
- dof_indices,
- *src[comp],
- values_dofs[comp],
- vector_selector);
+ if (dof_info->index_storage_variants[ind][cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ contiguous)
+ {
+ if (n_components == 1 || n_fe_components == 1)
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ operation.process_dofs_vectorized_transpose(
+ data->dofs_per_component_on_cell,
+ dof_indices,
+ *src[comp],
+ values_dofs[comp],
+ vector_selector);
+ else
+ operation.process_dofs_vectorized_transpose(
+ data->dofs_per_component_on_cell * n_components,
+ dof_indices,
+ *src[0],
+ &values_dofs[0][0],
+ vector_selector);
+ }
+ else if (dof_info->index_storage_variants[ind][cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous_strided)
+ {
+ if (n_components == 1 || n_fe_components == 1)
+ for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
+ {
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ operation.process_dof_gather(
+ dof_indices,
+ *src[comp],
+ i * VectorizedArray<Number>::n_array_elements,
+ values_dofs[comp][i],
+ vector_selector);
+ }
+ else
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+ ++i)
+ {
+ operation.process_dof_gather(
+ dof_indices,
+ *src[0],
+ (comp * data->dofs_per_component_on_cell + i) *
+ VectorizedArray<Number>::n_array_elements,
+ values_dofs[comp][i],
+ vector_selector);
+ }
+ }
else
- operation.process_dofs_vectorized_transpose(
- data->dofs_per_component_on_cell * n_components,
- dof_indices,
- *src[0],
- &values_dofs[0][0],
- vector_selector);
+ {
+ Assert(dof_info->index_storage_variants[ind][cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::
+ IndexStorageVariants::interleaved_contiguous_mixed_strides,
+ ExcNotImplemented());
+ const unsigned int *offsets =
+ &dof_info->dof_indices_interleave_strides
+ [ind][VectorizedArray<Number>::n_array_elements * cell];
+ if (n_components == 1 || n_fe_components == 1)
+ for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
+ {
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ operation.process_dof_gather(dof_indices,
+ *src[comp],
+ 0,
+ values_dofs[comp][i],
+ vector_selector);
+ DEAL_II_OPENMP_SIMD_PRAGMA
+ for (unsigned int v = 0;
+ v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ dof_indices[v] += offsets[v];
+ }
+ else
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+ ++i)
+ {
+ operation.process_dof_gather(dof_indices,
+ *src[0],
+ 0,
+ values_dofs[comp][i],
+ vector_selector);
+ DEAL_II_OPENMP_SIMD_PRAGMA
+ for (unsigned int v = 0;
+ v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ dof_indices[v] += offsets[v];
+ }
+ }
}
else
for (unsigned int comp = 0; comp < n_components; ++comp)
{
for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
operation.process_empty(values_dofs[comp][i]);
- if (n_components == 1 || n_fe_components == 1)
- for (unsigned int v = 0; v < vectorization_populated; ++v)
- for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
- operation.process_dof(dof_indices[v] + i,
- *src[comp],
- values_dofs[comp][i][v]);
+ if (dof_info->index_storage_variants[ind][cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ contiguous)
+ {
+ if (n_components == 1 || n_fe_components == 1)
+ for (unsigned int v = 0; v < vectorization_populated; ++v)
+ for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+ ++i)
+ operation.process_dof(dof_indices[v] + i,
+ *src[comp],
+ values_dofs[comp][i][v]);
+ else
+ for (unsigned int v = 0; v < vectorization_populated; ++v)
+ for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+ ++i)
+ operation.process_dof(dof_indices[v] + i +
+ comp *
+ data->dofs_per_component_on_cell,
+ *src[0],
+ values_dofs[comp][i][v]);
+ }
else
- for (unsigned int v = 0; v < vectorization_populated; ++v)
- for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
- operation.process_dof(dof_indices[v] + i +
- comp * data->dofs_per_component_on_cell,
- *src[0],
- values_dofs[comp][i][v]);
+ {
+ const unsigned int *offsets =
+ &dof_info->dof_indices_interleave_strides
+ [ind][VectorizedArray<Number>::n_array_elements * cell];
+ for (unsigned int v = 0; v < vectorization_populated; ++v)
+ AssertIndexRange(offsets[v],
+ VectorizedArray<Number>::n_array_elements + 1);
+ if (n_components == 1 || n_fe_components == 1)
+ for (unsigned int v = 0; v < vectorization_populated; ++v)
+ for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+ ++i)
+ operation.process_dof(dof_indices[v] + i * offsets[v],
+ *src[comp],
+ values_dofs[comp][i][v]);
+ else
+ for (unsigned int v = 0; v < vectorization_populated; ++v)
+ for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+ ++i)
+ operation.process_dof(
+ dof_indices[v] +
+ (i + comp * data->dofs_per_component_on_cell) *
+ offsets[v],
+ *src[0],
+ values_dofs[comp][i][v]);
+ }
}
}
const bool evaluate_gradients,
const bool evaluate_hessians)
{
- this->read_dof_values(input_vector);
- evaluate(this->begin_dof_values(),
- evaluate_values,
- evaluate_gradients,
- evaluate_hessians);
+ if (this->dof_info->index_storage_variants
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
+ IndexStorageVariants::interleaved_contiguous &&
+ reinterpret_cast<std::size_t>(
+ input_vector.begin() +
+ this->dof_info->dof_indices_contiguous
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [this->cell * VectorizedArray<Number>::n_array_elements]) %
+ sizeof(VectorizedArray<Number>) ==
+ 0)
+ {
+ const VectorizedArray<Number> *vec_values =
+ reinterpret_cast<const VectorizedArray<Number> *>(
+ input_vector.begin() +
+ this->dof_info->dof_indices_contiguous
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [this->cell * VectorizedArray<Number>::n_array_elements] +
+ this->dof_info
+ ->component_dof_indices_offset[this->active_fe_index]
+ [this->first_selected_component] *
+ VectorizedArray<Number>::n_array_elements);
+
+ evaluate(vec_values,
+ evaluate_values,
+ evaluate_gradients,
+ evaluate_hessians);
+ }
+ else
+ {
+ this->read_dof_values(input_vector);
+ evaluate(this->begin_dof_values(),
+ evaluate_values,
+ evaluate_gradients,
+ evaluate_hessians);
+ }
}
->gradients_quad[0][0],
this->scratch_data,
integrate_values,
- integrate_gradients);
+ integrate_gradients,
+ false);
# ifdef DEBUG
this->dof_values_initialized = true;
const bool integrate_gradients,
VectorType &destination)
{
- integrate(integrate_values, integrate_gradients, this->begin_dof_values());
- this->distribute_local_to_global(destination);
+ if (this->dof_info->index_storage_variants
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
+ IndexStorageVariants::interleaved_contiguous &&
+ reinterpret_cast<std::size_t>(
+ destination.begin() +
+ this->dof_info->dof_indices_contiguous
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [this->cell * VectorizedArray<Number>::n_array_elements]) %
+ sizeof(VectorizedArray<Number>) ==
+ 0)
+ {
+ VectorizedArray<Number> *vec_values =
+ reinterpret_cast<VectorizedArray<Number> *>(
+ destination.begin() +
+ this->dof_info->dof_indices_contiguous
+ [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+ [this->cell * VectorizedArray<Number>::n_array_elements] +
+ this->dof_info
+ ->component_dof_indices_offset[this->active_fe_index]
+ [this->first_selected_component] *
+ VectorizedArray<Number>::n_array_elements);
+ SelectEvaluator<
+ dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ VectorizedArray<Number>>::integrate(*this->data,
+ vec_values,
+ this->values_quad[0],
+ this->gradients_quad[0][0],
+ this->scratch_data,
+ integrate_values,
+ integrate_gradients,
+ true);
+ }
+ else
+ {
+ integrate(integrate_values,
+ integrate_gradients,
+ this->begin_dof_values());
+ this->distribute_local_to_global(destination);
+ }
}
temp1 = this->scratch_data;
internal::VectorReader<Number> reader;
+ std::integral_constant<
+ bool,
+ std::is_same<typename VectorType::value_type, Number>::value>
+ vector_selector;
- if (this->dof_info
- ->index_storage_variants[this->dof_access_index][this->cell] ==
- internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
- contiguous &&
- this->dof_info
- ->n_vectorization_lanes_filled[this->dof_access_index][this->cell] ==
- VectorizedArray<Number>::n_array_elements &&
- ((evaluate_gradients == false &&
+ // case 1: contiguous and interleaved indices
+ if (((evaluate_gradients == false &&
this->data->nodal_at_cell_boundaries == true) ||
(this->data->element_type ==
internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
- fe_degree > 1)))
+ fe_degree > 1)) &&
+ this->dof_info
+ ->index_storage_variants[this->dof_access_index][this->cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous)
+ {
+ AssertDimension(
+ this->dof_info
+ ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
+ VectorizedArray<Number>::n_array_elements);
+ const unsigned int dof_index =
+ this->dof_info
+ ->dof_indices_contiguous[this->dof_access_index]
+ [this->cell *
+ VectorizedArray<Number>::n_array_elements] +
+ this->dof_info
+ ->component_dof_indices_offset[this->active_fe_index]
+ [this->first_selected_component] *
+ VectorizedArray<Number>::n_array_elements;
+
+ if (fe_degree > 1 && evaluate_gradients == true)
+ {
+ // we know that the gradient weights for the Hermite case on the
+ // right (side==1) are the negative from the value at the left
+ // (side==0), so we only read out one of them.
+ const VectorizedArray<Number> grad_weight =
+ this->data->shape_data_on_face[0][fe_degree + 1 + side];
+ AssertDimension(this->data->face_to_cell_index_hermite.size(1),
+ 2 * dofs_per_face);
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_hermite(this->face_no, 0);
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind1 = index_array[2 * i];
+ const unsigned int ind2 = index_array[2 * i + 1];
+ AssertIndexRange(ind1, dofs_per_cell);
+ AssertIndexRange(ind2, dofs_per_cell);
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ {
+ reader.process_dofs_vectorized(
+ 1,
+ dof_index + (ind1 + comp * static_dofs_per_component) *
+ VectorizedArray<Number>::n_array_elements,
+ input_vector,
+ temp1 + i + 2 * comp * dofs_per_face,
+ vector_selector);
+ reader.process_dofs_vectorized(
+ 1,
+ dof_index + (ind2 + comp * static_dofs_per_component) *
+ VectorizedArray<Number>::n_array_elements,
+ input_vector,
+ temp1 + dofs_per_face + i + 2 * comp * dofs_per_face,
+ vector_selector);
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
+ grad_weight *
+ (temp1[i + 2 * comp * dofs_per_face] -
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face]);
+ }
+ }
+ }
+ else
+ {
+ AssertDimension(this->data->face_to_cell_index_nodal.size(1),
+ dofs_per_face);
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_nodal(this->face_no, 0);
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind = index_array[i];
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ reader.process_dofs_vectorized(
+ 1,
+ dof_index + (ind + comp * static_dofs_per_component) *
+ VectorizedArray<Number>::n_array_elements,
+ input_vector,
+ temp1 + i + 2 * comp * dofs_per_face,
+ vector_selector);
+ }
+ }
+ }
+
+ // case 2: contiguous and interleaved indices with fixed stride
+ else if (((evaluate_gradients == false &&
+ this->data->nodal_at_cell_boundaries == true) ||
+ (this->data->element_type ==
+ internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+ fe_degree > 1)) &&
+ this->dof_info
+ ->index_storage_variants[this->dof_access_index][this->cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous_strided)
+ {
+ AssertDimension(
+ this->dof_info
+ ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
+ VectorizedArray<Number>::n_array_elements);
+ const unsigned int *indices =
+ &this->dof_info
+ ->dof_indices_contiguous[this->dof_access_index]
+ [this->cell *
+ VectorizedArray<Number>::n_array_elements];
+ if (fe_degree > 1 && evaluate_gradients == true)
+ {
+ // we know that the gradient weights for the Hermite case on the
+ // right (side==1) are the negative from the value at the left
+ // (side==0), so we only read out one of them.
+ const VectorizedArray<Number> grad_weight =
+ this->data->shape_data_on_face[0][fe_degree + 1 + side];
+ AssertDimension(this->data->face_to_cell_index_hermite.size(1),
+ 2 * dofs_per_face);
+
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_hermite(this->face_no, 0);
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind1 =
+ index_array[2 * i] * VectorizedArray<Number>::n_array_elements;
+ const unsigned int ind2 =
+ index_array[2 * i + 1] *
+ VectorizedArray<Number>::n_array_elements;
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ {
+ reader.process_dof_gather(
+ indices,
+ input_vector,
+ ind1 +
+ comp * static_dofs_per_component *
+ VectorizedArray<Number>::n_array_elements +
+ this->dof_info->component_dof_indices_offset
+ [this->active_fe_index]
+ [this->first_selected_component] *
+ VectorizedArray<Number>::n_array_elements,
+ temp1[i + 2 * comp * dofs_per_face],
+ vector_selector);
+ VectorizedArray<Number> grad;
+ reader.process_dof_gather(
+ indices,
+ input_vector,
+ ind2 +
+ comp * static_dofs_per_component *
+ VectorizedArray<Number>::n_array_elements +
+ this->dof_info->component_dof_indices_offset
+ [this->active_fe_index]
+ [this->first_selected_component] *
+ VectorizedArray<Number>::n_array_elements,
+ grad,
+ vector_selector);
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
+ grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad);
+ }
+ }
+ }
+ else
+ {
+ AssertDimension(this->data->face_to_cell_index_nodal.size(1),
+ dofs_per_face);
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_nodal(this->face_no, 0);
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind =
+ index_array[i] * VectorizedArray<Number>::n_array_elements;
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ reader.process_dof_gather(
+ indices,
+ input_vector,
+ ind +
+ comp * static_dofs_per_component *
+ VectorizedArray<Number>::n_array_elements +
+ this->dof_info->component_dof_indices_offset
+ [this->active_fe_index]
+ [this->first_selected_component] *
+ VectorizedArray<Number>::n_array_elements,
+ temp1[i + 2 * comp * dofs_per_face],
+ vector_selector);
+ }
+ }
+ }
+
+ // case 3: contiguous and interleaved indices with mixed stride
+ else if (((evaluate_gradients == false &&
+ this->data->nodal_at_cell_boundaries == true) ||
+ (this->data->element_type ==
+ internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+ fe_degree > 1)) &&
+ this->dof_info
+ ->index_storage_variants[this->dof_access_index][this->cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous_mixed_strides)
+ {
+ const unsigned int *strides =
+ &this->dof_info->dof_indices_interleave_strides
+ [this->dof_access_index]
+ [this->cell * VectorizedArray<Number>::n_array_elements];
+ unsigned int indices[VectorizedArray<Number>::n_array_elements];
+ for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ indices[v] =
+ this->dof_info->dof_indices_contiguous
+ [this->dof_access_index]
+ [this->cell * VectorizedArray<Number>::n_array_elements + v] +
+ this->dof_info
+ ->component_dof_indices_offset[this->active_fe_index]
+ [this->first_selected_component] *
+ strides[v];
+ const unsigned int nvec =
+ this->dof_info
+ ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
+
+ if (fe_degree > 1 && evaluate_gradients == true)
+ {
+ // we know that the gradient weights for the Hermite case on the
+ // right (side==1) are the negative from the value at the left
+ // (side==0), so we only read out one of them.
+ const VectorizedArray<Number> grad_weight =
+ this->data->shape_data_on_face[0][fe_degree + 1 + side];
+ AssertDimension(this->data->face_to_cell_index_hermite.size(1),
+ 2 * dofs_per_face);
+
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_hermite(this->face_no, 0);
+ if (nvec == VectorizedArray<Number>::n_array_elements)
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ unsigned int ind1[VectorizedArray<Number>::n_array_elements];
+ DEAL_II_OPENMP_SIMD_PRAGMA
+ for (unsigned int v = 0;
+ v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ ind1[v] = indices[v] + (comp * static_dofs_per_component +
+ index_array[2 * i]) *
+ strides[v];
+ unsigned int ind2[VectorizedArray<Number>::n_array_elements];
+ DEAL_II_OPENMP_SIMD_PRAGMA
+ for (unsigned int v = 0;
+ v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ ind2[v] = indices[v] + (comp * static_dofs_per_component +
+ index_array[2 * i + 1]) *
+ strides[v];
+ reader.process_dof_gather(ind1,
+ input_vector,
+ 0,
+ temp1[i + 2 * comp * dofs_per_face],
+ vector_selector);
+ VectorizedArray<Number> grad;
+ reader.process_dof_gather(
+ ind2, input_vector, 0, grad, vector_selector);
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
+ grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad);
+ }
+ else
+ {
+ for (unsigned int i = 0; i < n_components_ * 2 * dofs_per_face;
+ ++i)
+ temp1[i] = VectorizedArray<Number>();
+ for (unsigned int v = 0; v < nvec; ++v)
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind1 =
+ indices[v] + (comp * static_dofs_per_component +
+ index_array[2 * i]) *
+ strides[v];
+ const unsigned int ind2 =
+ indices[v] + (comp * static_dofs_per_component +
+ index_array[2 * i + 1]) *
+ strides[v];
+ reader.process_dof(
+ ind1,
+ const_cast<VectorType &>(input_vector),
+ temp1[i + 2 * comp * dofs_per_face][v]);
+ Number grad;
+ reader.process_dof(ind2,
+ const_cast<VectorType &>(input_vector),
+ grad);
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face][v] =
+ grad_weight[0] *
+ (temp1[i + 2 * comp * dofs_per_face][v] - grad);
+ }
+ }
+ }
+ else
+ {
+ AssertDimension(this->data->face_to_cell_index_nodal.size(1),
+ dofs_per_face);
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_nodal(this->face_no, 0);
+ if (nvec == VectorizedArray<Number>::n_array_elements)
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ unsigned int ind[VectorizedArray<Number>::n_array_elements];
+ DEAL_II_OPENMP_SIMD_PRAGMA
+ for (unsigned int v = 0;
+ v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ ind[v] = indices[v] + (comp * static_dofs_per_component +
+ index_array[i]) *
+ strides[v];
+ reader.process_dof_gather(ind,
+ input_vector,
+ 0,
+ temp1[i + 2 * comp * dofs_per_face],
+ vector_selector);
+ }
+ else
+ {
+ for (unsigned int i = 0; i < n_components_ * dofs_per_face; ++i)
+ temp1[i] = VectorizedArray<Number>();
+ for (unsigned int v = 0; v < nvec; ++v)
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind1 =
+ indices[v] +
+ (comp * static_dofs_per_component + index_array[i]) *
+ strides[v];
+ reader.process_dof(
+ ind1,
+ const_cast<VectorType &>(input_vector),
+ temp1[i + 2 * comp * dofs_per_face][v]);
+ }
+ }
+ }
+ }
+
+ // case 4: contiguous indices without interleaving
+ else if (((evaluate_gradients == false &&
+ this->data->nodal_at_cell_boundaries == true) ||
+ (this->data->element_type ==
+ internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+ fe_degree > 1)) &&
+ this->dof_info
+ ->index_storage_variants[this->dof_access_index][this->cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ contiguous &&
+ this->dof_info->n_vectorization_lanes_filled[this->dof_access_index]
+ [this->cell] ==
+ VectorizedArray<Number>::n_array_elements)
{
const unsigned int *indices =
&this->dof_info
// we know that the gradient weights for the Hermite case on the
// right (side==1) are the negative from the value at the left
// (side==0), so we only read out one of them.
- const VectorizedArray<Number> grad_weight0 =
- (side ? -1. : 1.) *
- this->data->shape_data_on_face[0][fe_degree + 1];
- const VectorizedArray<Number> grad_weight1 =
- (side ? -1. : 1.) *
- this->data->shape_data_on_face[0][fe_degree + 2];
+ const VectorizedArray<Number> grad_weight =
+ this->data->shape_data_on_face[0][fe_degree + 1 + side];
AssertDimension(this->data->face_to_cell_index_hermite.size(1),
2 * dofs_per_face);
this->dof_info->component_dof_indices_offset
[this->active_fe_index][this->first_selected_component],
temp1[i + 2 * comp * dofs_per_face],
- std::integral_constant<
- bool,
- std::is_same<typename VectorType::value_type,
- Number>::value>());
+ vector_selector);
VectorizedArray<Number> grad;
reader.process_dof_gather(
indices,
this->dof_info->component_dof_indices_offset
[this->active_fe_index][this->first_selected_component],
grad,
- std::integral_constant<
- bool,
- std::is_same<typename VectorType::value_type,
- Number>::value>());
+ vector_selector);
temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
- grad_weight0 * temp1[i + 2 * comp * dofs_per_face] +
- grad_weight1 * grad;
+ grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad);
}
}
}
this->dof_info->component_dof_indices_offset
[this->active_fe_index][this->first_selected_component],
temp1[i + comp * 2 * dofs_per_face],
- std::integral_constant<
- bool,
- std::is_same<typename VectorType::value_type,
- Number>::value>());
+ vector_selector);
}
}
}
+
+ // case 5: default vector access
else
{
this->read_dof_values(input_vector);
# endif
internal::VectorDistributorLocalToGlobal<Number> writer;
+ std::integral_constant<
+ bool,
+ std::is_same<typename VectorType::value_type, Number>::value>
+ vector_selector;
- if (this->dof_info
- ->index_storage_variants[this->dof_access_index][this->cell] ==
- internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
- contiguous &&
- this->dof_info
- ->n_vectorization_lanes_filled[this->dof_access_index][this->cell] ==
- VectorizedArray<Number>::n_array_elements &&
- ((integrate_gradients == false &&
+ // case 1: contiguous and interleaved indices
+ if (((integrate_gradients == false &&
this->data->nodal_at_cell_boundaries == true) ||
(this->data->element_type ==
internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
- fe_degree > 1)))
+ fe_degree > 1)) &&
+ this->dof_info
+ ->index_storage_variants[this->dof_access_index][this->cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous)
+ {
+ AssertDimension(
+ this->dof_info
+ ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
+ VectorizedArray<Number>::n_array_elements);
+ const unsigned int dof_index =
+ this->dof_info
+ ->dof_indices_contiguous[this->dof_access_index]
+ [this->cell *
+ VectorizedArray<Number>::n_array_elements] +
+ this->dof_info
+ ->component_dof_indices_offset[this->active_fe_index]
+ [this->first_selected_component] *
+ VectorizedArray<Number>::n_array_elements;
+
+ if (fe_degree > 1 && integrate_gradients == true)
+ {
+ // we know that the gradient weights for the Hermite case on the
+ // right (side==1) are the negative from the value at the left
+ // (side==0), so we only read out one of them.
+ const VectorizedArray<Number> grad_weight =
+ this->data->shape_data_on_face[0][fe_degree + 2 - side];
+ AssertDimension(this->data->face_to_cell_index_hermite.size(1),
+ 2 * dofs_per_face);
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_hermite(this->face_no, 0);
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind1 = index_array[2 * i];
+ const unsigned int ind2 = index_array[2 * i + 1];
+ AssertIndexRange(ind1, dofs_per_cell);
+ AssertIndexRange(ind2, dofs_per_cell);
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ {
+ VectorizedArray<Number> val =
+ temp1[i + 2 * comp * dofs_per_face] -
+ grad_weight *
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+ VectorizedArray<Number> grad =
+ grad_weight *
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+ writer.process_dofs_vectorized(
+ 1,
+ dof_index + (ind1 + comp * static_dofs_per_component) *
+ VectorizedArray<Number>::n_array_elements,
+ destination,
+ &val,
+ vector_selector);
+ writer.process_dofs_vectorized(
+ 1,
+ dof_index + (ind2 + comp * static_dofs_per_component) *
+ VectorizedArray<Number>::n_array_elements,
+ destination,
+ &grad,
+ vector_selector);
+ }
+ }
+ }
+ else
+ {
+ AssertDimension(this->data->face_to_cell_index_nodal.size(1),
+ dofs_per_face);
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_nodal(this->face_no, 0);
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind = index_array[i];
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ writer.process_dofs_vectorized(
+ 1,
+ dof_index + (ind + comp * static_dofs_per_component) *
+ VectorizedArray<Number>::n_array_elements,
+ destination,
+ temp1 + i + 2 * comp * dofs_per_face,
+ vector_selector);
+ }
+ }
+ }
+
+ // case 2: contiguous and interleaved indices with fixed stride
+ else if (((integrate_gradients == false &&
+ this->data->nodal_at_cell_boundaries == true) ||
+ (this->data->element_type ==
+ internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+ fe_degree > 1)) &&
+ this->dof_info
+ ->index_storage_variants[this->dof_access_index][this->cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous_strided)
+ {
+ AssertDimension(
+ this->dof_info
+ ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
+ VectorizedArray<Number>::n_array_elements);
+ const unsigned int *indices =
+ &this->dof_info
+ ->dof_indices_contiguous[this->dof_access_index]
+ [this->cell *
+ VectorizedArray<Number>::n_array_elements];
+ if (fe_degree > 1 && integrate_gradients == true)
+ {
+ // we know that the gradient weights for the Hermite case on the
+ // right (side==1) are the negative from the value at the left
+ // (side==0), so we only read out one of them.
+ const VectorizedArray<Number> grad_weight =
+ this->data->shape_data_on_face[0][fe_degree + 2 - side];
+ AssertDimension(this->data->face_to_cell_index_hermite.size(1),
+ 2 * dofs_per_face);
+
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_hermite(this->face_no, 0);
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind1 =
+ index_array[2 * i] * VectorizedArray<Number>::n_array_elements;
+ const unsigned int ind2 =
+ index_array[2 * i + 1] *
+ VectorizedArray<Number>::n_array_elements;
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ {
+ VectorizedArray<Number> val =
+ temp1[i + 2 * comp * dofs_per_face] -
+ grad_weight *
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+ VectorizedArray<Number> grad =
+ grad_weight *
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+ writer.process_dof_gather(
+ indices,
+ destination,
+ ind1 +
+ comp * static_dofs_per_component *
+ VectorizedArray<Number>::n_array_elements +
+ this->dof_info->component_dof_indices_offset
+ [this->active_fe_index]
+ [this->first_selected_component] *
+ VectorizedArray<Number>::n_array_elements,
+ val,
+ vector_selector);
+ writer.process_dof_gather(
+ indices,
+ destination,
+ ind2 +
+ comp * static_dofs_per_component *
+ VectorizedArray<Number>::n_array_elements +
+ this->dof_info->component_dof_indices_offset
+ [this->active_fe_index]
+ [this->first_selected_component] *
+ VectorizedArray<Number>::n_array_elements,
+ grad,
+ vector_selector);
+ }
+ }
+ }
+ else
+ {
+ AssertDimension(this->data->face_to_cell_index_nodal.size(1),
+ dofs_per_face);
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_nodal(this->face_no, 0);
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind =
+ index_array[i] * VectorizedArray<Number>::n_array_elements;
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ writer.process_dof_gather(
+ indices,
+ destination,
+ ind +
+ comp * static_dofs_per_component *
+ VectorizedArray<Number>::n_array_elements +
+ this->dof_info->component_dof_indices_offset
+ [this->active_fe_index]
+ [this->first_selected_component] *
+ VectorizedArray<Number>::n_array_elements,
+ temp1[i + 2 * comp * dofs_per_face],
+ vector_selector);
+ }
+ }
+ }
+
+ // case 3: contiguous and interleaved indices with mixed stride
+ else if (((integrate_gradients == false &&
+ this->data->nodal_at_cell_boundaries == true) ||
+ (this->data->element_type ==
+ internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+ fe_degree > 1)) &&
+ this->dof_info
+ ->index_storage_variants[this->dof_access_index][this->cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ interleaved_contiguous_mixed_strides)
+ {
+ const unsigned int *strides =
+ &this->dof_info->dof_indices_interleave_strides
+ [this->dof_access_index]
+ [this->cell * VectorizedArray<Number>::n_array_elements];
+ unsigned int indices[VectorizedArray<Number>::n_array_elements];
+ for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ indices[v] =
+ this->dof_info->dof_indices_contiguous
+ [this->dof_access_index]
+ [this->cell * VectorizedArray<Number>::n_array_elements + v] +
+ this->dof_info
+ ->component_dof_indices_offset[this->active_fe_index]
+ [this->first_selected_component] *
+ strides[v];
+ const unsigned int nvec =
+ this->dof_info
+ ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
+
+ if (fe_degree > 1 && integrate_gradients == true)
+ {
+ // we know that the gradient weights for the Hermite case on the
+ // right (side==1) are the negative from the value at the left
+ // (side==0), so we only read out one of them.
+ const VectorizedArray<Number> grad_weight =
+ this->data->shape_data_on_face[0][fe_degree + 2 - side];
+ AssertDimension(this->data->face_to_cell_index_hermite.size(1),
+ 2 * dofs_per_face);
+
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_hermite(this->face_no, 0);
+ if (nvec == VectorizedArray<Number>::n_array_elements)
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ unsigned int ind1[VectorizedArray<Number>::n_array_elements];
+ DEAL_II_OPENMP_SIMD_PRAGMA
+ for (unsigned int v = 0;
+ v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ ind1[v] = indices[v] + (comp * static_dofs_per_component +
+ index_array[2 * i]) *
+ strides[v];
+ unsigned int ind2[VectorizedArray<Number>::n_array_elements];
+ DEAL_II_OPENMP_SIMD_PRAGMA
+ for (unsigned int v = 0;
+ v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ ind2[v] = indices[v] + (comp * static_dofs_per_component +
+ index_array[2 * i + 1]) *
+ strides[v];
+ VectorizedArray<Number> val =
+ temp1[i + 2 * comp * dofs_per_face] -
+ grad_weight *
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+ VectorizedArray<Number> grad =
+ grad_weight *
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+ writer.process_dof_gather(
+ ind1, destination, 0, val, vector_selector);
+ writer.process_dof_gather(
+ ind2, destination, 0, grad, vector_selector);
+ }
+ else
+ {
+ for (unsigned int v = 0; v < nvec; ++v)
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind1 =
+ indices[v] + (comp * static_dofs_per_component +
+ index_array[2 * i]) *
+ strides[v];
+ const unsigned int ind2 =
+ indices[v] + (comp * static_dofs_per_component +
+ index_array[2 * i + 1]) *
+ strides[v];
+ Number val =
+ temp1[i + 2 * comp * dofs_per_face][v] -
+ grad_weight[0] * temp1[i + dofs_per_face +
+ 2 * comp * dofs_per_face][v];
+ Number grad =
+ grad_weight[0] *
+ temp1[i + dofs_per_face + 2 * comp * dofs_per_face][v];
+ writer.process_dof(ind1, destination, val);
+ writer.process_dof(ind2, destination, grad);
+ }
+ }
+ }
+ else
+ {
+ AssertDimension(this->data->face_to_cell_index_nodal.size(1),
+ dofs_per_face);
+ const unsigned int *index_array =
+ &this->data->face_to_cell_index_nodal(this->face_no, 0);
+ if (nvec == VectorizedArray<Number>::n_array_elements)
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ unsigned int ind[VectorizedArray<Number>::n_array_elements];
+ DEAL_II_OPENMP_SIMD_PRAGMA
+ for (unsigned int v = 0;
+ v < VectorizedArray<Number>::n_array_elements;
+ ++v)
+ ind[v] = indices[v] + (comp * static_dofs_per_component +
+ index_array[i]) *
+ strides[v];
+ writer.process_dof_gather(ind,
+ destination,
+ 0,
+ temp1[i + 2 * comp * dofs_per_face],
+ vector_selector);
+ }
+ else
+ {
+ for (unsigned int v = 0; v < nvec; ++v)
+ for (unsigned int comp = 0; comp < n_components_; ++comp)
+ for (unsigned int i = 0; i < dofs_per_face; ++i)
+ {
+ const unsigned int ind1 =
+ indices[v] +
+ (comp * static_dofs_per_component + index_array[i]) *
+ strides[v];
+ writer.process_dof(
+ ind1,
+ destination,
+ temp1[i + 2 * comp * dofs_per_face][v]);
+ }
+ }
+ }
+ }
+
+ // case 4: contiguous indices without interleaving
+ else if (((integrate_gradients == false &&
+ this->data->nodal_at_cell_boundaries == true) ||
+ (this->data->element_type ==
+ internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+ fe_degree > 1)) &&
+ this->dof_info
+ ->index_storage_variants[this->dof_access_index][this->cell] ==
+ internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+ contiguous &&
+ this->dof_info->n_vectorization_lanes_filled[this->dof_access_index]
+ [this->cell] ==
+ VectorizedArray<Number>::n_array_elements)
{
const unsigned int *indices =
&this->dof_info
// we know that the gradient weights for the Hermite case on the
// right (side==1) are the negative from the value at the left
// (side==0), so we only read out one of them.
- const VectorizedArray<Number> grad_weight0 =
- (side ? -1. : 1.) *
- this->data->shape_data_on_face[0][fe_degree + 1];
- const VectorizedArray<Number> grad_weight1 =
- (side ? -1. : 1.) *
- this->data->shape_data_on_face[0][fe_degree + 2];
+ const VectorizedArray<Number> grad_weight =
+ this->data->shape_data_on_face[0][fe_degree + 2 - side];
AssertDimension(this->data->face_to_cell_index_hermite.size(1),
2 * dofs_per_face);
const unsigned int *index_array =
for (unsigned int comp = 0; comp < n_components_; ++comp)
{
VectorizedArray<Number> val =
- temp1[i + 2 * comp * dofs_per_face] +
- grad_weight0 *
+ temp1[i + 2 * comp * dofs_per_face] -
+ grad_weight *
temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
VectorizedArray<Number> grad =
- grad_weight1 *
+ grad_weight *
temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
writer.process_dof_gather(
indices,
this->dof_info->component_dof_indices_offset
[this->active_fe_index][this->first_selected_component],
val,
- std::integral_constant<
- bool,
- std::is_same<typename VectorType::value_type,
- Number>::value>());
+ vector_selector);
writer.process_dof_gather(
indices,
destination,
this->dof_info->component_dof_indices_offset
[this->active_fe_index][this->first_selected_component],
grad,
- std::integral_constant<
- bool,
- std::is_same<typename VectorType::value_type,
- Number>::value>());
+ vector_selector);
}
}
}
this->dof_info->component_dof_indices_offset
[this->active_fe_index][this->first_selected_component],
temp1[i + 2 * comp * dofs_per_face],
- std::integral_constant<
- bool,
- std::is_same<typename VectorType::value_type,
- Number>::value>());
+ vector_selector);
}
}
}
+
+ // case 5: default vector access
else
{
internal::FEFaceNormalEvaluationImpl<dim,
projects / dealii.git / commitdiff