From: Peter Munch Date: Sun, 9 Aug 2020 14:01:19 +0000 (+0200) Subject: Merge gather and adjust_for_face_orientation X-Git-Tag: v9.3.0-rc1~1199^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10809%2Fhead;p=dealii.git Merge gather and adjust_for_face_orientation --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 1a101a63d0..0cadb08e87 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -2355,7 +2355,9 @@ namespace internal { (void)subface_index; - if (integrate && (face_orientation > 0)) + if (integrate && + (face_orientation > 0 && + subface_index < GeometryInfo::max_children_per_cell)) adjust_for_face_orientation(face_orientation, orientation_map, true, @@ -2396,6 +2398,18 @@ namespace internal if (integrate) function_0(temp1, dofs_per_face); + const unsigned int dummy = 0; + const unsigned int *orientation = + (data.data.front().nodal_at_cell_boundaries == true) ? + &data.face_orientations[face_orientation][0] : + &dummy; + const auto reorientate = [&](const unsigned int i) { + return (dim < 3 || face_orientation == 0 || + subface_index < GeometryInfo::max_children_per_cell) ? + i : + orientation[i]; + }; + // case 1: contiguous and interleaved indices if (((do_gradients == false && data.data.front().nodal_at_cell_boundaries == true) || @@ -2434,10 +2448,11 @@ namespace internal 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); + const unsigned int i_ = reorientate(i); for (unsigned int comp = 0; comp < n_components; ++comp) function_1a( - temp1[i + 2 * comp * dofs_per_face], - temp1[i + dofs_per_face + 2 * comp * dofs_per_face], + temp1[i_ + 2 * comp * dofs_per_face], + temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face], global_vector_ptr + dof_index + (ind1 + comp * static_dofs_per_component) * VectorizedArrayType::size(), @@ -2455,9 +2470,10 @@ namespace internal &data.face_to_cell_index_nodal(face_no, 0); for (unsigned int i = 0; i < dofs_per_face; ++i) { + const unsigned int i_ = reorientate(i); const unsigned int ind = index_array[i]; for (unsigned int comp = 0; comp < n_components; ++comp) - function_1b(temp1[i + 2 * comp * dofs_per_face], + function_1b(temp1[i_ + 2 * comp * dofs_per_face], global_vector_ptr + dof_index + (ind + comp * static_dofs_per_component) * VectorizedArrayType::size()); @@ -2496,14 +2512,16 @@ namespace internal &data.face_to_cell_index_hermite(face_no, 0); for (unsigned int i = 0; i < dofs_per_face; ++i) { + const unsigned int i_ = reorientate(i); + const unsigned int ind1 = index_array[2 * i] * VectorizedArrayType::size(); const unsigned int ind2 = index_array[2 * i + 1] * VectorizedArrayType::size(); for (unsigned int comp = 0; comp < n_components; ++comp) function_2a( - temp1[i + 2 * comp * dofs_per_face], - temp1[i + dofs_per_face + 2 * comp * dofs_per_face], + temp1[i_ + 2 * comp * dofs_per_face], + temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face], global_vector_ptr + ind1 + comp * static_dofs_per_component * VectorizedArrayType::size() + @@ -2529,11 +2547,13 @@ namespace internal &data.face_to_cell_index_nodal(face_no, 0); for (unsigned int i = 0; i < dofs_per_face; ++i) { + const unsigned int i_ = reorientate(i); + const unsigned int ind = index_array[i] * VectorizedArrayType::size(); for (unsigned int comp = 0; comp < n_components; ++comp) function_2b( - temp1[i + 2 * comp * dofs_per_face], + temp1[i_ + 2 * comp * dofs_per_face], global_vector_ptr + ind + comp * static_dofs_per_component * VectorizedArrayType::size() + @@ -2585,7 +2605,8 @@ namespace internal for (unsigned int comp = 0; comp < n_components; ++comp) for (unsigned int i = 0; i < dofs_per_face; ++i) { - unsigned int ind1[VectorizedArrayType::size()]; + const unsigned int i_ = reorientate(i); + unsigned int ind1[VectorizedArrayType::size()]; DEAL_II_OPENMP_SIMD_PRAGMA for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) @@ -2602,8 +2623,8 @@ namespace internal index_array[2 * i + 1]) * strides[v]; function_2a( - temp1[i + 2 * comp * dofs_per_face], - temp1[i + dofs_per_face + 2 * comp * dofs_per_face], + temp1[i_ + 2 * comp * dofs_per_face], + temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face], global_vector_ptr, global_vector_ptr, grad_weight, @@ -2622,6 +2643,7 @@ namespace internal for (unsigned int comp = 0; comp < n_components; ++comp) for (unsigned int i = 0; i < dofs_per_face; ++i) { + const unsigned int i_ = reorientate(i); const unsigned int ind1 = indices[v] + (comp * static_dofs_per_component + index_array[2 * i]) * @@ -2630,8 +2652,8 @@ namespace internal indices[v] + (comp * static_dofs_per_component + index_array[2 * i + 1]) * strides[v]; - function_3a(temp1[i + 2 * comp * dofs_per_face][v], - temp1[i + dofs_per_face + + function_3a(temp1[i_ + 2 * comp * dofs_per_face][v], + temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face][v], global_vector_ptr[ind1], global_vector_ptr[ind2], @@ -2657,7 +2679,8 @@ namespace internal indices[v] + (comp * static_dofs_per_component + index_array[i]) * strides[v]; - function_2b(temp1[i + 2 * comp * dofs_per_face], + const unsigned int i_ = reorientate(i); + function_2b(temp1[i_ + 2 * comp * dofs_per_face], global_vector_ptr, ind); } @@ -2676,7 +2699,8 @@ namespace internal indices[v] + (comp * static_dofs_per_component + index_array[i]) * strides[v]; - function_3b(temp1[i + 2 * comp * dofs_per_face][v], + const unsigned int i_ = reorientate(i); + function_3b(temp1[i_ + 2 * comp * dofs_per_face][v], global_vector_ptr[ind1]); } } @@ -2718,10 +2742,12 @@ namespace internal { const unsigned int ind1 = index_array[2 * i]; const unsigned int ind2 = index_array[2 * i + 1]; + const unsigned int i_ = reorientate(i); + for (unsigned int comp = 0; comp < n_components; ++comp) function_2a( - temp1[i + 2 * comp * dofs_per_face], - temp1[i + dofs_per_face + 2 * comp * dofs_per_face], + temp1[i_ + 2 * comp * dofs_per_face], + temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face], global_vector_ptr + comp * static_dofs_per_component + ind1 + dof_info.component_dof_indices_offset @@ -2745,7 +2771,9 @@ namespace internal for (unsigned int comp = 0; comp < n_components; ++comp) { const unsigned int ind = index_array[i]; - function_2b(temp1[i + 2 * comp * dofs_per_face], + const unsigned int i_ = reorientate(i); + + function_2b(temp1[i_ + 2 * comp * dofs_per_face], global_vector_ptr + comp * static_dofs_per_component + ind + dof_info.component_dof_indices_offset @@ -2765,7 +2793,9 @@ namespace internal if (!integrate) function_0(temp1, dofs_per_face); - if (!integrate && (face_orientation > 0)) + if (!integrate && + (face_orientation > 0 && + subface_index < GeometryInfo::max_children_per_cell)) adjust_for_face_orientation(face_orientation, orientation_map, false, diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index ccfe490211..f20abbc98d 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -467,6 +467,14 @@ namespace internal */ dealii::Table<2, unsigned int> face_to_cell_index_hermite; + /** + * For degrees on faces, the basis functions are not + * in the correct order if a face is not in the standard orientation + * to a given element. This data structure is used to re-order the + * basis functions to represent the correct order. + */ + dealii::Table<2, unsigned int> face_orientations; + private: /** * Check whether we have symmetries in the shape values. In that case, diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 4bced72888..8b12453bab 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -481,6 +481,46 @@ namespace internal face_to_cell_index_nodal(f, l) = ind; } } + + // face orientation for faces in 3D + // (similar to MappingInfoStorage::QuadratureDescriptor::initialize) + if (dim == 3) + { + const unsigned int n = fe_degree + 1; + face_orientations.reinit(8, n * n); + for (unsigned int j = 0, i = 0; j < n; ++j) + for (unsigned int k = 0; k < n; ++k, ++i) + { + // face_orientation=true, face_flip=false, + // face_rotation=false + face_orientations[0][i] = i; + // face_orientation=false, face_flip=false, + // face_rotation=false + face_orientations[1][i] = j + k * n; + // face_orientation=true, face_flip=true, + // face_rotation=false + face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n; + // face_orientation=false, face_flip=true, + // face_rotation=false + face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n; + // face_orientation=true, face_flip=false, + // face_rotation=true + face_orientations[4][i] = j + (n - 1 - k) * n; + // face_orientation=false, face_flip=false, + // face_rotation=true + face_orientations[5][i] = k + (n - 1 - j) * n; + // face_orientation=true, face_flip=true, + // face_rotation=true + face_orientations[6][i] = (n - 1 - j) + k * n; + // face_orientation=false, face_flip=true, + // face_rotation=true + face_orientations[7][i] = (n - 1 - k) + j * n; + } + } + else + { + face_orientations.reinit(1, 1); + } } if (element_type == tensor_symmetric_hermite)