contiguous
};
+ /**
+ * Enum used to distinguish the data arrays for the vectorization type
+ * in cells and faces.
+ */
+ enum DoFAccessIndex : unsigned char
+ {
+ /**
+ * The data index for the faces designated as interior
+ */
+ dof_access_face_interior = 0,
+ /**
+ * The data index for the faces designated as exterior
+ */
+ dof_access_face_exterior = 1,
+ /**
+ * The data index for the cells
+ */
+ dof_access_cell = 2
+ };
+
/**
* Stores the dimension of the underlying DoFHandler. Since the indices
* are not templated, this is the variable that makes the dimension
unsigned int vectorization_length;
/**
- * Stores the index storage variant of all cell batches.
+ * Stores the index storage variant of all cell and face batches.
*
* The three arrays given here address the types for the faces decorated
* as interior (0), the faces decorated with as exterior (1), and the
- * cells (2).
+ * cells (2) according to CellOrFaceAccess.
*/
std::vector<IndexStorageVariants> index_storage_variants[3];
*
* The three arrays given here address the types for the faces decorated
* as interior (0), the faces decorated with as exterior (1), and the
- * cells (2).
+ * cells (2) according to CellOrFaceAccess.
*/
std::vector<unsigned int> dof_indices_contiguous[3];
*
* The three arrays given here address the types for the faces decorated
* as interior (0), the faces decorated with as exterior (1), and the
- * cells (2).
+ * cells (2) according to CellOrFaceAccess.
*/
std::vector<unsigned char> n_vectorization_lanes_filled[3];
if (vectorization_length > 1)
AssertDimension(row_starts.size()/vectorization_length/n_components,
irregular_cells.size());
- index_storage_variants[2].resize(irregular_cells.size(),
- IndexStorageVariants::full);
- n_vectorization_lanes_filled[2].resize(irregular_cells.size());
+ index_storage_variants[dof_access_cell].resize(irregular_cells.size(),
+ IndexStorageVariants::full);
+ n_vectorization_lanes_filled[dof_access_cell].resize(irregular_cells.size());
for (unsigned int i=0; i<irregular_cells.size(); ++i)
if (irregular_cells[i] > 0)
- n_vectorization_lanes_filled[2][i] = irregular_cells[i];
+ n_vectorization_lanes_filled[dof_access_cell][i] = irregular_cells[i];
else
- n_vectorization_lanes_filled[2][i] = vectorization_length;
+ n_vectorization_lanes_filled[dof_access_cell][i] = vectorization_length;
- dof_indices_contiguous[2].resize(irregular_cells.size()*vectorization_length,
- numbers::invalid_unsigned_int);
+ dof_indices_contiguous[dof_access_cell].resize(irregular_cells.size()*vectorization_length,
+ numbers::invalid_unsigned_int);
dof_indices_interleaved.resize(dof_indices.size(), numbers::invalid_unsigned_int);
std::vector<unsigned int> index_kinds(static_cast<unsigned int>(IndexStorageVariants::contiguous)+1);
for (unsigned int i=0; i<irregular_cells.size(); ++i)
{
const unsigned int ndofs = dofs_per_cell[have_hp ? cell_active_fe_index[i] : 0];
- const unsigned int n_comp = n_vectorization_lanes_filled[2][i];
+ const unsigned int n_comp = n_vectorization_lanes_filled[dof_access_cell][i];
// check 1: Check if there are constraints -> no compression possible
bool has_constraints = false;
}
}
if (has_constraints)
- index_storage_variants[2][i] = IndexStorageVariants::full;
+ index_storage_variants[dof_access_cell][i] = IndexStorageVariants::full;
else
{
bool indices_are_contiguous = true;
indices_are_interleaved_and_contiguous)
{
for (unsigned int j=0; j<n_comp; ++j)
- dof_indices_contiguous[2][i*vectorization_length+j] =
+ dof_indices_contiguous[dof_access_cell][i*vectorization_length+j] =
this->dof_indices[row_starts[(i*vectorization_length+j)*n_components].first];
}
if (indices_are_contiguous)
{
- index_storage_variants[2][i] = IndexStorageVariants::contiguous;
+ index_storage_variants[dof_access_cell][i] = IndexStorageVariants::contiguous;
}
else
{
const unsigned int *dof_indices =
&this->dof_indices[row_starts[i*vectorization_length*n_components].first];
if (n_comp == vectorization_length)
- index_storage_variants[2][i] = IndexStorageVariants::interleaved;
+ index_storage_variants[dof_access_cell][i] = IndexStorageVariants::interleaved;
else
- index_storage_variants[2][i] = IndexStorageVariants::full;
+ 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)
+ for (unsigned int j=1; j<n_comp; ++j)
if (dof_indices[j*ndofs+k] == dof_indices[k])
{
- index_storage_variants[2][i] = IndexStorageVariants::full;
+ index_storage_variants[dof_access_cell][i] = IndexStorageVariants::full;
break;
}
- if (index_storage_variants[2][i] != IndexStorageVariants::full)
+ if (index_storage_variants[dof_access_cell][i] != IndexStorageVariants::full)
{
unsigned int *interleaved_dof_indices =
&this->dof_indices_interleaved[row_starts[i*vectorization_length*n_components].first];
}
}
}
- index_kinds[static_cast<unsigned int>(index_storage_variants[2][i])]++;
+ index_kinds[static_cast<unsigned int>(index_storage_variants[dof_access_cell][i])]++;
}
}
{
AssertDimension(length, vectorization_length);
- index_storage_variants[0].resize(faces.size(), IndexStorageVariants::full);
- dof_indices_contiguous[0].resize(faces.size()*length, numbers::invalid_unsigned_int);
- n_vectorization_lanes_filled[0].resize(faces.size());
+ index_storage_variants[dof_access_face_interior]
+ .resize(faces.size(), IndexStorageVariants::full);
+ dof_indices_contiguous[dof_access_face_interior]
+ .resize(faces.size()*length, numbers::invalid_unsigned_int);
+ n_vectorization_lanes_filled[dof_access_face_interior].resize(faces.size());
// all interior faces come before the boundary faces
- unsigned int n_plus_faces = 0;
- for (; n_plus_faces < faces.size(); ++n_plus_faces)
- if (faces[n_plus_faces].cells_exterior[0] == numbers::invalid_unsigned_int)
+ unsigned int n_exterior_faces = 0;
+ for (; n_exterior_faces < faces.size(); ++n_exterior_faces)
+ if (faces[n_exterior_faces].cells_exterior[0] == numbers::invalid_unsigned_int)
break;
- index_storage_variants[1].resize(n_plus_faces, IndexStorageVariants::full);
- dof_indices_contiguous[1].resize(n_plus_faces*length, numbers::invalid_unsigned_int);
- n_vectorization_lanes_filled[1].resize(n_plus_faces);
+ index_storage_variants[dof_access_face_exterior].resize(n_exterior_faces, IndexStorageVariants::full);
+ dof_indices_contiguous[dof_access_face_exterior].resize(n_exterior_faces*length, numbers::invalid_unsigned_int);
+ n_vectorization_lanes_filled[dof_access_face_exterior].resize(n_exterior_faces);
for (unsigned int face=0; face<faces.size(); ++face)
{
- auto face_computation = [&](const unsigned int face_index,
+ auto face_computation = [&](const DoFAccessIndex face_index,
const unsigned int *cell_indices_face)
{
bool is_contiguous = false;
numbers::invalid_unsigned_int; ++v)
{
n_vectorization_lanes_filled[face_index][face]++;
- if (index_storage_variants[2][cell_indices_face[v]/length] ==
+ if (index_storage_variants[dof_access_cell][cell_indices_face[v]/length] ==
IndexStorageVariants::contiguous)
is_contiguous = true;
- if (index_storage_variants[2][cell_indices_face[v]/length] <
+ if (index_storage_variants[dof_access_cell][cell_indices_face[v]/length] <
IndexStorageVariants::contiguous)
needs_full_storage = true;
}
if (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[2][cell_indices_face[v]];
+ dof_indices_contiguous[dof_access_cell][cell_indices_face[v]];
if (is_contiguous && !needs_full_storage)
index_storage_variants[face_index][face] = IndexStorageVariants::contiguous;
else
index_storage_variants[face_index][face] = IndexStorageVariants::full;
};
- face_computation(0, faces[face].cells_interior);
- if (face < n_plus_faces)
- face_computation(1, faces[face].cells_exterior);
+ face_computation(dof_access_face_interior, faces[face].cells_interior);
+ if (face < n_exterior_faces)
+ face_computation(dof_access_face_exterior, faces[face].cells_exterior);
}
}
types::global_dof_index counter = 0;
const unsigned int n_components = start_components.back();
- const unsigned int n_macro_cells = n_vectorization_lanes_filled[2].size();
+ const unsigned int n_macro_cells = n_vectorization_lanes_filled[dof_access_cell].size();
Assert(n_macro_cells <= (row_starts.size()-1)/vectorization_length/n_components,
ExcInternalError());
for (unsigned int cell_no=0; cell_no<n_macro_cells; ++cell_no)
cell_active_fe_index[cell_no] : 0]);
const unsigned int *dof_ind = &dof_indices[row_starts[cell_no*n_components*vectorization_length].first];
for (unsigned int i=0; i<ndofs; ++i)
- for (unsigned int j=0; j<n_vectorization_lanes_filled[2][cell_no]; ++j)
+ for (unsigned int j=0; j<n_vectorization_lanes_filled[dof_access_cell][cell_no]; ++j)
if (dof_ind[j*ndofs+i]<local_size)
if (renumbering[dof_ind[j*ndofs+i]] == numbers::invalid_dof_index)
renumbering[dof_ind[j*ndofs+i]] = counter++;
* Stores the index an FEFaceEvaluation object is currently pointing into
* (interior face, exterior face, data associated with cell).
*/
- unsigned int face_vector_access_index;
+ internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index;
/**
* Stores the current number of a face within the given cell in case
quadrature_weights (mapping_data->descriptor[active_quad_index].quadrature_weights.begin()),
cell (numbers::invalid_unsigned_int),
is_interior_face (is_interior_face),
- face_vector_access_index (is_face ? (is_interior_face ? 0 : 1) : 2),
+ dof_access_index (is_face ?
+ (is_interior_face ?
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior
+ :
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior)
+ :
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell),
cell_type (internal::MatrixFreeFunctions::general),
dof_values_initialized (false),
values_quad_initialized (false),
cell (0),
cell_type (internal::MatrixFreeFunctions::general),
is_interior_face (true),
- face_vector_access_index (numbers::invalid_unsigned_int),
+ dof_access_index (internal::MatrixFreeFunctions::DoFInfo::dof_access_cell),
dof_values_initialized (false),
values_quad_initialized (false),
gradients_quad_initialized(false),
cell (numbers::invalid_unsigned_int),
cell_type (internal::MatrixFreeFunctions::general),
is_interior_face (other.is_interior_face),
- face_vector_access_index (other.face_vector_access_index),
+ dof_access_index (other.dof_access_index),
dof_values_initialized (false),
values_quad_initialized (false),
gradients_quad_initialized(false),
cell = numbers::invalid_unsigned_int;
cell_type = internal::MatrixFreeFunctions::general;
is_interior_face = other.is_interior_face;
- face_vector_access_index = other.face_vector_access_index;
+ dof_access_index = other.dof_access_index;
// Create deep copy of mapped geometry for use in parallel...
if (other.mapped_geometry.get() != nullptr)
// Case 2: contiguous indices which use reduced storage of indices and can
// use vectorized load/store operations -> go to separate function
AssertIndexRange(cell,
- dof_info->index_storage_variants[face_vector_access_index].size());
- if (dof_info->index_storage_variants[is_face ? face_vector_access_index : 2][cell] >=
+ dof_info->index_storage_variants[dof_access_index].size());
+ if (dof_info->index_storage_variants[is_face ? dof_access_index :
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] >=
internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
{
read_write_operation_contiguous(operation, src);
constexpr unsigned int n_vectorization = VectorizedArray<Number>::n_array_elements;
const unsigned int dofs_per_component = this->data->dofs_per_component_on_cell;
- if (dof_info->index_storage_variants[is_face ? face_vector_access_index : 2][cell] ==
+ if (dof_info->index_storage_variants[is_face ? dof_access_index :
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] ==
internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::interleaved)
{
const unsigned int *dof_indices =
- &dof_info->dof_indices_interleaved[dof_info->row_starts[cell*n_vectorization*n_fe_components+first_selected_component].first];
+ &dof_info->dof_indices_interleaved[dof_info->row_starts[cell*n_fe_components*n_vectorization].first]+
+ dof_info->component_dof_indices_offset[active_fe_index][first_selected_component]*n_vectorization;
if (n_components == 1 || n_fe_components == 1)
for (unsigned int i=0; i<dofs_per_component; ++i, dof_indices += n_vectorization)
for (unsigned int comp=0; comp<n_components; ++comp)
unsigned int cells_copied[VectorizedArray<Number>::n_array_elements];
const unsigned int *cells;
unsigned int n_vectorization_actual =
- dof_info->n_vectorization_lanes_filled[face_vector_access_index][cell];
+ dof_info->n_vectorization_lanes_filled[dof_access_index][cell];
bool has_constraints = false;
if (is_face)
{
- if (face_vector_access_index == 2)
+ if (dof_access_index == internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
for (unsigned int v=0; v<n_vectorization_actual; ++v)
cells_copied[v] = cell*VectorizedArray<Number>::n_array_elements+v;
cells =
- face_vector_access_index == 2 ?
+ dof_access_index == internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ?
&cells_copied[0]
:
(is_interior_face ?
std::integral_constant<bool,std::is_same<typename VectorType::value_type,Number>::value>
vector_selector;
- const unsigned int ind = is_face ? face_vector_access_index : 2;
+ const internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex ind =
+ is_face ? dof_access_index : internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
const std::vector<unsigned int> &dof_indices_cont
= dof_info->dof_indices_contiguous[ind];
return;
this->cell = face_index;
- this->face_vector_access_index = this->is_interior_face ? 0 : 1;
- if (face_index >= this->matrix_info->get_task_info().refinement_edge_face_partition_data[0])
- this->face_vector_access_index = 0;
+ this->dof_access_index = this->is_interior_face ?
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior;
Assert (this->mapping_data != nullptr, ExcNotInitialized());
const unsigned int n_vectors = VectorizedArray<Number>::n_array_elements;
const internal::MatrixFreeFunctions::FaceToCellTopology<n_vectors> &faces =
this->face_orientation = 0;
this->subface_index = GeometryInfo<dim>::max_children_per_cell;
this->face_no = face_number;
- this->face_vector_access_index = 2;
+ this->dof_access_index = internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
const unsigned int offsets =
this->matrix_info->get_mapping_info().face_data_by_cells[this->quad_no].
internal::VectorReader<Number> reader;
- if (this->dof_info->index_storage_variants[this->face_vector_access_index][this->cell] ==
+ 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->face_vector_access_index][this->cell] ==
+ this->dof_info->n_vectorization_lanes_filled[this->dof_access_index][this->cell] ==
VectorizedArray<Number>::n_array_elements
&&
((evaluate_gradients == false && this->data->nodal_at_cell_boundaries == true) ||
fe_degree > 1)))
{
const unsigned int *indices = &this->dof_info->dof_indices_contiguous
- [this->face_vector_access_index][this->cell*VectorizedArray<Number>::n_array_elements];
+ [this->dof_access_index][this->cell*VectorizedArray<Number>::n_array_elements];
if (evaluate_gradients == true &&
this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite)
{
internal::VectorDistributorLocalToGlobal<Number> writer;
- if (this->dof_info->index_storage_variants[this->face_vector_access_index][this->cell] ==
+ 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->face_vector_access_index][this->cell] ==
+ this->dof_info->n_vectorization_lanes_filled[this->dof_access_index][this->cell] ==
VectorizedArray<Number>::n_array_elements
&&
((integrate_gradients == false && this->data->nodal_at_cell_boundaries == true) ||
fe_degree > 1)))
{
const unsigned int *indices = &this->dof_info->dof_indices_contiguous
- [this->face_vector_access_index][this->cell*VectorizedArray<Number>::n_array_elements];
+ [this->dof_access_index][this->cell*VectorizedArray<Number>::n_array_elements];
if (integrate_gradients == true &&
this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite)
::quadrature_point (const unsigned int q) const
{
AssertIndexRange (q, n_q_points);
- if (this->face_vector_access_index < 2)
+ if (this->dof_access_index < 2)
{
Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
ExcNotImplemented());
AssertIndexRange(face_info.faces[f].cells_interior[v],
n_macro_cells_before*
VectorizedArray<Number>::n_array_elements);
- if (dof_info[no].index_storage_variants[1][f] >=
+ if (dof_info[no].index_storage_variants[internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior][f] >=
internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous
&&
- dof_info[no].dof_indices_contiguous[2][face_info.faces[f].cells_exterior[v]] >=
+ dof_info[no].dof_indices_contiguous[internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][face_info.faces[f].cells_exterior[v]] >=
part.local_size())
{
const unsigned int p = face_info.faces[f].cells_exterior[v];
= shape_info(dof_info[no].global_base_element_offset+e,0,0,0);
for (unsigned int j=0; j<shape.dofs_per_component_on_face; ++j)
ghost_indices.push_back(part.local_to_global
- (dof_info[no].dof_indices_contiguous[2][p]+
+ (dof_info[no].dof_indices_contiguous[internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][p]+
i+
shape.face_to_cell_index_nodal
(face_info.faces[f].exterior_face_no,j)*
}
AssertDimension(i, dof_info[no].dofs_per_cell[0]*stride);
}
- else if (dof_info[no].index_storage_variants[1][f] <
+ else if (dof_info[no].index_storage_variants[internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior][f] <
internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
has_noncontiguous_cell = true;
}
{
AssertIndexRange(face_info.faces[f].cells_interior[v],
n_macro_cells_before*VectorizedArray<Number>::n_array_elements);
- if (dof_info[no].index_storage_variants[1][f] >=
+ if (dof_info[no].index_storage_variants[internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior][f] >=
internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous
&&
- dof_info[no].dof_indices_contiguous[2][face_info.faces[f].cells_exterior[v]] >=
+ dof_info[no].dof_indices_contiguous[internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][face_info.faces[f].cells_exterior[v]] >=
part.local_size())
{
const unsigned int p = face_info.faces[f].cells_exterior[v];
= shape_info(dof_info[no].global_base_element_offset+e,0,0,0);
for (unsigned int j=0; j<2*shape.dofs_per_component_on_face; ++j)
ghost_indices.
- push_back(part.local_to_global(dof_info[no].dof_indices_contiguous[2][p]
+ push_back(part.local_to_global(dof_info[no].dof_indices_contiguous[internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][p]
+i+
shape.face_to_cell_index_hermite
(face_info.faces[f].exterior_face_no,j)*stride));