From: Martin Kronbichler Date: Wed, 2 Aug 2017 10:19:50 +0000 (+0200) Subject: Generate vectorized data on 1D faces for MF::shape_info. X-Git-Tag: v9.0.0-rc1~1364^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=db65fbc39276c0a087a2795df7e6dd40efcdf603;p=dealii.git Generate vectorized data on 1D faces for MF::shape_info. --- diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index e386e1a5ad..c942b1fe6c 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -780,7 +780,7 @@ namespace MatrixFreeOperators FullMatrix shapes_1d(fe_degree+1, fe_degree+1); for (unsigned int i=0, c=0; i #include -#include - DEAL_II_NAMESPACE_OPEN @@ -190,40 +188,29 @@ namespace internal AlignedVector > shape_hessians_collocation_eo; /** - * Stores the indices from cell DoFs to face DoFs. The rows go through - * the 2*dim faces, and the columns the DoFs on the faces. - */ - dealii::Table<2,unsigned int> face_indices; - - /** - * Stores one-dimensional values of shape functions evaluated in zero - * and one, i.e., on the one-dimensional faces. Not vectorized. + * Collects all data of 1D shape values evaluated at the point 0 and 1 + * (the vertices) in one data structure. Sorting is first the values, + * then gradients, then second derivatives. */ - std::vector face_value[2]; - - /** - * Stores one-dimensional gradients of shape functions evaluated in zero - * and one, i.e., on the one-dimensional faces. Not vectorized. - */ - std::vector face_gradient[2]; + AlignedVector > shape_data_on_face[2]; /** * Stores one-dimensional values of shape functions on subface. Since - * there are two subfaces, store two variants. Not vectorized. + * there are two subfaces, store two variants. */ - std::vector subface_value[2]; + AlignedVector > values_within_subface[2]; /** - * Non-vectorized version of shape values. Needed when evaluating face - * info. + * Stores one-dimensional gradients of shape functions on subface. Since + * there are two subfaces, store two variants. */ - std::vector shape_values_number; + AlignedVector > gradients_within_subface[2]; /** - * Non-vectorized version of shape gradients. Needed when evaluating - * face info. + * Stores one-dimensional gradients of shape functions on subface. Since + * there are two subfaces, store two variants. */ - std::vector shape_gradient_number; + AlignedVector > hessians_within_subface[2]; /** * Renumbering from deal.II's numbering of cell degrees of freedom to @@ -265,6 +252,43 @@ namespace internal */ unsigned int dofs_per_face; + /** + * Indicates whether the basis functions are nodal in 0 and 1, i.e., the + * end points of the unit cell. + */ + bool nodal_at_cell_boundaries; + + /** + * For nodal cells, we might get around by simply loading the indices to + * the degrees of freedom that act on a particular face, rather than the + * whole set of indices that is then interpolated down to the + * element. This array stores this indirect addressing. + * + * The first table index runs through the faces of a cell, and the + * second runs through the nodal degrees of freedom of the face, using + * @p dofs_per_face entries. + * + * @note This object is only filled in case @p nodal_at_cell_boundaries + * evaluates to @p true. + */ + dealii::Table<2,unsigned int> face_to_cell_index_nodal; + + /** + * For Hermite-type basis functions, the @p face_to_cell_index_nodal for + * the values on a face of the cell is used together with a respective + * slot in the derivative. In the lexicographic ordering, this index is + * in the next "layer" of degrees of freedom. This array stores the + * indirect addressing of both the values and the gradient. + * + * The first table index runs through the faces of a cell, and the + * second runs through the pairs of the nodal degrees of freedom of the + * face and the derivatives, using 2*dofs_per_face entries. + * + * @note This object is only filled in case @p element_type evaluates to + * @p tensor_symmetric_hermite. + */ + dealii::Table<2,unsigned int> face_to_cell_index_hermite; + /** * Check whether we have symmetries in the shape values. In that case, * also fill the shape_???_eo fields. @@ -297,7 +321,8 @@ namespace internal n_q_points (0), dofs_per_cell (0), n_q_points_face (0), - dofs_per_face (0) + dofs_per_face (0), + nodal_at_cell_boundaries (false) { reinit (quad, fe_in, base_element_number); } diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index d29e75dd50..cefbd43e06 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -49,7 +49,8 @@ namespace internal n_q_points (0), dofs_per_cell (0), n_q_points_face (0), - dofs_per_face (0) + dofs_per_face (0), + nodal_at_cell_boundaries (false) {} @@ -69,8 +70,7 @@ namespace internal fe_degree = fe->degree; n_q_points_1d = quad.size(); - const unsigned int n_dofs_1d = fe_degree+1, - n_q_points_1d = quad.size(); + const unsigned int n_dofs_1d = std::min(fe->dofs_per_cell, fe_degree+1); // renumber (this is necessary for FE_Q, for example, since there the // vertex DoFs come first, which is incompatible with the lexicographic @@ -111,6 +111,10 @@ namespace internal scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse(); element_type = tensor_symmetric_plus_dg0; } + else if (fe->dofs_per_cell == 0) + { + // FE_Nothing case -> nothing to do here + } else Assert(false, ExcNotImplemented()); @@ -138,7 +142,7 @@ namespace internal // invert numbering again. Need to do it manually because we might // have undefined blocks - lexicographic_numbering.resize(fe_in.element_multiplicity(base_element_number)*fe->dofs_per_cell); + lexicographic_numbering.resize(fe_in.element_multiplicity(base_element_number)*fe->dofs_per_cell, numbers::invalid_unsigned_int); for (unsigned int i=0; ihas_support_points()) unit_point = fe->get_unit_support_points()[scalar_lexicographic[0]]; - Assert(std::fabs(fe->shape_value(scalar_lexicographic[0], + Assert(fe->dofs_per_cell == 0 || + std::fabs(fe->shape_value(scalar_lexicographic[0], unit_point)-1) < 1e-13, - ExcInternalError()); + ExcInternalError("Could not decode 1D shape functions for the " + "element " + fe->get_name())); } n_q_points = Utilities::fixed_power(n_q_points_1d); dofs_per_cell = fe->dofs_per_cell; n_q_points_face = dim>1?Utilities::fixed_power(n_q_points_1d):1; - dofs_per_face = fe->dofs_per_face; + dofs_per_face = dim>1?Utilities::fixed_power(fe_degree+1):1; const unsigned int array_size = n_dofs_1d*n_q_points_1d; this->shape_gradients.resize_fast (array_size); this->shape_values.resize_fast (array_size); this->shape_hessians.resize_fast (array_size); - this->face_value[0].resize(n_dofs_1d); - this->face_gradient[0].resize(n_dofs_1d); - this->subface_value[0].resize(array_size); - this->face_value[1].resize(n_dofs_1d); - this->face_gradient[1].resize(n_dofs_1d); - this->subface_value[1].resize(array_size); - this->shape_values_number.resize (array_size); - this->shape_gradient_number.resize (array_size); + this->shape_data_on_face[0].resize(3*n_dofs_1d); + this->shape_data_on_face[1].resize(3*n_dofs_1d); + this->values_within_subface[0].resize(array_size); + this->values_within_subface[1].resize(array_size); + this->gradients_within_subface[0].resize(array_size); + this->gradients_within_subface[1].resize(array_size); + this->hessians_within_subface[0].resize(array_size); + this->hessians_within_subface[1].resize(array_size); for (unsigned int i=0; i::n_array_elements - // copies for the shape information and - // non-vectorized fields + // VectorizedArray::n_array_elements copies for the + // shape information and non-vectorized fields Point q_point = unit_point; q_point[0] = quad.get_points()[q][0]; - shape_values_number[i*n_q_points_1d+q] = fe->shape_value(my_i,q_point); - shape_gradient_number[i*n_q_points_1d+q] = fe->shape_grad (my_i,q_point)[0]; - shape_values [i*n_q_points_1d+q] = - shape_values_number [i*n_q_points_1d+q]; - shape_gradients[i*n_q_points_1d+q] = - shape_gradient_number[i*n_q_points_1d+q]; - shape_hessians[i*n_q_points_1d+q] = - fe->shape_grad_grad(my_i,q_point)[0][0]; + + shape_values [i*n_q_points_1d+q] = fe->shape_value(my_i,q_point); + shape_gradients[i*n_q_points_1d+q] = fe->shape_grad(my_i,q_point)[0]; + shape_hessians [i*n_q_points_1d+q] = fe->shape_grad_grad(my_i,q_point)[0][0]; + + // evaluate basis functions on the two 1D subfaces (i.e., at the + // positions divided by one half and shifted by one half, + // respectively) q_point[0] *= 0.5; - subface_value[0][i*n_q_points_1d+q] = fe->shape_value(my_i,q_point); + values_within_subface[0][i*n_q_points_1d+q] = fe->shape_value(my_i,q_point); + gradients_within_subface[0][i*n_q_points_1d+q] = fe->shape_grad(my_i,q_point)[0]; + hessians_within_subface[0][i*n_q_points_1d+q] = fe->shape_grad_grad(my_i,q_point)[0][0]; q_point[0] += 0.5; - subface_value[1][i*n_q_points_1d+q] = fe->shape_value(my_i,q_point); + values_within_subface[1][i*n_q_points_1d+q] = fe->shape_value(my_i,q_point); + gradients_within_subface[1][i*n_q_points_1d+q] = fe->shape_grad(my_i,q_point)[0]; + hessians_within_subface[1][i*n_q_points_1d+q] = fe->shape_grad_grad(my_i,q_point)[0][0]; } - Point q_point; - this->face_value[0][i] = fe->shape_value(my_i,q_point); - this->face_gradient[0][i] = fe->shape_grad(my_i,q_point)[0]; + + // evaluate basis functions on the 1D faces, i.e., in zero and one + Point q_point = unit_point; + q_point[0] = 0; + this->shape_data_on_face[0][i] = fe->shape_value(my_i,q_point); + this->shape_data_on_face[0][i+n_dofs_1d] = fe->shape_grad(my_i,q_point)[0]; + this->shape_data_on_face[0][i+2*n_dofs_1d] = fe->shape_grad_grad(my_i,q_point)[0][0]; q_point[0] = 1; - this->face_value[1][i] = fe->shape_value(my_i,q_point); - this->face_gradient[1][i] = fe->shape_grad(my_i,q_point)[0]; + this->shape_data_on_face[1][i] = fe->shape_value(my_i,q_point); + this->shape_data_on_face[1][i+n_dofs_1d] = fe->shape_grad(my_i,q_point)[0]; + this->shape_data_on_face[1][i+2*n_dofs_1d] = fe->shape_grad_grad(my_i,q_point)[0][0]; } + // get gradient and Hessian transformation matrix for the polynomial + // space associated with the quadrature rule (collocation space) + { + const unsigned int stride = (n_q_points_1d+1)/2; + shape_gradients_collocation_eo.resize(n_q_points_1d*stride); + shape_hessians_collocation_eo.resize(n_q_points_1d*stride); + FE_DGQArbitraryNodes<1> fe(quad.get_points()); + for (unsigned int i=0; i 3 && element_type == tensor_symmetric) { - element_type = tensor_symmetric; - // get gradient and Hessian transformation matrix for the - // polynomial space associated with the quadrature rule - // (collocation space) - if (fe_degree+1 == n_q_points_1d) - { - const unsigned int stride = fe_degree/2+1; - shape_gradients_collocation_eo.resize((fe_degree+1)*stride); - shape_hessians_collocation_eo.resize((fe_degree+1)*stride); - FE_DGQArbitraryNodes<1> fe(quad.get_points()); - for (unsigned int i=0; i<(fe_degree+1)/2; ++i) - for (unsigned int q=0; qshape_data_on_face[0][i][0]) > 1e-12) + element_type = tensor_symmetric; + for (unsigned int i=2; ishape_data_on_face[0][n_dofs_1d+i][0]) > 1e-12) + element_type = tensor_symmetric; } } else if (element_type == tensor_symmetric_plus_dg0) check_1d_shapes_symmetric(n_q_points_1d); - // face information - unsigned int n_faces = GeometryInfo::faces_per_cell; - this->face_indices.reinit(n_faces, this->dofs_per_face); - switch (dim) - { - case 3: - { - for (unsigned int i=0; idofs_per_face; i++) - { - const unsigned int jump_term = - this->dofs_per_face*((i*n_dofs_1d)/this->dofs_per_face); - this->face_indices(0,i) = i*n_dofs_1d; - this->face_indices(1,i) = i*n_dofs_1d + n_dofs_1d-1; - this->face_indices(2,i) = i%n_dofs_1d + jump_term; - this->face_indices(3,i) = (i%n_dofs_1d + jump_term + - (n_dofs_1d-1)*n_dofs_1d); - this->face_indices(4,i) = i; - this->face_indices(5,i) = (n_dofs_1d-1)*this->dofs_per_face+i; - } - break; - } - case 2: + nodal_at_cell_boundaries = true; + for (unsigned int i=1; ishape_data_on_face[0][i][0]) > 1e-13 || + std::abs(this->shape_data_on_face[1][i-1][0]) > 1e-13) + nodal_at_cell_boundaries = false; + + if (nodal_at_cell_boundaries == true) { - for (unsigned int i=0; idofs_per_face; i++) + face_to_cell_index_nodal.reinit(GeometryInfo::faces_per_cell, + dofs_per_face); + for (unsigned int f=0; f::faces_per_cell; ++f) { - this->face_indices(0,i) = n_dofs_1d*i; - this->face_indices(1,i) = n_dofs_1d*i + n_dofs_1d-1; - this->face_indices(2,i) = i; - this->face_indices(3,i) = (n_dofs_1d-1)*n_dofs_1d+i; + const unsigned int direction = f/2; + const unsigned int stride = direction < dim-1 ? (fe_degree+1) : 1; + int shift = 1; + for (unsigned int d=0; d adjust + // that here + for (unsigned int j=0; j<=fe_degree; ++j) + for (unsigned int i=0; i<=fe_degree; ++i) + { + const unsigned int ind = offset + j*dofs_per_face + i; + AssertIndexRange(ind, dofs_per_cell); + const unsigned int l = i*(fe_degree+1)+j; + face_to_cell_index_nodal(f,l) = ind; + } } - break; } - case 1: + + if (element_type == tensor_symmetric_hermite) { - if (this->dofs_per_face>0) + face_to_cell_index_hermite.reinit(GeometryInfo::faces_per_cell, + 2*dofs_per_face); + for (unsigned int f=0; f::faces_per_cell; ++f) { - this->face_indices(0,0) = 0; - this->face_indices(1,0) = n_dofs_1d-1; + const unsigned int direction = f/2; + const unsigned int stride = direction < dim-1 ? (fe_degree+1) : 1; + int shift = 1; + for (unsigned int d=0; d adjust + // that here + for (unsigned int j=0; j<=fe_degree; ++j) + for (unsigned int i=0; i<=fe_degree; ++i) + { + const unsigned int ind = offset + j*dofs_per_face + i; + AssertIndexRange(ind, dofs_per_cell); + const unsigned int l = i*(fe_degree+1)+j; + face_to_cell_index_hermite(f,2*l) = ind; + face_to_cell_index_hermite(f,2*l+1) = ind+shift; + } } - break; - } - default: - Assert (false, ExcNotImplemented()); } } @@ -311,15 +362,20 @@ namespace internal bool ShapeInfo::check_1d_shapes_symmetric(const unsigned int n_q_points_1d) { + if (dofs_per_cell == 0) + return false; + const double zero_tol = - types_are_equal::value==true?1e-10:1e-7; + types_are_equal::value==true?1e-12:1e-7; // symmetry for values const unsigned int n_dofs_1d = fe_degree + 1; for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) for (unsigned int j=0; j zero_tol) + -j-1][0]) > + std::max(zero_tol, zero_tol* + std::abs(shape_values[i*n_q_points_1d+j][0]))) return false; // shape values should be zero at x=0.5 for all basis functions except @@ -412,10 +468,8 @@ namespace internal return false; const double zero_tol = - types_are_equal::value==true?1e-10:1e-7; - // check: - identity operation for shape values - // - zero diagonal at interior points for gradients - // - gradient equal to unity at element boundary + types_are_equal::value==true?1e-12:1e-7; + // check: identity operation for shape values const unsigned int n_points_1d = fe_degree+1; for (unsigned int i=0; i