From: Martin Kronbichler Date: Mon, 9 Sep 2019 08:35:55 +0000 (+0200) Subject: Put implementation details of FEFaceEvaluation into separate file X-Git-Tag: v9.2.0-rc1~1045^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F8713%2Fhead;p=dealii.git Put implementation details of FEFaceEvaluation into separate file --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index de9e359c82..e77ae68e11 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -22,6 +22,7 @@ #include #include +#include #include #include @@ -1730,6 +1731,1420 @@ namespace internal } } }; + + + + // internal helper function for reading data; base version of different types + template + void + do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst) + { + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + dst[v] = src_ptr[v]; + } + + + + // internal helper function for reading data; specialized version where we + // can use a dedicated load function + template + void + do_vectorized_read(const Number *src_ptr, VectorizedArray &dst) + { + dst.load(src_ptr); + } + + + + // internal helper function for reading data; base version of different types + template + void + do_vectorized_gather(const Number2 * src_ptr, + const unsigned int * indices, + VectorizedArrayType &dst) + { + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + dst[v] = src_ptr[indices[v]]; + } + + + + // internal helper function for reading data; specialized version where we + // can use a dedicated gather function + template + void + do_vectorized_gather(const Number * src_ptr, + const unsigned int * indices, + VectorizedArray &dst) + { + dst.gather(src_ptr, indices); + } + + + + // internal helper function for reading data; base version of different types + template + void + do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr) + { + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + dst_ptr[v] += src[v]; + } + + + + // internal helper function for reading data; specialized version where we + // can use a dedicated load function + template + void + do_vectorized_add(const VectorizedArray src, Number *dst_ptr) + { + VectorizedArray tmp; + tmp.load(dst_ptr); + (tmp + src).store(dst_ptr); + } + + + + // internal helper function for reading data; base version of different types + template + void + do_vectorized_scatter_add(const VectorizedArrayType src, + const unsigned int * indices, + Number2 * dst_ptr) + { + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + dst_ptr[indices[v]] += src[v]; + } + + + + // internal helper function for reading data; specialized version where we + // can use a dedicated gather function + template + void + do_vectorized_scatter_add(const VectorizedArray src, + const unsigned int * indices, + Number * dst_ptr) + { +#if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3 + for (unsigned int v = 0; v < width; ++v) + dst_ptr[indices[v]] += src[v]; +#else + VectorizedArray tmp; + tmp.gather(dst_ptr, indices); + (tmp + src).scatter(indices, dst_ptr); +#endif + } + + + + template + struct FEFaceEvaluationSelector + { + static void + evaluate(const MatrixFreeFunctions::ShapeInfo &data, + const VectorizedArrayType * values_array, + VectorizedArrayType * values_quad, + VectorizedArrayType * gradients_quad, + VectorizedArrayType * scratch_data, + const bool evaluate_values, + const bool evaluate_gradients, + const unsigned int face_no, + const unsigned int subface_index, + const unsigned int face_orientation, + const Table<2, unsigned int> &orientation_map) + { + constexpr unsigned int static_dofs_per_face = + fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : + numbers::invalid_unsigned_int; + const unsigned int dofs_per_face = + fe_degree > -1 ? static_dofs_per_face : + Utilities::pow(data.fe_degree + 1, dim - 1); + + // we allocate small amounts of data on the stack to signal the compiler + // that this temporary data is only needed for the calculations but the + // final results can be discarded and need not be written back to + // memory. For large sizes or when the dofs per face is not a compile-time + // constant, however, we want to go to the heap in the `scratch_data` + // variable to not risk a stack overflow. + constexpr unsigned int stack_array_size_threshold = 100; + + VectorizedArrayType + temp_data[static_dofs_per_face < stack_array_size_threshold ? + n_components * 2 * static_dofs_per_face : + 1]; + VectorizedArrayType *temp1; + if (static_dofs_per_face < stack_array_size_threshold) + temp1 = &temp_data[0]; + else + temp1 = scratch_data; + + FEFaceNormalEvaluationImpl:: + template interpolate( + data, values_array, temp1, evaluate_gradients, face_no); + + const unsigned int n_q_points_1d_actual = + fe_degree > -1 ? n_q_points_1d : 0; + if (fe_degree > -1 && + subface_index >= GeometryInfo::max_children_per_cell && + data.element_type <= MatrixFreeFunctions::tensor_symmetric) + FEFaceEvaluationImpl< + true, + dim, + fe_degree, + n_q_points_1d_actual, + n_components, + VectorizedArrayType>::evaluate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + 2 * + n_components * + dofs_per_face, + evaluate_values, + evaluate_gradients, + subface_index); + else + FEFaceEvaluationImpl< + false, + dim, + fe_degree, + n_q_points_1d_actual, + n_components, + VectorizedArrayType>::evaluate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + 2 * + n_components * + dofs_per_face, + evaluate_values, + evaluate_gradients, + subface_index); + + if (face_orientation) + adjust_for_face_orientation(face_orientation, + orientation_map, + false, + evaluate_values, + evaluate_gradients, + data.n_q_points_face, + scratch_data, + values_quad, + gradients_quad); + } + + static void + integrate(const MatrixFreeFunctions::ShapeInfo &data, + VectorizedArrayType * values_array, + VectorizedArrayType * values_quad, + VectorizedArrayType * gradients_quad, + VectorizedArrayType * scratch_data, + const bool integrate_values, + const bool integrate_gradients, + const unsigned int face_no, + const unsigned int subface_index, + const unsigned int face_orientation, + const Table<2, unsigned int> &orientation_map) + { + if (face_orientation) + adjust_for_face_orientation(face_orientation, + orientation_map, + true, + integrate_values, + integrate_gradients, + data.n_q_points_face, + scratch_data, + values_quad, + gradients_quad); + + constexpr unsigned int static_dofs_per_face = + fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : + numbers::invalid_unsigned_int; + const unsigned int dofs_per_face = + fe_degree > -1 ? static_dofs_per_face : + Utilities::pow(data.fe_degree + 1, dim - 1); + + constexpr unsigned int stack_array_size_threshold = 100; + + VectorizedArrayType + temp_data[static_dofs_per_face < stack_array_size_threshold ? + n_components * 2 * static_dofs_per_face : + 1]; + VectorizedArrayType *temp1; + if (static_dofs_per_face < stack_array_size_threshold) + temp1 = &temp_data[0]; + else + temp1 = scratch_data; + + const unsigned int n_q_points_1d_actual = + fe_degree > -1 ? n_q_points_1d : 0; + if (fe_degree > -1 && + subface_index >= GeometryInfo::max_children_per_cell && + data.element_type <= MatrixFreeFunctions::tensor_symmetric) + FEFaceEvaluationImpl< + true, + dim, + fe_degree, + n_q_points_1d_actual, + n_components, + VectorizedArrayType>::integrate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + + 2 * n_components * + dofs_per_face, + integrate_values, + integrate_gradients, + subface_index); + else + FEFaceEvaluationImpl< + false, + dim, + fe_degree, + n_q_points_1d_actual, + n_components, + VectorizedArrayType>::integrate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + + 2 * n_components * + dofs_per_face, + integrate_values, + integrate_gradients, + subface_index); + + FEFaceNormalEvaluationImpl:: + template interpolate( + data, temp1, values_array, integrate_gradients, face_no); + } + + static bool + gather_evaluate( + const Number2 * src_ptr, + const MatrixFreeFunctions::ShapeInfo &data, + const MatrixFreeFunctions::DoFInfo & dof_info, + VectorizedArrayType * values_quad, + VectorizedArrayType * gradients_quad, + VectorizedArrayType * scratch_data, + const bool evaluate_values, + const bool evaluate_gradients, + const unsigned int active_fe_index, + const unsigned int first_selected_component, + const unsigned int cell, + const unsigned int face_no, + const unsigned int subface_index, + const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, + const unsigned int face_orientation, + const Table<2, unsigned int> & orientation_map) + { + const unsigned int side = face_no % 2; + + constexpr unsigned int static_dofs_per_component = + fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) : + numbers::invalid_unsigned_int; + constexpr unsigned int static_dofs_per_face = + fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : + numbers::invalid_unsigned_int; + const unsigned int dofs_per_face = + fe_degree > -1 ? static_dofs_per_face : + Utilities::pow(data.fe_degree + 1, dim - 1); + + constexpr unsigned int stack_array_size_threshold = 100; + + VectorizedArrayType + temp_data[static_dofs_per_face < stack_array_size_threshold ? + n_components * 2 * dofs_per_face : + 1]; + VectorizedArrayType *__restrict temp1; + if (static_dofs_per_face < stack_array_size_threshold) + temp1 = &temp_data[0]; + else + temp1 = scratch_data; + + // case 1: contiguous and interleaved indices + if (((evaluate_gradients == false && + data.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) + { + AssertDimension( + dof_info.n_vectorization_lanes_filled[dof_access_index][cell], + VectorizedArrayType::n_array_elements); + const unsigned int dof_index = + dof_info.dof_indices_contiguous + [dof_access_index][cell * VectorizedArrayType::n_array_elements] + + dof_info.component_dof_indices_offset[active_fe_index] + [first_selected_component] * + VectorizedArrayType::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 VectorizedArrayType grad_weight = + data.shape_data_on_face[0][fe_degree + 1 + side]; + AssertDimension(data.face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + const unsigned int *index_array = + &data.face_to_cell_index_hermite(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, data.dofs_per_component_on_cell); + AssertIndexRange(ind2, data.dofs_per_component_on_cell); + for (unsigned int comp = 0; comp < n_components; ++comp) + { + do_vectorized_read( + src_ptr + dof_index + + (ind1 + comp * static_dofs_per_component) * + VectorizedArrayType::n_array_elements, + temp1[i + 2 * comp * dofs_per_face]); + do_vectorized_read( + src_ptr + dof_index + + (ind2 + comp * static_dofs_per_component) * + VectorizedArrayType::n_array_elements, + temp1[dofs_per_face + i + 2 * comp * dofs_per_face]); + 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(data.face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &data.face_to_cell_index_nodal(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) + do_vectorized_read( + src_ptr + dof_index + + (ind + comp * static_dofs_per_component) * + VectorizedArrayType::n_array_elements, + temp1[i + 2 * comp * dofs_per_face]); + } + } + } + + // case 2: contiguous and interleaved indices with fixed stride + else if (((evaluate_gradients == false && + data.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::n_array_elements); + const unsigned int *indices = + &dof_info.dof_indices_contiguous + [dof_access_index][cell * VectorizedArrayType::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 VectorizedArrayType grad_weight = + data.shape_data_on_face[0][fe_degree + 1 + side]; + AssertDimension(data.face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + + const unsigned int *index_array = + &data.face_to_cell_index_hermite(face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = + index_array[2 * i] * VectorizedArrayType::n_array_elements; + const unsigned int ind2 = + index_array[2 * i + 1] * + VectorizedArrayType::n_array_elements; + for (unsigned int comp = 0; comp < n_components; ++comp) + { + do_vectorized_gather( + src_ptr + ind1 + + comp * static_dofs_per_component * + VectorizedArrayType::n_array_elements + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component] * + VectorizedArrayType::n_array_elements, + indices, + temp1[i + 2 * comp * dofs_per_face]); + VectorizedArrayType grad; + do_vectorized_gather( + src_ptr + ind2 + + comp * static_dofs_per_component * + VectorizedArrayType::n_array_elements + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component] * + VectorizedArrayType::n_array_elements, + indices, + grad); + temp1[i + dofs_per_face + 2 * comp * dofs_per_face] = + grad_weight * + (temp1[i + 2 * comp * dofs_per_face] - grad); + } + } + } + else + { + AssertDimension(data.face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &data.face_to_cell_index_nodal(face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind = + index_array[i] * VectorizedArrayType::n_array_elements; + for (unsigned int comp = 0; comp < n_components; ++comp) + do_vectorized_gather( + src_ptr + ind + + comp * static_dofs_per_component * + VectorizedArrayType::n_array_elements + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component] * + VectorizedArrayType::n_array_elements, + indices, + temp1[i + 2 * comp * dofs_per_face]); + } + } + } + + // case 3: contiguous and interleaved indices with mixed stride + else if (((evaluate_gradients == false && + data.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::n_array_elements]; + unsigned int indices[VectorizedArrayType::n_array_elements]; + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; + ++v) + indices[v] = + dof_info.dof_indices_contiguous + [dof_access_index] + [cell * VectorizedArrayType::n_array_elements + v] + + dof_info.component_dof_indices_offset[active_fe_index] + [first_selected_component] * + strides[v]; + const unsigned int nvec = + dof_info.n_vectorization_lanes_filled[dof_access_index][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 VectorizedArrayType grad_weight = + data.shape_data_on_face[0][fe_degree + 1 + side]; + AssertDimension(data.face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + + const unsigned int *index_array = + &data.face_to_cell_index_hermite(face_no, 0); + if (nvec == VectorizedArrayType::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[VectorizedArrayType::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::n_array_elements; + ++v) + ind1[v] = + indices[v] + (comp * static_dofs_per_component + + index_array[2 * i]) * + strides[v]; + unsigned int ind2[VectorizedArrayType::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::n_array_elements; + ++v) + ind2[v] = + indices[v] + (comp * static_dofs_per_component + + index_array[2 * i + 1]) * + strides[v]; + do_vectorized_gather(src_ptr, + ind1, + temp1[i + 2 * comp * dofs_per_face]); + VectorizedArrayType grad; + do_vectorized_gather(src_ptr, ind2, grad); + 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] = VectorizedArrayType(); + 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]; + temp1[i + 2 * comp * dofs_per_face][v] = + src_ptr[ind1]; + const Number grad = src_ptr[ind2]; + 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(data.face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &data.face_to_cell_index_nodal(face_no, 0); + if (nvec == VectorizedArrayType::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[VectorizedArrayType::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::n_array_elements; + ++v) + ind[v] = + indices[v] + + (comp * static_dofs_per_component + index_array[i]) * + strides[v]; + do_vectorized_gather(src_ptr, + ind, + temp1[i + 2 * comp * dofs_per_face]); + } + else + { + for (unsigned int i = 0; i < n_components * dofs_per_face; + ++i) + temp1[i] = VectorizedArrayType(); + 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]; + temp1[i + 2 * comp * dofs_per_face][v] = + src_ptr[ind1]; + } + } + } + } + + // case 4: contiguous indices without interleaving + else if (((evaluate_gradients == false && + data.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::n_array_elements) + { + const unsigned int *indices = + &dof_info.dof_indices_contiguous + [dof_access_index][cell * VectorizedArrayType::n_array_elements]; + if (evaluate_gradients == true && + data.element_type == + MatrixFreeFunctions::tensor_symmetric_hermite) + { + // 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 VectorizedArrayType grad_weight = + data.shape_data_on_face[0][fe_degree + 1 + side]; + AssertDimension(data.face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + + const unsigned int *index_array = + &data.face_to_cell_index_hermite(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]; + for (unsigned int comp = 0; comp < n_components; ++comp) + { + do_vectorized_gather( + src_ptr + ind1 + comp * static_dofs_per_component + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component], + indices, + temp1[i + 2 * comp * dofs_per_face]); + VectorizedArrayType grad; + do_vectorized_gather( + src_ptr + ind2 + comp * static_dofs_per_component + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component], + indices, + grad); + temp1[i + dofs_per_face + 2 * comp * dofs_per_face] = + grad_weight * + (temp1[i + 2 * comp * dofs_per_face] - grad); + } + } + } + else + { + AssertDimension(data.face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &data.face_to_cell_index_nodal(face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + for (unsigned int comp = 0; comp < n_components; ++comp) + { + const unsigned int ind = index_array[i]; + do_vectorized_gather( + src_ptr + ind + comp * static_dofs_per_component + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component], + indices, + temp1[i + comp * 2 * dofs_per_face]); + } + } + } + + // case 5: default vector access + else + { + return false; + } + + if (fe_degree > -1 && + subface_index >= GeometryInfo::max_children_per_cell && + data.element_type <= MatrixFreeFunctions::tensor_symmetric) + FEFaceEvaluationImpl< + true, + dim, + fe_degree, + n_q_points_1d, + n_components, + VectorizedArrayType>::evaluate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + 2 * + n_components * + dofs_per_face, + evaluate_values, + evaluate_gradients, + subface_index); + else + FEFaceEvaluationImpl< + false, + dim, + fe_degree, + n_q_points_1d, + n_components, + VectorizedArrayType>::evaluate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + 2 * + n_components * + dofs_per_face, + evaluate_values, + evaluate_gradients, + subface_index); + + if (face_orientation) + adjust_for_face_orientation(face_orientation, + orientation_map, + false, + evaluate_values, + evaluate_gradients, + data.n_q_points_face, + scratch_data, + values_quad, + gradients_quad); + + return true; + } + + static bool + integrate_scatter( + Number2 * dst_ptr, + const MatrixFreeFunctions::ShapeInfo &data, + const MatrixFreeFunctions::DoFInfo & dof_info, + VectorizedArrayType * values_array, + VectorizedArrayType * values_quad, + VectorizedArrayType * gradients_quad, + VectorizedArrayType * scratch_data, + const bool integrate_values, + const bool integrate_gradients, + const unsigned int active_fe_index, + const unsigned int first_selected_component, + const unsigned int cell, + const unsigned int face_no, + const unsigned int subface_index, + const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, + const unsigned int face_orientation, + const Table<2, unsigned int> & orientation_map) + { + if (face_orientation) + adjust_for_face_orientation(face_orientation, + orientation_map, + true, + integrate_values, + integrate_gradients, + data.n_q_points_face, + scratch_data, + values_quad, + gradients_quad); + + const unsigned int side = face_no % 2; + constexpr unsigned int static_dofs_per_component = + fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) : + numbers::invalid_unsigned_int; + const unsigned int dofs_per_face = + fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::pow(data.fe_degree + 1, dim - 1); + + constexpr unsigned int stack_array_size_threshold = 100; + + VectorizedArrayType temp_data[dofs_per_face < stack_array_size_threshold ? + n_components * 2 * dofs_per_face : + 1]; + VectorizedArrayType *__restrict temp1; + if (dofs_per_face < stack_array_size_threshold) + temp1 = &temp_data[0]; + else + temp1 = scratch_data; + + if (fe_degree > -1 && + subface_index >= GeometryInfo::max_children_per_cell && + data.element_type <= internal::MatrixFreeFunctions::tensor_symmetric) + internal::FEFaceEvaluationImpl< + true, + dim, + fe_degree, + n_q_points_1d, + n_components, + VectorizedArrayType>::integrate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + + 2 * n_components * + dofs_per_face, + integrate_values, + integrate_gradients, + subface_index); + else + internal::FEFaceEvaluationImpl< + false, + dim, + fe_degree, + n_q_points_1d, + n_components, + VectorizedArrayType>::integrate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + + 2 * n_components * + dofs_per_face, + integrate_values, + integrate_gradients, + subface_index); + + // case 1: contiguous and interleaved indices + if (((integrate_gradients == false && + data.nodal_at_cell_boundaries == true) || + (data.element_type == + internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) && + dof_info.index_storage_variants[dof_access_index][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous) + { + AssertDimension( + dof_info.n_vectorization_lanes_filled[dof_access_index][cell], + VectorizedArrayType::n_array_elements); + const unsigned int dof_index = + dof_info.dof_indices_contiguous + [dof_access_index][cell * VectorizedArrayType::n_array_elements] + + dof_info.component_dof_indices_offset[active_fe_index] + [first_selected_component] * + VectorizedArrayType::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 VectorizedArrayType grad_weight = + data.shape_data_on_face[0][fe_degree + 2 - side]; + AssertDimension(data.face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + const unsigned int *index_array = + &data.face_to_cell_index_hermite(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, data.dofs_per_component_on_cell); + AssertIndexRange(ind2, data.dofs_per_component_on_cell); + for (unsigned int comp = 0; comp < n_components; ++comp) + { + VectorizedArrayType val = + temp1[i + 2 * comp * dofs_per_face] - + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + VectorizedArrayType grad = + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + do_vectorized_add( + val, + dst_ptr + dof_index + + (ind1 + comp * static_dofs_per_component) * + VectorizedArrayType::n_array_elements); + do_vectorized_add( + grad, + dst_ptr + dof_index + + (ind2 + comp * static_dofs_per_component) * + VectorizedArrayType::n_array_elements); + } + } + } + else + { + AssertDimension(data.face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &data.face_to_cell_index_nodal(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) + do_vectorized_add( + temp1[i + 2 * comp * dofs_per_face], + dst_ptr + dof_index + + (ind + comp * static_dofs_per_component) * + VectorizedArrayType::n_array_elements); + } + } + } + + // case 2: contiguous and interleaved indices with fixed stride + else if (((integrate_gradients == false && + data.nodal_at_cell_boundaries == true) || + (data.element_type == + internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) && + dof_info.index_storage_variants[dof_access_index][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_strided) + { + AssertDimension( + dof_info.n_vectorization_lanes_filled[dof_access_index][cell], + VectorizedArrayType::n_array_elements); + const unsigned int *indices = + &dof_info.dof_indices_contiguous + [dof_access_index][cell * VectorizedArrayType::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 VectorizedArrayType grad_weight = + data.shape_data_on_face[0][fe_degree + 2 - side]; + AssertDimension(data.face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + + const unsigned int *index_array = + &data.face_to_cell_index_hermite(face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = + index_array[2 * i] * VectorizedArrayType::n_array_elements; + const unsigned int ind2 = + index_array[2 * i + 1] * + VectorizedArrayType::n_array_elements; + for (unsigned int comp = 0; comp < n_components; ++comp) + { + VectorizedArrayType val = + temp1[i + 2 * comp * dofs_per_face] - + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + VectorizedArrayType grad = + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + do_vectorized_scatter_add( + val, + indices, + dst_ptr + ind1 + + comp * static_dofs_per_component * + VectorizedArrayType::n_array_elements + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component] * + VectorizedArrayType::n_array_elements); + do_vectorized_scatter_add( + grad, + indices, + dst_ptr + ind2 + + comp * static_dofs_per_component * + VectorizedArrayType::n_array_elements + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component] * + VectorizedArrayType::n_array_elements); + } + } + } + else + { + AssertDimension(data.face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &data.face_to_cell_index_nodal(face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind = + index_array[i] * VectorizedArrayType::n_array_elements; + for (unsigned int comp = 0; comp < n_components; ++comp) + do_vectorized_scatter_add( + temp1[i + 2 * comp * dofs_per_face], + indices, + dst_ptr + ind + + comp * static_dofs_per_component * + VectorizedArrayType::n_array_elements + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component] * + VectorizedArrayType::n_array_elements); + } + } + } + + // case 3: contiguous and interleaved indices with mixed stride + else if (((integrate_gradients == false && + data.nodal_at_cell_boundaries == true) || + (data.element_type == + internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) && + dof_info.index_storage_variants[dof_access_index][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_mixed_strides) + { + const unsigned int *strides = + &dof_info.dof_indices_interleave_strides + [dof_access_index][cell * VectorizedArrayType::n_array_elements]; + unsigned int indices[VectorizedArrayType::n_array_elements]; + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; + ++v) + indices[v] = + dof_info.dof_indices_contiguous + [dof_access_index] + [cell * VectorizedArrayType::n_array_elements + v] + + dof_info.component_dof_indices_offset[active_fe_index] + [first_selected_component] * + strides[v]; + const unsigned int nvec = + dof_info.n_vectorization_lanes_filled[dof_access_index][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 VectorizedArrayType grad_weight = + data.shape_data_on_face[0][fe_degree + 2 - side]; + AssertDimension(data.face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + + const unsigned int *index_array = + &data.face_to_cell_index_hermite(face_no, 0); + if (nvec == VectorizedArrayType::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[VectorizedArrayType::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::n_array_elements; + ++v) + ind1[v] = + indices[v] + (comp * static_dofs_per_component + + index_array[2 * i]) * + strides[v]; + unsigned int ind2[VectorizedArrayType::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::n_array_elements; + ++v) + ind2[v] = + indices[v] + (comp * static_dofs_per_component + + index_array[2 * i + 1]) * + strides[v]; + VectorizedArrayType val = + temp1[i + 2 * comp * dofs_per_face] - + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + VectorizedArrayType grad = + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + do_vectorized_scatter_add(val, ind1, dst_ptr); + do_vectorized_scatter_add(grad, ind2, dst_ptr); + } + 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]; + dst_ptr[ind1] += val; + dst_ptr[ind2] += grad; + } + } + } + else + { + AssertDimension(data.face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &data.face_to_cell_index_nodal(face_no, 0); + if (nvec == VectorizedArrayType::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[VectorizedArrayType::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::n_array_elements; + ++v) + ind[v] = + indices[v] + + (comp * static_dofs_per_component + index_array[i]) * + strides[v]; + do_vectorized_scatter_add( + temp1[i + 2 * comp * dofs_per_face], ind, dst_ptr); + } + 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]; + dst_ptr[ind1] += + temp1[i + 2 * comp * dofs_per_face][v]; + } + } + } + } + + // case 4: contiguous indices without interleaving + else if (((integrate_gradients == false && + data.nodal_at_cell_boundaries == true) || + (data.element_type == + internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) && + dof_info.index_storage_variants[dof_access_index][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + contiguous && + dof_info.n_vectorization_lanes_filled[dof_access_index][cell] == + VectorizedArrayType::n_array_elements) + { + const unsigned int *indices = + &dof_info.dof_indices_contiguous + [dof_access_index][cell * VectorizedArrayType::n_array_elements]; + + if (integrate_gradients == true && + data.element_type == + internal::MatrixFreeFunctions::tensor_symmetric_hermite) + { + // 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 VectorizedArrayType grad_weight = + data.shape_data_on_face[0][fe_degree + 2 - side]; + AssertDimension(data.face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + const unsigned int *index_array = + &data.face_to_cell_index_hermite(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]; + for (unsigned int comp = 0; comp < n_components; ++comp) + { + VectorizedArrayType val = + temp1[i + 2 * comp * dofs_per_face] - + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + VectorizedArrayType grad = + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + do_vectorized_scatter_add( + val, + indices, + dst_ptr + comp * static_dofs_per_component + ind1 + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component]); + do_vectorized_scatter_add( + grad, + indices, + dst_ptr + comp * static_dofs_per_component + ind2 + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component]); + } + } + } + else + { + AssertDimension(data.face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &data.face_to_cell_index_nodal(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) + do_vectorized_scatter_add( + temp1[i + 2 * comp * dofs_per_face], + indices, + dst_ptr + comp * static_dofs_per_component + ind + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component]); + } + } + } + + // case 5: default vector access, must be handled separately, just do + // the face-normal interpolation + else + { + FEFaceNormalEvaluationImpl:: + template interpolate( + data, temp1, values_array, integrate_gradients, face_no); + return false; + } + + return true; + } + + static void + adjust_for_face_orientation(const unsigned int face_orientation, + const Table<2, unsigned int> &orientation_map, + const bool integrate, + const bool values, + const bool gradients, + const unsigned int n_q_points, + VectorizedArrayType * tmp_values, + VectorizedArrayType * values_quad, + VectorizedArrayType * gradients_quad) + { + Assert(face_orientation, ExcInternalError()); + const unsigned int *orientation = &orientation_map[face_orientation][0]; + for (unsigned int c = 0; c < n_components; ++c) + { + if (values == true) + { + if (integrate) + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[q] = values_quad[c * n_q_points + orientation[q]]; + else + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[orientation[q]] = values_quad[c * n_q_points + q]; + for (unsigned int q = 0; q < n_q_points; ++q) + values_quad[c * n_q_points + q] = tmp_values[q]; + } + if (gradients == true) + for (unsigned int d = 0; d < dim; ++d) + { + if (integrate) + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points + + orientation[q]]; + else + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[orientation[q]] = + gradients_quad[(c * dim + d) * n_q_points + q]; + for (unsigned int q = 0; q < n_q_points; ++q) + gradients_quad[(c * dim + d) * n_q_points + q] = + tmp_values[q]; + } + } + } + }; + + + + template + struct CellwiseInverseMassMatrixImpl + { + static void + apply(const AlignedVector &inverse_shape, + const AlignedVector &inverse_coefficients, + const unsigned int n_desired_components, + const Number * in_array, + Number * out_array) + { + constexpr unsigned int dofs_per_component = + Utilities::pow(fe_degree + 1, dim); + Assert(inverse_coefficients.size() > 0 && + inverse_coefficients.size() % dofs_per_component == 0, + ExcMessage( + "Expected diagonal to be a multiple of scalar dof per cells")); + if (inverse_coefficients.size() != dofs_per_component) + AssertDimension(n_desired_components * dofs_per_component, + inverse_coefficients.size()); + + Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); + + internal::EvaluatorTensorProduct + evaluator(AlignedVector(), + AlignedVector(), + inverse_shape); + + const unsigned int shift_coefficient = + inverse_coefficients.size() > dofs_per_component ? dofs_per_component : + 0; + const Number *inv_coefficient = inverse_coefficients.data(); + for (unsigned int d = 0; d < n_desired_components; ++d) + { + const Number *in = in_array + d * dofs_per_component; + Number * out = out_array + d * dofs_per_component; + // Need to select 'apply' method with hessian slot because values + // assume symmetries that do not exist in the inverse shapes + evaluator.template hessians<0, false, false>(in, out); + if (dim > 1) + { + evaluator.template hessians<1, false, false>(out, out); + + if (dim == 3) + { + evaluator.template hessians<2, false, false>(out, out); + for (unsigned int q = 0; q < dofs_per_component; ++q) + out[q] *= inv_coefficient[q]; + evaluator.template hessians<2, true, false>(out, out); + } + else if (dim == 2) + for (unsigned int q = 0; q < dofs_per_component; ++q) + out[q] *= inv_coefficient[q]; + + evaluator.template hessians<1, true, false>(out, out); + } + else + { + for (unsigned int q = 0; q < dofs_per_component; ++q) + out[q] *= inv_coefficient[q]; + } + evaluator.template hessians<0, true, false>(out, out); + + inv_coefficient += shift_coefficient; + } + } + + static void + transform_from_q_points_to_basis(const AlignedVector &inverse_shape, + const unsigned int n_desired_components, + const Number * in_array, + Number * out_array) + { + constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim); + internal::EvaluatorTensorProduct + evaluator(AlignedVector(), + AlignedVector(), + inverse_shape); + + for (unsigned int d = 0; d < n_desired_components; ++d) + { + const Number *in = in_array + d * dofs_per_cell; + Number * out = out_array + d * dofs_per_cell; + + if (dim == 3) + { + evaluator.template hessians<2, true, false>(in, out); + evaluator.template hessians<1, true, false>(out, out); + evaluator.template hessians<0, true, false>(out, out); + } + if (dim == 2) + { + evaluator.template hessians<1, true, false>(in, out); + evaluator.template hessians<0, true, false>(out, out); + } + if (dim == 1) + evaluator.template hessians<0, true, false>(in, out); + } + } + }; + } // end of namespace internal diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 753fb8fc57..9716922d9d 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -36,6 +36,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -2918,17 +2919,6 @@ public: * points is inaccurate and this value must be used instead. */ const unsigned int n_q_points; - -protected: - /** - * For faces not oriented in the standard way, this method applies - * re-indexing on quadrature points. Called at the end of evaluate() and at - * the beginning of integrate(). - */ - void - adjust_for_face_orientation(const bool integrate, - const bool values, - const bool gradients); }; @@ -3590,654 +3580,6 @@ FEEvaluationBase:: namespace internal { - // below we use type-traits from matrix-free/type_traits.h - - // access to generic const vectors that have operator (). - // FIXME: this is wrong for Trilinos/Petsc MPI vectors - // where we should first do Partitioner::local_to_global() - template ::value, - VectorType>::type * = nullptr> - inline typename VectorType::value_type - vector_access(const VectorType &vec, const unsigned int entry) - { - return vec(entry); - } - - - - // access to generic non-const vectors that have operator (). - // FIXME: this is wrong for Trilinos/Petsc MPI vectors - // where we should first do Partitioner::local_to_global() - template ::value, - VectorType>::type * = nullptr> - inline typename VectorType::value_type & - vector_access(VectorType &vec, const unsigned int entry) - { - return vec(entry); - } - - - - // access to distributed MPI vectors that have a local_element(uint) - // method to access data in local index space, which is what we use in - // DoFInfo and hence in read_dof_values etc. - template ::value, - VectorType>::type * = nullptr> - inline typename VectorType::value_type & - vector_access(VectorType &vec, const unsigned int entry) - { - return vec.local_element(entry); - } - - - - // same for const access - template ::value, - VectorType>::type * = nullptr> - inline typename VectorType::value_type - vector_access(const VectorType &vec, const unsigned int entry) - { - return vec.local_element(entry); - } - - - - template ::value, - VectorType>::type * = nullptr> - inline void - vector_access_add(VectorType & vec, - const unsigned int entry, - const typename VectorType::value_type &val) - { - vec.add_local_element(entry, val); - } - - - - template ::value, - VectorType>::type * = nullptr> - inline void - vector_access_add(VectorType & vec, - const unsigned int entry, - const typename VectorType::value_type &val) - { - vector_access(vec, entry) += val; - } - - - - template ::value, - VectorType>::type * = nullptr> - inline void - vector_access_add_global(VectorType & vec, - const types::global_dof_index entry, - const typename VectorType::value_type &val) - { - vec.add(entry, val); - } - - - - template ::value, - VectorType>::type * = nullptr> - inline void - vector_access_add_global(VectorType & vec, - const types::global_dof_index entry, - const typename VectorType::value_type &val) - { - vec(entry) += val; - } - - - - template ::value, - VectorType>::type * = nullptr> - inline void - vector_access_set(VectorType & vec, - const unsigned int entry, - const typename VectorType::value_type &val) - { - vec.set_local_element(entry, val); - } - - - - template ::value, - VectorType>::type * = nullptr> - inline void - vector_access_set(VectorType & vec, - const unsigned int entry, - const typename VectorType::value_type &val) - { - vector_access(vec, entry) = val; - } - - - - // this is to make sure that the parallel partitioning in VectorType - // is really the same as stored in MatrixFree. - // version below is when has_partitioners_are_compatible == false - // FIXME: this is incorrect for PETSc/Trilinos MPI vectors - template < - typename VectorType, - typename std::enable_if::value, - VectorType>::type * = nullptr> - inline void - check_vector_compatibility( - const VectorType & vec, - const internal::MatrixFreeFunctions::DoFInfo &dof_info) - { - (void)vec; - (void)dof_info; - - AssertDimension(vec.size(), dof_info.vector_partitioner->size()); - } - - - - // same as above for has_partitioners_are_compatible == true - template < - typename VectorType, - typename std::enable_if::value, - VectorType>::type * = nullptr> - inline void - check_vector_compatibility( - const VectorType & vec, - const internal::MatrixFreeFunctions::DoFInfo &dof_info) - { - (void)vec; - (void)dof_info; - Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner), - ExcMessage( - "The parallel layout of the given vector is not " - "compatible with the parallel partitioning in MatrixFree. " - "Use MatrixFree::initialize_dof_vector to get a " - "compatible vector.")); - } - - - - // Below, three classes (VectorReader, VectorSetter, - // VectorDistributorLocalToGlobal) implement the same interface and can be - // used to to read from vector, set elements of a vector and add to elements - // of the vector. - - // 1. A class to read data from vector - template - struct VectorReader - { - template - void - process_dof(const unsigned int index, - const VectorType & vec, - Number & res) const - { - res = vector_access(vec, index); - } - - - - template - void - process_dofs_vectorized(const unsigned int dofs_per_cell, - const unsigned int dof_index, - VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - const Number *vec_ptr = vec.begin() + dof_index; - for (unsigned int i = 0; i < dofs_per_cell; - ++i, vec_ptr += VectorizedArrayType::n_array_elements) - dof_values[i].load(vec_ptr); - } - - - - template - void - process_dofs_vectorized(const unsigned int dofs_per_cell, - const unsigned int dof_index, - const VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - dof_values[i][v] = vector_access( - vec, dof_index + v + i * VectorizedArrayType::n_array_elements); - } - - - - template - void - process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, - const unsigned int * dof_indices, - VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - dealii::vectorized_load_and_transpose(dofs_per_cell, - vec.begin(), - dof_indices, - dof_values); - } - - - - template - void - process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, - const unsigned int * dof_indices, - const VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - for (unsigned int d = 0; d < dofs_per_cell; ++d) - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - dof_values[d][v] = vector_access(vec, dof_indices[v] + d); - } - - - - // variant where VectorType::value_type is the same as Number -> can call - // gather - template - void - process_dof_gather(const unsigned int * indices, - VectorType & vec, - const unsigned int constant_offset, - VectorizedArrayType &res, - std::integral_constant) const - { - res.gather(vec.begin() + constant_offset, indices); - } - - - - // variant where VectorType::value_type is not the same as Number -> must - // manually load the data - template - void - process_dof_gather(const unsigned int * indices, - const VectorType & vec, - const unsigned int constant_offset, - VectorizedArrayType &res, - std::integral_constant) const - { - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - res[v] = vector_access(vec, indices[v] + constant_offset); - } - - - - template - void - process_dof_global(const types::global_dof_index index, - const VectorType & vec, - Number & res) const - { - res = vec(index); - } - - - - void - pre_constraints(const Number &, Number &res) const - { - res = Number(); - } - - - - template - void - process_constraint(const unsigned int index, - const Number weight, - const VectorType & vec, - Number & res) const - { - res += weight * vector_access(vec, index); - } - - - - void - post_constraints(const Number &sum, Number &write_pos) const - { - write_pos = sum; - } - - - - void - process_empty(VectorizedArrayType &res) const - { - res = VectorizedArrayType(); - } - }; - - - - // 2. A class to add values to the vector during - // FEEvaluation::distribute_local_to_global() call - template - struct VectorDistributorLocalToGlobal - { - template - void - process_dof(const unsigned int index, VectorType &vec, Number &res) const - { - vector_access_add(vec, index, res); - } - - - - template - void - process_dofs_vectorized(const unsigned int dofs_per_cell, - const unsigned int dof_index, - VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - Number *vec_ptr = vec.begin() + dof_index; - for (unsigned int i = 0; i < dofs_per_cell; - ++i, vec_ptr += VectorizedArrayType::n_array_elements) - { - VectorizedArrayType tmp; - tmp.load(vec_ptr); - tmp += dof_values[i]; - tmp.store(vec_ptr); - } - } - - - - template - void - process_dofs_vectorized(const unsigned int dofs_per_cell, - const unsigned int dof_index, - VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - vector_access_add(vec, - dof_index + v + - i * VectorizedArrayType::n_array_elements, - dof_values[i][v]); - } - - - - template - void - process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, - const unsigned int * dof_indices, - VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - vectorized_transpose_and_store( - true, dofs_per_cell, dof_values, dof_indices, vec.begin()); - } - - - - template - void - process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, - const unsigned int * dof_indices, - VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - for (unsigned int d = 0; d < dofs_per_cell; ++d) - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - vector_access_add(vec, dof_indices[v] + d, dof_values[d][v]); - } - - - - // variant where VectorType::value_type is the same as Number -> can call - // scatter - template - void - process_dof_gather(const unsigned int * indices, - VectorType & vec, - const unsigned int constant_offset, - VectorizedArrayType &res, - std::integral_constant) const - { -# if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3 - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - vector_access(vec, indices[v] + constant_offset) += res[v]; -# else - // only use gather in case there is also scatter. - VectorizedArrayType tmp; - tmp.gather(vec.begin() + constant_offset, indices); - tmp += res; - tmp.scatter(indices, vec.begin() + constant_offset); -# endif - } - - - - // variant where VectorType::value_type is not the same as Number -> must - // manually append all data - template - void - process_dof_gather(const unsigned int * indices, - VectorType & vec, - const unsigned int constant_offset, - VectorizedArrayType &res, - std::integral_constant) const - { - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - vector_access_add(vec, indices[v] + constant_offset, res[v]); - } - - - - template - void - process_dof_global(const types::global_dof_index index, - VectorType & vec, - Number & res) const - { - vector_access_add_global(vec, index, res); - } - - - - void - pre_constraints(const Number &input, Number &res) const - { - res = input; - } - - - - template - void - process_constraint(const unsigned int index, - const Number weight, - VectorType & vec, - Number & res) const - { - vector_access_add(vec, index, weight * res); - } - - - - void - post_constraints(const Number &, Number &) const - {} - - - - void - process_empty(VectorizedArrayType &) const - {} - }; - - - - // 3. A class to set elements of the vector - template - struct VectorSetter - { - template - void - process_dof(const unsigned int index, VectorType &vec, Number &res) const - { - vector_access(vec, index) = res; - } - - - - template - void - process_dofs_vectorized(const unsigned int dofs_per_cell, - const unsigned int dof_index, - VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - Number *vec_ptr = vec.begin() + dof_index; - for (unsigned int i = 0; i < dofs_per_cell; - ++i, vec_ptr += VectorizedArrayType::n_array_elements) - dof_values[i].store(vec_ptr); - } - - - - template - void - process_dofs_vectorized(const unsigned int dofs_per_cell, - const unsigned int dof_index, - VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - vector_access(vec, - dof_index + v + - i * VectorizedArrayType::n_array_elements) = - dof_values[i][v]; - } - - - - template - void - process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, - const unsigned int * dof_indices, - VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - vectorized_transpose_and_store( - false, dofs_per_cell, dof_values, dof_indices, vec.begin()); - } - - - - template - void - process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, - const unsigned int * dof_indices, - VectorType & vec, - VectorizedArrayType *dof_values, - std::integral_constant) const - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - vector_access(vec, dof_indices[v] + i) = dof_values[i][v]; - } - - - - template - void - process_dof_gather(const unsigned int * indices, - VectorType & vec, - const unsigned int constant_offset, - VectorizedArrayType &res, - std::integral_constant) const - { - res.scatter(indices, vec.begin() + constant_offset); - } - - - - template - void - process_dof_gather(const unsigned int * indices, - VectorType & vec, - const unsigned int constant_offset, - VectorizedArrayType &res, - std::integral_constant) const - { - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - vector_access(vec, indices[v] + constant_offset) = res[v]; - } - - - - template - void - process_dof_global(const types::global_dof_index index, - VectorType & vec, - Number & res) const - { - vec(index) = res; - } - - - - void - pre_constraints(const Number &, Number &) const - {} - - - - template - void - process_constraint(const unsigned int, - const Number, - VectorType &, - Number &) const - {} - - - - void - post_constraints(const Number &, Number &) const - {} - - - - void - process_empty(VectorizedArrayType &) const - {} - }; - - - // allows to select between block vectors and non-block vectors, which // allows to use a unified interface for extracting blocks on block vectors // and doing nothing on usual vectors @@ -8017,80 +7359,25 @@ FEFaceEvaluation -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - numbers::invalid_unsigned_int; - const unsigned int dofs_per_face = - fe_degree > -1 ? static_dofs_per_face : - Utilities::pow(this->data->fe_degree + 1, dim - 1); - - // we allocate small amounts of data on the stack to signal the compiler - // that this temporary data is only needed for the calculations but the - // final results can be discarded and need not be written back to - // memory. For large sizes or when the dofs per face is not a compile-time - // constant, however, we want to go to the heap in the `scratch_data` - // variable to not risk a stack overflow. - constexpr unsigned int stack_array_size_threshold = 100; - - VectorizedArrayType - temp_data[static_dofs_per_face < stack_array_size_threshold ? - n_components * 2 * static_dofs_per_face : - 1]; - VectorizedArrayType *temp1; - if (static_dofs_per_face < stack_array_size_threshold) - temp1 = &temp_data[0]; - else - temp1 = this->scratch_data; - - internal::FEFaceNormalEvaluationImpl< + internal::FEFaceEvaluationSelector< dim, fe_degree, + n_q_points_1d, n_components, - VectorizedArrayType>::template interpolate(*this->data, - values_array, - temp1, - evaluate_gradients, - this->face_no); - - const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0; - if (fe_degree > -1 && - this->subface_index >= GeometryInfo::max_children_per_cell && - this->data->element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) - internal::FEFaceEvaluationImpl< - true, - dim, - fe_degree, - n_q_points_1d_actual, - n_components, - VectorizedArrayType>::evaluate_in_face(*this->data, - temp1, - this->begin_values(), - this->begin_gradients(), - this->scratch_data + - 2 * n_components * dofs_per_face, - evaluate_values, - evaluate_gradients, - this->subface_index); - else - internal::FEFaceEvaluationImpl< - false, - dim, - fe_degree, - n_q_points_1d_actual, - n_components, - VectorizedArrayType>::evaluate_in_face(*this->data, - temp1, - this->begin_values(), - this->begin_gradients(), - this->scratch_data + - 2 * n_components * dofs_per_face, - evaluate_values, - evaluate_gradients, - this->subface_index); - - if (this->face_orientation) - adjust_for_face_orientation(false, evaluate_values, evaluate_gradients); + Number, + VectorizedArrayType>::evaluate(*this->data, + values_array, + this->begin_values(), + this->begin_gradients(), + this->scratch_data, + evaluate_values, + evaluate_gradients, + this->face_no, + this->subface_index, + this->face_orientation, + this->mapping_data + ->descriptor[this->active_fe_index] + .face_orientations); # ifdef DEBUG if (evaluate_values == true) @@ -8146,73 +7433,25 @@ FEFaceEvaluationface_orientation) - adjust_for_face_orientation(true, integrate_values, integrate_gradients); - - constexpr unsigned int static_dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - numbers::invalid_unsigned_int; - const unsigned int dofs_per_face = - fe_degree > -1 ? static_dofs_per_face : - Utilities::pow(this->data->fe_degree + 1, dim - 1); - - constexpr unsigned int stack_array_size_threshold = 100; - - VectorizedArrayType - temp_data[static_dofs_per_face < stack_array_size_threshold ? - n_components * 2 * static_dofs_per_face : - 1]; - VectorizedArrayType *temp1; - if (static_dofs_per_face < stack_array_size_threshold) - temp1 = &temp_data[0]; - else - temp1 = this->scratch_data; - - const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0; - if (fe_degree > -1 && - this->subface_index >= GeometryInfo::max_children_per_cell && - this->data->element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) - internal::FEFaceEvaluationImpl< - true, - dim, - fe_degree, - n_q_points_1d_actual, - n_components, - VectorizedArrayType>::integrate_in_face(*this->data, - temp1, - this->begin_values(), - this->begin_gradients(), - this->scratch_data + - 2 * n_components * - dofs_per_face, - integrate_values, - integrate_gradients, - this->subface_index); - else - internal::FEFaceEvaluationImpl< - false, - dim, - fe_degree, - n_q_points_1d_actual, - n_components, - VectorizedArrayType>::integrate_in_face(*this->data, - temp1, - this->begin_values(), - this->begin_gradients(), - this->scratch_data + - 2 * n_components * - dofs_per_face, - integrate_values, - integrate_gradients, - this->subface_index); - - internal::FEFaceNormalEvaluationImpl:: - template interpolate( - *this->data, temp1, values_array, integrate_gradients, this->face_no); + internal::FEFaceEvaluationSelector< + dim, + fe_degree, + n_q_points_1d, + n_components, + Number, + VectorizedArrayType>::integrate(*this->data, + values_array, + this->begin_values(), + this->begin_gradients(), + this->scratch_data, + integrate_values, + integrate_gradients, + this->face_no, + this->subface_index, + this->face_orientation, + this->mapping_data + ->descriptor[this->active_fe_index] + .face_orientations); } @@ -8235,502 +7474,43 @@ FEFaceEvaluation< const bool evaluate_values, const bool evaluate_gradients) { - const unsigned int side = this->face_no % 2; - - constexpr unsigned int static_dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - numbers::invalid_unsigned_int; - const unsigned int dofs_per_face = - fe_degree > -1 ? static_dofs_per_face : - Utilities::pow(this->data->fe_degree + 1, dim - 1); - - constexpr unsigned int stack_array_size_threshold = 100; - - VectorizedArrayType - temp_data[static_dofs_per_face < stack_array_size_threshold ? - n_components_ * 2 * dofs_per_face : - 1]; - VectorizedArrayType *__restrict temp1; - if (static_dofs_per_face < stack_array_size_threshold) - temp1 = &temp_data[0]; - else - temp1 = this->scratch_data; - - internal::VectorReader reader; - std::integral_constant::value> - vector_selector; - - // 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)) && - 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], - VectorizedArrayType::n_array_elements); - const unsigned int dof_index = - this->dof_info - ->dof_indices_contiguous[this->dof_access_index] - [this->cell * - VectorizedArrayType::n_array_elements] + - this->dof_info - ->component_dof_indices_offset[this->active_fe_index] - [this->first_selected_component] * - VectorizedArrayType::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 VectorizedArrayType 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) * - VectorizedArrayType::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) * - VectorizedArrayType::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) * - VectorizedArrayType::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], - VectorizedArrayType::n_array_elements); - const unsigned int *indices = - &this->dof_info - ->dof_indices_contiguous[this->dof_access_index] - [this->cell * - VectorizedArrayType::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 VectorizedArrayType 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] * VectorizedArrayType::n_array_elements; - const unsigned int ind2 = - index_array[2 * i + 1] * VectorizedArrayType::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 * - VectorizedArrayType::n_array_elements + - this->dof_info->component_dof_indices_offset - [this->active_fe_index] - [this->first_selected_component] * - VectorizedArrayType::n_array_elements, - temp1[i + 2 * comp * dofs_per_face], - vector_selector); - VectorizedArrayType grad; - reader.process_dof_gather( - indices, - input_vector, - ind2 + - comp * static_dofs_per_component * - VectorizedArrayType::n_array_elements + - this->dof_info->component_dof_indices_offset - [this->active_fe_index] - [this->first_selected_component] * - VectorizedArrayType::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] * VectorizedArrayType::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 * - VectorizedArrayType::n_array_elements + - this->dof_info->component_dof_indices_offset - [this->active_fe_index] - [this->first_selected_component] * - VectorizedArrayType::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 * VectorizedArrayType::n_array_elements]; - unsigned int indices[VectorizedArrayType::n_array_elements]; - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - indices[v] = - this->dof_info->dof_indices_contiguous - [this->dof_access_index] - [this->cell * VectorizedArrayType::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 VectorizedArrayType 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 == VectorizedArrayType::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[VectorizedArrayType::n_array_elements]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::n_array_elements; - ++v) - ind1[v] = indices[v] + (comp * static_dofs_per_component + - index_array[2 * i]) * - strides[v]; - unsigned int ind2[VectorizedArrayType::n_array_elements]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::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); - VectorizedArrayType 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] = VectorizedArrayType(); - 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(input_vector), - temp1[i + 2 * comp * dofs_per_face][v]); - Number grad; - reader.process_dof(ind2, - const_cast(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 == VectorizedArrayType::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[VectorizedArrayType::n_array_elements]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::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] = VectorizedArrayType(); - 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(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] == - VectorizedArrayType::n_array_elements) - { - const unsigned int *indices = - &this->dof_info - ->dof_indices_contiguous[this->dof_access_index] - [this->cell * - VectorizedArrayType::n_array_elements]; - if (evaluate_gradients == true && - this->data->element_type == - internal::MatrixFreeFunctions::tensor_symmetric_hermite) - { - // 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 VectorizedArrayType 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]; - for (unsigned int comp = 0; comp < n_components_; ++comp) - { - reader.process_dof_gather( - indices, - input_vector, - ind1 + comp * static_dofs_per_component + - this->dof_info->component_dof_indices_offset - [this->active_fe_index][this->first_selected_component], - temp1[i + 2 * comp * dofs_per_face], - vector_selector); - VectorizedArrayType grad; - reader.process_dof_gather( - indices, - input_vector, - ind2 + comp * static_dofs_per_component + - this->dof_info->component_dof_indices_offset - [this->active_fe_index][this->first_selected_component], - 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) - for (unsigned int comp = 0; comp < n_components_; ++comp) - { - const unsigned int ind = index_array[i]; - reader.process_dof_gather( - indices, - input_vector, - ind + comp * static_dofs_per_component + - this->dof_info->component_dof_indices_offset - [this->active_fe_index][this->first_selected_component], - temp1[i + comp * 2 * dofs_per_face], - vector_selector); - } - } - } - - // case 5: default vector access - else + static_assert(internal::has_begin::value && + (std::is_same().begin()), + double *>::value || + std::is_same().begin()), + float *>::value), + "This function requires a vector type with begin() function " + "evaluating to a pointer to basic number (float,double). " + "Use read_dof_values() followed by evaluate() instead."); + + if (!internal::FEFaceEvaluationSelector:: + gather_evaluate(input_vector.begin(), + *this->data, + *this->dof_info, + this->begin_values(), + this->begin_gradients(), + this->scratch_data, + evaluate_values, + evaluate_gradients, + this->active_fe_index, + this->first_selected_component, + this->cell, + this->face_no, + this->subface_index, + this->dof_access_index, + this->face_orientation, + this->mapping_data->descriptor[this->active_fe_index] + .face_orientations)) { this->read_dof_values(input_vector); - internal::FEFaceNormalEvaluationImpl:: - template interpolate(*this->data, - this->values_dofs[0], - temp1, - evaluate_gradients, - this->face_no); + this->evaluate(evaluate_values, evaluate_gradients); } - if (fe_degree > -1 && - this->subface_index >= GeometryInfo::max_children_per_cell && - this->data->element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) - internal::FEFaceEvaluationImpl< - true, - dim, - fe_degree, - n_q_points_1d, - n_components_, - VectorizedArrayType>::evaluate_in_face(*this->data, - temp1, - this->values_quad[0], - this->gradients_quad[0][0], - this->scratch_data + - 2 * n_components_ * - dofs_per_face, - evaluate_values, - evaluate_gradients, - this->subface_index); - else - internal::FEFaceEvaluationImpl< - false, - dim, - fe_degree, - n_q_points_1d, - n_components_, - VectorizedArrayType>::evaluate_in_face(*this->data, - temp1, - this->values_quad[0], - this->gradients_quad[0][0], - this->scratch_data + - 2 * n_components_ * - dofs_per_face, - evaluate_values, - evaluate_gradients, - this->subface_index); - - if (this->face_orientation) - adjust_for_face_orientation(false, evaluate_values, evaluate_gradients); - # ifdef DEBUG if (evaluate_values == true) this->values_quad_initialized = true; @@ -8759,559 +7539,51 @@ FEFaceEvaluation< const bool integrate_gradients, VectorType &destination) { - const unsigned int side = this->face_no % 2; - const unsigned int dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - Utilities::pow(this->data->fe_degree + 1, dim - 1); - - constexpr unsigned int stack_array_size_threshold = 100; - - VectorizedArrayType temp_data[dofs_per_face < stack_array_size_threshold ? - n_components_ * 2 * dofs_per_face : - 1]; - VectorizedArrayType *__restrict temp1; - if (dofs_per_face < stack_array_size_threshold) - temp1 = &temp_data[0]; - else - temp1 = this->scratch_data; - - if (this->face_orientation) - adjust_for_face_orientation(true, integrate_values, integrate_gradients); - if (fe_degree > -1 && - this->subface_index >= GeometryInfo::max_children_per_cell && - this->data->element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) - internal::FEFaceEvaluationImpl< - true, - dim, - fe_degree, - n_q_points_1d, - n_components_, - VectorizedArrayType>::integrate_in_face(*this->data, - temp1, - this->values_quad[0], - this->gradients_quad[0][0], - this->scratch_data + - 2 * n_components_ * - dofs_per_face, - integrate_values, - integrate_gradients, - this->subface_index); - else - internal::FEFaceEvaluationImpl< - false, - dim, - fe_degree, - n_q_points_1d, - n_components_, - VectorizedArrayType>::integrate_in_face(*this->data, - temp1, - this->values_quad[0], - this->gradients_quad[0][0], - this->scratch_data + - 2 * n_components_ * - dofs_per_face, - integrate_values, - integrate_gradients, - this->subface_index); - -# ifdef DEBUG - this->dof_values_initialized = true; -# endif - - internal::VectorDistributorLocalToGlobal writer; - std::integral_constant::value> - vector_selector; - - // 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)) && - 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], - VectorizedArrayType::n_array_elements); - const unsigned int dof_index = - this->dof_info - ->dof_indices_contiguous[this->dof_access_index] - [this->cell * - VectorizedArrayType::n_array_elements] + - this->dof_info - ->component_dof_indices_offset[this->active_fe_index] - [this->first_selected_component] * - VectorizedArrayType::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 VectorizedArrayType 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) - { - VectorizedArrayType val = - temp1[i + 2 * comp * dofs_per_face] - - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - VectorizedArrayType 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) * - VectorizedArrayType::n_array_elements, - destination, - &val, - vector_selector); - writer.process_dofs_vectorized( - 1, - dof_index + (ind2 + comp * static_dofs_per_component) * - VectorizedArrayType::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) * - VectorizedArrayType::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], - VectorizedArrayType::n_array_elements); - const unsigned int *indices = - &this->dof_info - ->dof_indices_contiguous[this->dof_access_index] - [this->cell * - VectorizedArrayType::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 VectorizedArrayType 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] * VectorizedArrayType::n_array_elements; - const unsigned int ind2 = - index_array[2 * i + 1] * VectorizedArrayType::n_array_elements; - for (unsigned int comp = 0; comp < n_components_; ++comp) - { - VectorizedArrayType val = - temp1[i + 2 * comp * dofs_per_face] - - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - VectorizedArrayType 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 * - VectorizedArrayType::n_array_elements + - this->dof_info->component_dof_indices_offset - [this->active_fe_index] - [this->first_selected_component] * - VectorizedArrayType::n_array_elements, - val, - vector_selector); - writer.process_dof_gather( - indices, - destination, - ind2 + - comp * static_dofs_per_component * - VectorizedArrayType::n_array_elements + - this->dof_info->component_dof_indices_offset - [this->active_fe_index] - [this->first_selected_component] * - VectorizedArrayType::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] * VectorizedArrayType::n_array_elements; - for (unsigned int comp = 0; comp < n_components_; ++comp) - writer.process_dof_gather( - indices, - destination, - ind + - comp * static_dofs_per_component * - VectorizedArrayType::n_array_elements + - this->dof_info->component_dof_indices_offset - [this->active_fe_index] - [this->first_selected_component] * - VectorizedArrayType::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 * VectorizedArrayType::n_array_elements]; - unsigned int indices[VectorizedArrayType::n_array_elements]; - for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) - indices[v] = - this->dof_info->dof_indices_contiguous - [this->dof_access_index] - [this->cell * VectorizedArrayType::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 VectorizedArrayType 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 == VectorizedArrayType::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[VectorizedArrayType::n_array_elements]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::n_array_elements; - ++v) - ind1[v] = indices[v] + (comp * static_dofs_per_component + - index_array[2 * i]) * - strides[v]; - unsigned int ind2[VectorizedArrayType::n_array_elements]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::n_array_elements; - ++v) - ind2[v] = indices[v] + (comp * static_dofs_per_component + - index_array[2 * i + 1]) * - strides[v]; - VectorizedArrayType val = - temp1[i + 2 * comp * dofs_per_face] - - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - VectorizedArrayType 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 == VectorizedArrayType::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[VectorizedArrayType::n_array_elements]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::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] == - VectorizedArrayType::n_array_elements) - { - const unsigned int *indices = - &this->dof_info - ->dof_indices_contiguous[this->dof_access_index] - [this->cell * - VectorizedArrayType::n_array_elements]; - - if (integrate_gradients == true && - this->data->element_type == - internal::MatrixFreeFunctions::tensor_symmetric_hermite) - { - // 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 VectorizedArrayType 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]; - for (unsigned int comp = 0; comp < n_components_; ++comp) - { - VectorizedArrayType val = - temp1[i + 2 * comp * dofs_per_face] - - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - VectorizedArrayType grad = - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - writer.process_dof_gather( - indices, - destination, - comp * static_dofs_per_component + ind1 + - this->dof_info->component_dof_indices_offset - [this->active_fe_index][this->first_selected_component], - val, - vector_selector); - writer.process_dof_gather( - indices, - destination, - comp * static_dofs_per_component + ind2 + - this->dof_info->component_dof_indices_offset - [this->active_fe_index][this->first_selected_component], - 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_dof_gather( - indices, - destination, - comp * static_dofs_per_component + ind + - this->dof_info->component_dof_indices_offset - [this->active_fe_index][this->first_selected_component], - temp1[i + 2 * comp * dofs_per_face], - vector_selector); - } - } - } - - // case 5: default vector access - else + static_assert(internal::has_begin::value && + (std::is_same().begin()), + double *>::value || + std::is_same().begin()), + float *>::value), + "This function requires a vector type with begin() function " + "evaluating to a pointer to basic number (float,double). " + "Use integrate() followed by distribute_local_to_global() " + "instead."); + + if (!internal::FEFaceEvaluationSelector:: + integrate_scatter(destination.begin(), + *this->data, + *this->dof_info, + this->begin_dof_values(), + this->begin_values(), + this->begin_gradients(), + this->scratch_data, + integrate_values, + integrate_gradients, + this->active_fe_index, + this->first_selected_component, + this->cell, + this->face_no, + this->subface_index, + this->dof_access_index, + this->face_orientation, + this->mapping_data->descriptor[this->active_fe_index] + .face_orientations)) { - internal::FEFaceNormalEvaluationImpl:: - template interpolate(*this->data, - temp1, - this->values_dofs[0], - integrate_gradients, - this->face_no); + // if we arrive here, writing into the destination vector did not succeed + // because some of the assumptions in integrate_scatter were not + // fulfilled (e.g. an element or degree that does not support direct + // writing), so we must do it here this->distribute_local_to_global(destination); } } -template -inline void -FEFaceEvaluation< - dim, - fe_degree, - n_q_points_1d, - n_components, - Number, - VectorizedArrayType>::adjust_for_face_orientation(const bool integrate, - const bool values, - const bool gradients) -{ - VectorizedArrayType *tmp_values = this->scratch_data; - const unsigned int * orientations = - &this->mapping_data->descriptor[this->active_fe_index] - .face_orientations[this->face_orientation][0]; - for (unsigned int c = 0; c < n_components; ++c) - { - if (values == true) - { - if (integrate) - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[q] = this->values_quad[c][orientations[q]]; - else - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[orientations[q]] = this->values_quad[c][q]; - for (unsigned int q = 0; q < n_q_points; ++q) - this->values_quad[c][q] = tmp_values[q]; - } - if (gradients == true) - for (unsigned int d = 0; d < dim; ++d) - { - if (integrate) - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[q] = this->gradients_quad[c][d][orientations[q]]; - else - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[orientations[q]] = this->gradients_quad[c][d][q]; - for (unsigned int q = 0; q < n_q_points; ++q) - this->gradients_quad[c][d][q] = tmp_values[q]; - } - } -} - - - template 0 && - inverse_coefficients.size() % dofs_per_component == 0, - ExcMessage( - "Expected diagonal to be a multiple of scalar dof per cells")); - if (inverse_coefficients.size() != dofs_per_component) - AssertDimension(n_actual_components * dofs_per_component, - inverse_coefficients.size()); - - Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); - - internal::EvaluatorTensorProduct - evaluator(inverse_shape, inverse_shape, inverse_shape); - - const unsigned int shift_coefficient = - inverse_coefficients.size() > dofs_per_component ? dofs_per_component : 0; - const VectorizedArrayType *inv_coefficient = inverse_coefficients.data(); - VectorizedArrayType temp_data_field[dofs_per_component]; - for (unsigned int d = 0; d < n_actual_components; ++d) - { - const VectorizedArrayType *in = in_array + d * dofs_per_component; - VectorizedArrayType * out = out_array + d * dofs_per_component; - // Need to select 'apply' method with hessian slot because values - // assume symmetries that do not exist in the inverse shapes - evaluator.template hessians<0, false, false>(in, temp_data_field); - if (dim > 1) - { - evaluator.template hessians<1, false, false>(temp_data_field, out); - - if (dim == 3) - { - evaluator.template hessians<2, false, false>(out, - temp_data_field); - for (unsigned int q = 0; q < dofs_per_component; ++q) - temp_data_field[q] *= inv_coefficient[q]; - evaluator.template hessians<2, true, false>(temp_data_field, - out); - } - else if (dim == 2) - for (unsigned int q = 0; q < dofs_per_component; ++q) - out[q] *= inv_coefficient[q]; - - evaluator.template hessians<1, true, false>(out, temp_data_field); - } - else - { - for (unsigned int q = 0; q < dofs_per_component; ++q) - temp_data_field[q] *= inv_coefficient[q]; - } - evaluator.template hessians<0, true, false>(temp_data_field, out); - - inv_coefficient += shift_coefficient; - } + internal::CellwiseInverseMassMatrixImpl< + dim, + fe_degree, + n_components, + VectorizedArrayType>::apply(inverse_shape, + inverse_coefficients, + n_actual_components, + in_array, + out_array); } @@ -1091,35 +1042,14 @@ namespace MatrixFreeOperators const VectorizedArrayType *in_array, VectorizedArrayType * out_array) const { - constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim); - internal::EvaluatorTensorProduct - evaluator(AlignedVector(), - AlignedVector(), - inverse_shape); - - for (unsigned int d = 0; d < n_actual_components; ++d) - { - const VectorizedArrayType *in = in_array + d * dofs_per_cell; - VectorizedArrayType * out = out_array + d * dofs_per_cell; - - if (dim == 3) - { - evaluator.template hessians<2, true, false>(in, out); - evaluator.template hessians<1, true, false>(out, out); - evaluator.template hessians<0, true, false>(out, out); - } - if (dim == 2) - { - evaluator.template hessians<1, true, false>(in, out); - evaluator.template hessians<0, true, false>(out, out); - } - if (dim == 1) - evaluator.template hessians<0, true, false>(in, out); - } + internal::CellwiseInverseMassMatrixImpl:: + transform_from_q_points_to_basis(inverse_shape, + n_actual_components, + in_array, + out_array); } diff --git a/include/deal.II/matrix_free/vector_access_internal.h b/include/deal.II/matrix_free/vector_access_internal.h new file mode 100644 index 0000000000..5275f90aa5 --- /dev/null +++ b/include/deal.II/matrix_free/vector_access_internal.h @@ -0,0 +1,683 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_vector_access_internal_h +#define dealii_matrix_free_vector_access_internal_h + +#include + +#include + +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace internal +{ + // below we use type-traits from matrix-free/type_traits.h + + // access to generic const vectors that have operator (). + // FIXME: this is wrong for Trilinos/Petsc MPI vectors + // where we should first do Partitioner::local_to_global() + template ::value, + VectorType>::type * = nullptr> + inline typename VectorType::value_type + vector_access(const VectorType &vec, const unsigned int entry) + { + return vec(entry); + } + + + + // access to generic non-const vectors that have operator (). + // FIXME: this is wrong for Trilinos/Petsc MPI vectors + // where we should first do Partitioner::local_to_global() + template ::value, + VectorType>::type * = nullptr> + inline typename VectorType::value_type & + vector_access(VectorType &vec, const unsigned int entry) + { + return vec(entry); + } + + + + // access to distributed MPI vectors that have a local_element(uint) + // method to access data in local index space, which is what we use in + // DoFInfo and hence in read_dof_values etc. + template ::value, + VectorType>::type * = nullptr> + inline typename VectorType::value_type & + vector_access(VectorType &vec, const unsigned int entry) + { + return vec.local_element(entry); + } + + + + // same for const access + template ::value, + VectorType>::type * = nullptr> + inline typename VectorType::value_type + vector_access(const VectorType &vec, const unsigned int entry) + { + return vec.local_element(entry); + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_add(VectorType & vec, + const unsigned int entry, + const typename VectorType::value_type &val) + { + vec.add_local_element(entry, val); + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_add(VectorType & vec, + const unsigned int entry, + const typename VectorType::value_type &val) + { + vector_access(vec, entry) += val; + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_add_global(VectorType & vec, + const types::global_dof_index entry, + const typename VectorType::value_type &val) + { + vec.add(entry, val); + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_add_global(VectorType & vec, + const types::global_dof_index entry, + const typename VectorType::value_type &val) + { + vec(entry) += val; + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_set(VectorType & vec, + const unsigned int entry, + const typename VectorType::value_type &val) + { + vec.set_local_element(entry, val); + } + + + + template ::value, + VectorType>::type * = nullptr> + inline void + vector_access_set(VectorType & vec, + const unsigned int entry, + const typename VectorType::value_type &val) + { + vector_access(vec, entry) = val; + } + + + + // this is to make sure that the parallel partitioning in VectorType + // is really the same as stored in MatrixFree. + // version below is when has_partitioners_are_compatible == false + // FIXME: this is incorrect for PETSc/Trilinos MPI vectors + template < + typename VectorType, + typename std::enable_if::value, + VectorType>::type * = nullptr> + inline void + check_vector_compatibility( + const VectorType & vec, + const internal::MatrixFreeFunctions::DoFInfo &dof_info) + { + (void)vec; + (void)dof_info; + + AssertDimension(vec.size(), dof_info.vector_partitioner->size()); + } + + + + // same as above for has_partitioners_are_compatible == true + template < + typename VectorType, + typename std::enable_if::value, + VectorType>::type * = nullptr> + inline void + check_vector_compatibility( + const VectorType & vec, + const internal::MatrixFreeFunctions::DoFInfo &dof_info) + { + (void)vec; + (void)dof_info; + Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner), + ExcMessage( + "The parallel layout of the given vector is not " + "compatible with the parallel partitioning in MatrixFree. " + "Use MatrixFree::initialize_dof_vector to get a " + "compatible vector.")); + } + + + + // Below, three classes (VectorReader, VectorSetter, + // VectorDistributorLocalToGlobal) implement the same interface and can be + // used to to read from vector, set elements of a vector and add to elements + // of the vector. + + // 1. A class to read data from vector + template + struct VectorReader + { + template + void + process_dof(const unsigned int index, + const VectorType & vec, + Number & res) const + { + res = vector_access(vec, index); + } + + + + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + const Number *vec_ptr = vec.begin() + dof_index; + for (unsigned int i = 0; i < dofs_per_cell; + ++i, vec_ptr += VectorizedArrayType::n_array_elements) + dof_values[i].load(vec_ptr); + } + + + + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + const VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + dof_values[i][v] = vector_access( + vec, dof_index + v + i * VectorizedArrayType::n_array_elements); + } + + + + template + void + process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, + const unsigned int * dof_indices, + VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + dealii::vectorized_load_and_transpose(dofs_per_cell, + vec.begin(), + dof_indices, + dof_values); + } + + + + template + void + process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, + const unsigned int * dof_indices, + const VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + for (unsigned int d = 0; d < dofs_per_cell; ++d) + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + dof_values[d][v] = vector_access(vec, dof_indices[v] + d); + } + + + + // variant where VectorType::value_type is the same as Number -> can call + // gather + template + void + process_dof_gather(const unsigned int * indices, + VectorType & vec, + const unsigned int constant_offset, + VectorizedArrayType &res, + std::integral_constant) const + { + res.gather(vec.begin() + constant_offset, indices); + } + + + + // variant where VectorType::value_type is not the same as Number -> must + // manually load the data + template + void + process_dof_gather(const unsigned int * indices, + const VectorType & vec, + const unsigned int constant_offset, + VectorizedArrayType &res, + std::integral_constant) const + { + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + res[v] = vector_access(vec, indices[v] + constant_offset); + } + + + + template + void + process_dof_global(const types::global_dof_index index, + const VectorType & vec, + Number & res) const + { + res = vec(index); + } + + + + void + pre_constraints(const Number &, Number &res) const + { + res = Number(); + } + + + + template + void + process_constraint(const unsigned int index, + const Number weight, + const VectorType & vec, + Number & res) const + { + res += weight * vector_access(vec, index); + } + + + + void + post_constraints(const Number &sum, Number &write_pos) const + { + write_pos = sum; + } + + + + void + process_empty(VectorizedArrayType &res) const + { + res = VectorizedArrayType(); + } + }; + + + + // 2. A class to add values to the vector during + // FEEvaluation::distribute_local_to_global() call + template + struct VectorDistributorLocalToGlobal + { + template + void + process_dof(const unsigned int index, VectorType &vec, Number &res) const + { + vector_access_add(vec, index, res); + } + + + + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + Number *vec_ptr = vec.begin() + dof_index; + for (unsigned int i = 0; i < dofs_per_cell; + ++i, vec_ptr += VectorizedArrayType::n_array_elements) + { + VectorizedArrayType tmp; + tmp.load(vec_ptr); + tmp += dof_values[i]; + tmp.store(vec_ptr); + } + } + + + + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + vector_access_add(vec, + dof_index + v + + i * VectorizedArrayType::n_array_elements, + dof_values[i][v]); + } + + + + template + void + process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, + const unsigned int * dof_indices, + VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + vectorized_transpose_and_store( + true, dofs_per_cell, dof_values, dof_indices, vec.begin()); + } + + + + template + void + process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, + const unsigned int * dof_indices, + VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + for (unsigned int d = 0; d < dofs_per_cell; ++d) + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + vector_access_add(vec, dof_indices[v] + d, dof_values[d][v]); + } + + + + // variant where VectorType::value_type is the same as Number -> can call + // scatter + template + void + process_dof_gather(const unsigned int * indices, + VectorType & vec, + const unsigned int constant_offset, + VectorizedArrayType &res, + std::integral_constant) const + { +#if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3 + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + vector_access(vec, indices[v] + constant_offset) += res[v]; +#else + // only use gather in case there is also scatter. + VectorizedArrayType tmp; + tmp.gather(vec.begin() + constant_offset, indices); + tmp += res; + tmp.scatter(indices, vec.begin() + constant_offset); +#endif + } + + + + // variant where VectorType::value_type is not the same as Number -> must + // manually append all data + template + void + process_dof_gather(const unsigned int * indices, + VectorType & vec, + const unsigned int constant_offset, + VectorizedArrayType &res, + std::integral_constant) const + { + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + vector_access_add(vec, indices[v] + constant_offset, res[v]); + } + + + + template + void + process_dof_global(const types::global_dof_index index, + VectorType & vec, + Number & res) const + { + vector_access_add_global(vec, index, res); + } + + + + void + pre_constraints(const Number &input, Number &res) const + { + res = input; + } + + + + template + void + process_constraint(const unsigned int index, + const Number weight, + VectorType & vec, + Number & res) const + { + vector_access_add(vec, index, weight * res); + } + + + + void + post_constraints(const Number &, Number &) const + {} + + + + void + process_empty(VectorizedArrayType &) const + {} + }; + + + + // 3. A class to set elements of the vector + template + struct VectorSetter + { + template + void + process_dof(const unsigned int index, VectorType &vec, Number &res) const + { + vector_access(vec, index) = res; + } + + + + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + Number *vec_ptr = vec.begin() + dof_index; + for (unsigned int i = 0; i < dofs_per_cell; + ++i, vec_ptr += VectorizedArrayType::n_array_elements) + dof_values[i].store(vec_ptr); + } + + + + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + vector_access(vec, + dof_index + v + + i * VectorizedArrayType::n_array_elements) = + dof_values[i][v]; + } + + + + template + void + process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, + const unsigned int * dof_indices, + VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + vectorized_transpose_and_store( + false, dofs_per_cell, dof_values, dof_indices, vec.begin()); + } + + + + template + void + process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, + const unsigned int * dof_indices, + VectorType & vec, + VectorizedArrayType *dof_values, + std::integral_constant) const + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + vector_access(vec, dof_indices[v] + i) = dof_values[i][v]; + } + + + + template + void + process_dof_gather(const unsigned int * indices, + VectorType & vec, + const unsigned int constant_offset, + VectorizedArrayType &res, + std::integral_constant) const + { + res.scatter(indices, vec.begin() + constant_offset); + } + + + + template + void + process_dof_gather(const unsigned int * indices, + VectorType & vec, + const unsigned int constant_offset, + VectorizedArrayType &res, + std::integral_constant) const + { + for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v) + vector_access(vec, indices[v] + constant_offset) = res[v]; + } + + + + template + void + process_dof_global(const types::global_dof_index index, + VectorType & vec, + Number & res) const + { + vec(index) = res; + } + + + + void + pre_constraints(const Number &, Number &) const + {} + + + + template + void + process_constraint(const unsigned int, + const Number, + VectorType &, + Number &) const + {} + + + + void + post_constraints(const Number &, Number &) const + {} + + + + void + process_empty(VectorizedArrayType &) const + {} + }; +} // namespace internal + + +DEAL_II_NAMESPACE_CLOSE + +#endif