From: Peter Munch Date: Wed, 25 Sep 2019 10:55:34 +0000 (+0200) Subject: Compute shape_gradients_collocation and shape_hessians_collocation in ShapeInfo X-Git-Tag: v9.2.0-rc1~1017^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=dd5f8d992eede912002fcafeb7af37e7019dae5d;p=dealii.git Compute shape_gradients_collocation and shape_hessians_collocation in ShapeInfo --- diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index a601d2a591..694008faa0 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -172,6 +172,18 @@ namespace internal */ AlignedVector shape_hessians; + /** + * Stores the shape gradients of the shape function space associated to + * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). + */ + AlignedVector shape_gradients_collocation; + + /** + * Stores the shape hessians of the shape function space associated to + * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). + */ + AlignedVector shape_hessians_collocation; + /** * Stores the shape values in a different format, namely the so-called * even-odd scheme where the symmetries in shape_values are used for @@ -195,15 +207,17 @@ namespace internal /** * Stores the shape gradients of the shape function space associated to - * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For - * faster evaluation only the even-odd format is necessary. + * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). This + * array provides an alternative representation of the + * shape_gradients_collocation field in the even-odd format. */ AlignedVector shape_gradients_collocation_eo; /** * Stores the shape hessians of the shape function space associated to - * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For - * faster evaluation only the even-odd format is necessary. + * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). This + * array provides an alternative representation of the + * shape_hessians_collocation field in the even-odd format. */ AlignedVector shape_hessians_collocation_eo; diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 8459115ac5..a3e38677ae 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -272,43 +272,16 @@ namespace internal // FE_DGQArbitraryNodes underflow. if (n_q_points_1d < 200) { - 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); + shape_gradients_collocation.resize(n_q_points_1d * n_q_points_1d); + shape_hessians_collocation.resize(n_q_points_1d * n_q_points_1d); FE_DGQArbitraryNodes<1> fe(quad.get_points()); - for (unsigned int i = 0; i < n_q_points_1d / 2; ++i) - for (unsigned int q = 0; q < stride; ++q) + for (unsigned int i = 0; i < n_q_points_1d; ++i) + for (unsigned int q = 0; q < n_q_points_1d; ++q) { - shape_gradients_collocation_eo[i * stride + q] = - 0.5 * - (fe.shape_grad(i, quad.get_points()[q])[0] + - fe.shape_grad(i, - quad.get_points()[n_q_points_1d - 1 - q])[0]); - shape_gradients_collocation_eo[(n_q_points_1d - 1 - i) * - stride + - q] = - 0.5 * - (fe.shape_grad(i, quad.get_points()[q])[0] - - fe.shape_grad(i, - quad.get_points()[n_q_points_1d - 1 - q])[0]); - shape_hessians_collocation_eo[i * stride + q] = - 0.5 * (fe.shape_grad_grad(i, quad.get_points()[q])[0][0] + - fe.shape_grad_grad( - i, quad.get_points()[n_q_points_1d - 1 - q])[0][0]); - shape_hessians_collocation_eo[(n_q_points_1d - 1 - i) * stride + - q] = - 0.5 * (fe.shape_grad_grad(i, quad.get_points()[q])[0][0] - - fe.shape_grad_grad( - i, quad.get_points()[n_q_points_1d - 1 - q])[0][0]); - } - if (n_q_points_1d % 2 == 1) - for (unsigned int q = 0; q < stride; ++q) - { - shape_gradients_collocation_eo[n_q_points_1d / 2 * stride + q] = - fe.shape_grad(n_q_points_1d / 2, quad.get_points()[q])[0]; - shape_hessians_collocation_eo[n_q_points_1d / 2 * stride + q] = - fe.shape_grad_grad(n_q_points_1d / 2, - quad.get_points()[q])[0][0]; + shape_gradients_collocation[i * n_q_points_1d + q] = + fe.shape_grad(i, quad.get_points()[q])[0]; + shape_hessians_collocation[i * n_q_points_1d + q] = + fe.shape_grad_grad(i, quad.get_points()[q])[0][0]; } } @@ -483,46 +456,51 @@ namespace internal zero_tol_hessian) return false; - const unsigned int stride = (n_q_points_1d + 1) / 2; - shape_values_eo.resize((fe_degree + 1) * stride); - shape_gradients_eo.resize((fe_degree + 1) * stride); - shape_hessians_eo.resize((fe_degree + 1) * stride); - for (unsigned int i = 0; i < (fe_degree + 1) / 2; ++i) - for (unsigned int q = 0; q < stride; ++q) - { - shape_values_eo[i * stride + q] = - 0.5 * (shape_values[i * n_q_points_1d + q] + - shape_values[i * n_q_points_1d + n_q_points_1d - 1 - q]); - shape_values_eo[(fe_degree - i) * stride + q] = - 0.5 * (shape_values[i * n_q_points_1d + q] - - shape_values[i * n_q_points_1d + n_q_points_1d - 1 - q]); - - shape_gradients_eo[i * stride + q] = - 0.5 * - (shape_gradients[i * n_q_points_1d + q] + - shape_gradients[i * n_q_points_1d + n_q_points_1d - 1 - q]); - shape_gradients_eo[(fe_degree - i) * stride + q] = - 0.5 * - (shape_gradients[i * n_q_points_1d + q] - - shape_gradients[i * n_q_points_1d + n_q_points_1d - 1 - q]); - - shape_hessians_eo[i * stride + q] = - 0.5 * (shape_hessians[i * n_q_points_1d + q] + - shape_hessians[i * n_q_points_1d + n_q_points_1d - 1 - q]); - shape_hessians_eo[(fe_degree - i) * stride + q] = - 0.5 * (shape_hessians[i * n_q_points_1d + q] - - shape_hessians[i * n_q_points_1d + n_q_points_1d - 1 - q]); - } - if (fe_degree % 2 == 0) - for (unsigned int q = 0; q < stride; ++q) - { - shape_values_eo[fe_degree / 2 * stride + q] = - shape_values[(fe_degree / 2) * n_q_points_1d + q]; - shape_gradients_eo[fe_degree / 2 * stride + q] = - shape_gradients[(fe_degree / 2) * n_q_points_1d + q]; - shape_hessians_eo[fe_degree / 2 * stride + q] = - shape_hessians[(fe_degree / 2) * n_q_points_1d + q]; - } + auto convert_to_eo = [](const AlignedVector &array, + const unsigned n_points_dst, + const unsigned n_points_src) { + const unsigned int stride = (n_points_dst + 1) / 2; + AlignedVector array_eo(n_points_src * stride); + for (unsigned int i = 0; i < n_points_src / 2; ++i) + for (unsigned int q = 0; q < stride; ++q) + { + array_eo[i * stride + q] = + 0.5 * (array[i * n_points_dst + q] + + array[i * n_points_dst + n_points_dst - 1 - q]); + array_eo[(n_points_src - 1 - i) * stride + q] = + 0.5 * (array[i * n_points_dst + q] - + array[i * n_points_dst + n_points_dst - 1 - q]); + } + if ((n_points_src - 1) % 2 == 0) + for (unsigned int q = 0; q < stride; ++q) + { + array_eo[(n_points_src - 1) / 2 * stride + q] = + array[((n_points_src - 1) / 2) * n_points_dst + q]; + } + + return array_eo; + }; + + shape_values_eo = + convert_to_eo(shape_values, n_q_points_1d, fe_degree + 1); + shape_gradients_eo = + convert_to_eo(shape_gradients, n_q_points_1d, fe_degree + 1); + shape_hessians_eo = + convert_to_eo(shape_hessians, n_q_points_1d, fe_degree + 1); + + // FE_DGQArbitraryNodes underflow (see also above where + // shape_gradients_collocation and shape_hessians_collocation is set up). + if (n_q_points_1d < 200) + { + shape_gradients_collocation_eo = + convert_to_eo(shape_gradients_collocation, + n_q_points_1d, + n_q_points_1d); + shape_hessians_collocation_eo = + convert_to_eo(shape_hessians_collocation, + n_q_points_1d, + n_q_points_1d); + } return true; } @@ -568,6 +546,10 @@ namespace internal memory += MemoryConsumption::memory_consumption(shape_values); memory += MemoryConsumption::memory_consumption(shape_gradients); memory += MemoryConsumption::memory_consumption(shape_hessians); + memory += + MemoryConsumption::memory_consumption(shape_gradients_collocation); + memory += + MemoryConsumption::memory_consumption(shape_hessians_collocation); memory += MemoryConsumption::memory_consumption(shape_values_eo); memory += MemoryConsumption::memory_consumption(shape_gradients_eo); memory += MemoryConsumption::memory_consumption(shape_hessians_eo);