1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2011 - 2020 by the deal.II authors 4 // 5 // This file is part of the deal.II library. 6 // 7 // The deal.II library is free software; you can use it, redistribute 8 // it, and/or modify it under the terms of the GNU Lesser General 9 // Public License as published by the Free Software Foundation; either 10 // version 2.1 of the License, or (at your option) any later version. 11 // The full text of the license can be found in the file LICENSE.md at 12 // the top level directory of deal.II. 13 // 14 // --------------------------------------------------------------------- 15 16 17 #ifndef dealii_matrix_free_shape_info_h 18 #define dealii_matrix_free_shape_info_h 19 20 21 #include <deal.II/base/config.h> 22 23 #include <deal.II/base/aligned_vector.h> 24 #include <deal.II/base/exceptions.h> 25 #include <deal.II/base/quadrature_lib.h> 26 27 #include <deal.II/fe/fe.h> 28 29 30 DEAL_II_NAMESPACE_OPEN 31 32 namespace internal 33 { 34 namespace MatrixFreeFunctions 35 { 36 /** 37 * An enum that encodes the type of element detected during 38 * initialization. FEEvaluation will select the most efficient algorithm 39 * based on the given element type. 40 * 41 * There is an implied ordering in the type ElementType::tensor_symmetric 42 * in the sense that both ElementType::tensor_symmetric_collocation and 43 * ElementType::tensor_symmetric_hermite are also of type 44 * ElementType::tensor_symmetric. Likewise, a configuration of type 45 * ElementType::tensor_symmetric is also of type 46 * ElementType::tensor_general. As a consequence, we support `<=` 47 * operations between the types with this sorting, but not against the 48 * even higher indexed types such as ElementType::truncated_tensor. 49 * 50 * @ingroup matrixfree 51 */ 52 enum ElementType 53 { 54 /** 55 * Tensor product shape function where the shape value evaluation in the 56 * quadrature point corresponds to the identity operation and no 57 * interpolation needs to be performed (collocation approach, also 58 * called spectral evaluation). This is for example the case for an 59 * element with nodes in the Gauss-Lobatto support points and 60 * integration in the Gauss-Lobatto quadrature points of the same order. 61 */ 62 tensor_symmetric_collocation = 0, 63 64 /** 65 * Symmetric tensor product shape functions fulfilling a Hermite 66 * identity with values and first derivatives zero at the element end 67 * points in 1D. 68 */ 69 tensor_symmetric_hermite = 1, 70 71 /** 72 * Usual tensor product shape functions whose shape values and 73 * quadrature points are symmetric about the midpoint of the unit 74 * interval 0.5 75 */ 76 tensor_symmetric = 2, 77 78 /** 79 * Tensor product shape functions without further particular properties 80 */ 81 tensor_general = 3, 82 83 /** 84 * Polynomials of complete degree rather than tensor degree which can be 85 * described by a truncated tensor product 86 */ 87 truncated_tensor = 4, 88 89 /** 90 * Tensor product shape functions that are symmetric about the midpoint 91 * of the unit interval 0.5 that additionally add a constant shape 92 * function according to FE_Q_DG0. 93 */ 94 tensor_symmetric_plus_dg0 = 5 95 }; 96 97 98 99 /** 100 * This struct stores the shape functions, their gradients and Hessians 101 * evaluated for a one-dimensional section of a tensor product finite 102 * element and tensor product quadrature formula in reference 103 * coordinates. This data structure also includes the evaluation of 104 * quantities at the cell boundary and on the sub-interval $(0, 0.5)$ and 105 * $(0.5, 1)$ for face integrals. 106 */ 107 template <typename Number> 108 struct UnivariateShapeData 109 { 110 /** 111 * Empty constructor. Sets default configuration. 112 */ 113 UnivariateShapeData(); 114 115 /** 116 * Return the memory consumption of this class in bytes. 117 */ 118 std::size_t 119 memory_consumption() const; 120 121 /** 122 * Encodes the type of element detected at construction. FEEvaluation 123 * will select the most efficient algorithm based on the given element 124 * type. 125 */ 126 ElementType element_type; 127 128 /** 129 * Stores the shape values of the 1D finite element evaluated on all 1D 130 * quadrature points in vectorized format, i.e., as an array of 131 * VectorizedArray<dim>::size equal elements. The length of 132 * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature 133 * points are the index running fastest. 134 */ 135 AlignedVector<Number> shape_values; 136 137 /** 138 * Stores the shape gradients of the 1D finite element evaluated on all 139 * 1D quadrature points in vectorized format, i.e., as an array of 140 * VectorizedArray<dim>::size equal elements. The length of 141 * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature 142 * points are the index running fastest. 143 */ 144 AlignedVector<Number> shape_gradients; 145 146 /** 147 * Stores the shape Hessians of the 1D finite element evaluated on all 148 * 1D quadrature points in vectorized format, i.e., as an array of 149 * VectorizedArray<dim>::size equal elements. The length of 150 * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature 151 * points are the index running fastest. 152 */ 153 AlignedVector<Number> shape_hessians; 154 155 /** 156 * Stores the shape gradients of the shape function space associated to 157 * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). 158 */ 159 AlignedVector<Number> shape_gradients_collocation; 160 161 /** 162 * Stores the shape hessians of the shape function space associated to 163 * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). 164 */ 165 AlignedVector<Number> shape_hessians_collocation; 166 167 /** 168 * Stores the shape values in a different format, namely the so-called 169 * even-odd scheme where the symmetries in shape_values are used for 170 * faster evaluation. 171 */ 172 AlignedVector<Number> shape_values_eo; 173 174 /** 175 * Stores the shape gradients in a different format, namely the so- 176 * called even-odd scheme where the symmetries in shape_gradients are 177 * used for faster evaluation. 178 */ 179 AlignedVector<Number> shape_gradients_eo; 180 181 /** 182 * Stores the shape second derivatives in a different format, namely the 183 * so-called even-odd scheme where the symmetries in shape_hessians are 184 * used for faster evaluation. 185 */ 186 AlignedVector<Number> shape_hessians_eo; 187 188 /** 189 * Stores the shape gradients of the shape function space associated to 190 * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). This 191 * array provides an alternative representation of the 192 * shape_gradients_collocation field in the even-odd format. 193 */ 194 AlignedVector<Number> shape_gradients_collocation_eo; 195 196 /** 197 * Stores the shape hessians of the shape function space associated to 198 * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). This 199 * array provides an alternative representation of the 200 * shape_hessians_collocation field in the even-odd format. 201 */ 202 AlignedVector<Number> shape_hessians_collocation_eo; 203 204 /** 205 * Stores the inverse transformation from the data at quadrature points 206 * to the basis defined by the shape_values fields. The data at 207 * quadrature points is interpreted either implicitly by its polynomial 208 * interpolation, or explicitly in terms of separate polynomials such as 209 * with the `_collocation` fields. The size of the array equals the 210 * layout of the `shape_values` array, and it is combined with the shape 211 * values array such that this matrix is the pseudo inverse of 212 * shape_values. In case the number of 1D quadrature points equals the 213 * size of the basis, this array is exactly the inverse of the 214 * shape_values array. The length of this array is <tt>n_dofs_1d * 215 * n_q_points_1d</tt> and quadrature points are the index running 216 * fastest. 217 */ 218 AlignedVector<Number> inverse_shape_values; 219 220 /** 221 * Stores the even-odd variant of the `inverse_shape_values` field. 222 */ 223 AlignedVector<Number> inverse_shape_values_eo; 224 225 /** 226 * Collects all data of 1D shape values evaluated at the point 0 and 1 227 * (the vertices) in one data structure. Sorting is first the values, 228 * then gradients, then second derivatives. 229 */ 230 std::array<AlignedVector<Number>, 2> shape_data_on_face; 231 232 /** 233 * Collects all data of 1D nodal shape values (defined by the Lagrange 234 * polynomials in the points of the quadrature rule) evaluated at the 235 * point 0 and 1 (the vertices) in one data structure. 236 * 237 * This data structure can be used to interpolate from the cell to the 238 * face quadrature points. 239 * 240 * @note In contrast to shape_data_on_face, only the vales are evaluated. 241 */ 242 std::array<AlignedVector<Number>, 2> quadrature_data_on_face; 243 244 /** 245 * Stores one-dimensional values of shape functions on subface. Since 246 * there are two subfaces, store two variants. 247 */ 248 std::array<AlignedVector<Number>, 2> values_within_subface; 249 250 /** 251 * Stores one-dimensional gradients of shape functions on subface. Since 252 * there are two subfaces, store two variants. 253 */ 254 std::array<AlignedVector<Number>, 2> gradients_within_subface; 255 256 /** 257 * Stores one-dimensional gradients of shape functions on subface. Since 258 * there are two subfaces, store two variants. 259 */ 260 std::array<AlignedVector<Number>, 2> hessians_within_subface; 261 262 /** 263 * We store a copy of the one-dimensional quadrature formula 264 * used for initialization. 265 */ 266 Quadrature<1> quadrature; 267 268 /** 269 * Stores the degree of the element. 270 */ 271 unsigned int fe_degree; 272 273 /** 274 * Stores the number of quadrature points per dimension. 275 */ 276 unsigned int n_q_points_1d; 277 278 /** 279 * Indicates whether the basis functions are nodal in 0 and 1, i.e., the 280 * end points of the unit cell. 281 */ 282 bool nodal_at_cell_boundaries; 283 }; 284 285 286 287 /** 288 * This struct stores a tensor (Kronecker) product view of the finite 289 * element and quadrature formula used for evaluation. It is based on a 290 * single or a collection of UnivariateShapeData object(s) that describe 291 * one-dimensional ingredients, plus some additional information about how 292 * these are combined and how indices are laid out in the multi-dimensional 293 * case such as the hierarchical -> lexicographic ordering of FE_Q. 294 * 295 * @ingroup matrixfree 296 */ 297 template <typename Number> 298 struct ShapeInfo 299 { 300 /** 301 * Encodes the type of element detected at construction. FEEvaluation 302 * will select the most efficient algorithm based on the given element 303 * type. 304 */ 305 ElementType element_type; 306 307 /** 308 * Empty constructor. Does nothing. 309 */ 310 ShapeInfo(); 311 312 /** 313 * Constructor that initializes the data fields using the reinit method. 314 */ 315 template <int dim> 316 ShapeInfo(const Quadrature<1> & quad, 317 const FiniteElement<dim> &fe, 318 const unsigned int base_element = 0); 319 320 /** 321 * Initializes the data fields. Takes a one-dimensional quadrature 322 * formula and a finite element as arguments and evaluates the shape 323 * functions, gradients and Hessians on the one-dimensional unit cell. 324 * This function assumes that the finite element is derived from a one- 325 * dimensional element by a tensor product and that the zeroth shape 326 * function in zero evaluates to one. 327 */ 328 template <int dim> 329 void 330 reinit(const Quadrature<1> & quad, 331 const FiniteElement<dim> &fe_dim, 332 const unsigned int base_element = 0); 333 334 /** 335 * Return data of univariate shape functions which defines the 336 * dimension @p dimension of tensor product shape functions 337 * regarding vector component @p component of the underlying 338 * finite element. 339 */ 340 const UnivariateShapeData<Number> & 341 get_shape_data(const unsigned int dimension = 0, 342 const unsigned int component = 0) const; 343 344 /** 345 * Return the memory consumption of this class in bytes. 346 */ 347 std::size_t 348 memory_consumption() const; 349 350 /** 351 * Renumbering from deal.II's numbering of cell degrees of freedom to 352 * lexicographic numbering used inside the FEEvaluation schemes of the 353 * underlying element in the DoFHandler. For vector-valued elements, the 354 * renumbering starts with a lexicographic numbering of the first 355 * component, then everything of the second component, and so on. 356 */ 357 std::vector<unsigned int> lexicographic_numbering; 358 359 /** 360 * Stores data of univariate shape functions defining the 361 * underlying tensor product finite element. 362 */ 363 std::vector<UnivariateShapeData<Number>> data; 364 365 /** 366 * Grants access to univariate shape function data of given 367 * dimension and vector component. Rows identify dimensions and 368 * columns identify vector components. 369 */ 370 dealii::Table<2, UnivariateShapeData<Number> *> data_access; 371 372 /** 373 * Stores the number of space dimensions. 374 */ 375 unsigned int n_dimensions; 376 377 /** 378 * Stores the number of vector components of the underlying 379 * vector-valued finite element. 380 */ 381 unsigned int n_components; 382 383 /** 384 * Stores the number of quadrature points in @p dim dimensions for a 385 * cell. 386 */ 387 unsigned int n_q_points; 388 389 /** 390 * Stores the number of DoFs per cell of the scalar element in @p dim 391 * dimensions. 392 */ 393 unsigned int dofs_per_component_on_cell; 394 395 /** 396 * Stores the number of quadrature points per face in @p dim dimensions. 397 */ 398 unsigned int n_q_points_face; 399 400 /** 401 * Stores the number of DoFs per face in @p dim dimensions. 402 */ 403 unsigned int dofs_per_component_on_face; 404 405 /** 406 * For nodal basis functions with nodes located at the boundary of the 407 * unit cell, face integrals that involve only the values of the shape 408 * functions (approximations of first derivatives in DG) do not need to 409 * load all degrees of freedom of the cell but rather only the degrees 410 * of freedom located on the face. While it would also be possible to 411 * compute these indices on the fly, we choose to simplify the code and 412 * store the indirect addressing in this class. 413 * 414 * The first table index runs through the faces of a cell, and the 415 * second runs through the nodal degrees of freedom of the face, using 416 * @p dofs_per_face entries. 417 * 418 * The indices stored in this member variable are as follows. Consider 419 * for example a 2D element of degree 3 with the following degrees of 420 * freedom in lexicographic numbering: 421 * @code 422 * 12 13 14 15 423 * 8 9 10 11 424 * 4 5 6 7 425 * 0 1 2 3 426 * @endcode 427 * 428 * The first row stores the indices on the face with index 0, i.e., the 429 * numbers <code>0, 4, 8, 12</code>, the second row holds the indices 430 * <code>3, 7, 11, 15</code> for face 1, the third row holds the indices 431 * <code>0, 1, 2, 3</code> for face 2, and the last (fourth) row holds 432 * the indices <code>12, 13, 14, 15</code>. Similarly, the indices are 433 * stored in 3D. (Note that the y faces in 3D use indices reversed in 434 * terms of the lexicographic numbers due to the orientation of the 435 * coordinate system.) 436 * 437 * @note This object is only filled in case @p nodal_at_cell_boundaries 438 * evaluates to @p true. 439 */ 440 dealii::Table<2, unsigned int> face_to_cell_index_nodal; 441 442 /** 443 * The @p face_to_cell_index_nodal provides a shortcut for the 444 * evaluation of values on the faces. For Hermite-type basis functions, 445 * one can go one step further and also use shortcuts to get derivatives 446 * more cheaply where only two layers of degrees of freedom contribute 447 * to the derivative on the face. In the lexicographic ordering, the 448 * respective indices is in the next "layer" of degrees of freedom as 449 * compared to the nodal values. This array stores the indirect 450 * addressing of both the values and the gradient. 451 * 452 * The first table index runs through the faces of a cell, and the 453 * second runs through the pairs of the nodal degrees of freedom of the 454 * face and the derivatives, using <code>2*dofs_per_face</code> entries. 455 * 456 * The indices stored in this member variable are as follows. Consider 457 * for example a 2D element of degree 3 with the following degrees of 458 * freedom in lexicographic numbering: 459 * @code 460 * 20 21 22 23 24 461 * 15 16 17 18 19 462 * 10 11 12 13 14 463 * 5 6 7 8 9 464 * 0 1 2 3 4 465 * @endcode 466 * 467 * The first row stores the indices for values and gradients on the face 468 * with index 0, i.e., the numbers <code>0, 1, 5, 6, 10, 11, 15, 16, 20, 469 * 21</code>, the second row holds the indices <code>4, 3, 9, 8, 14, 13, 470 * 19, 18, 24, 23</code> for face 1, the third row holds the indices 471 * <code>0, 5, 1, 6, 2, 7, 3, 8, 4, 9</code> for face 2, and the last 472 * (fourth) row holds the indices <code>20, 15, 21, 16, 22, 17, 23, 18, 473 * 24, 19</code>. Similarly, the indices are stored in 3D. (Note that 474 * the y faces in 3D use indices reversed in terms of the lexicographic 475 * numbers due to the orientation of the coordinate system.) 476 * 477 * @note This object is only filled in case @p element_type evaluates to 478 * @p tensor_symmetric_hermite. 479 */ 480 dealii::Table<2, unsigned int> face_to_cell_index_hermite; 481 482 /** 483 * For degrees on faces, the basis functions are not 484 * in the correct order if a face is not in the standard orientation 485 * to a given element. This data structure is used to re-order the 486 * basis functions to represent the correct order. 487 */ 488 dealii::Table<2, unsigned int> face_orientations; 489 490 private: 491 /** 492 * Check whether we have symmetries in the shape values. In that case, 493 * also fill the shape_???_eo fields. 494 */ 495 bool 496 check_1d_shapes_symmetric( 497 UnivariateShapeData<Number> &univariate_shape_data); 498 499 /** 500 * Check whether symmetric 1D basis functions are such that the shape 501 * values form a diagonal matrix, i.e., the nodal points are collocated 502 * with the quadrature points. This allows for specialized algorithms 503 * that save some operations in the evaluation. 504 */ 505 bool 506 check_1d_shapes_collocation( 507 const UnivariateShapeData<Number> &univariate_shape_data) const; 508 }; 509 510 511 512 // ------------------------------------------ inline functions 513 514 template <typename Number> 515 template <int dim> ShapeInfo(const Quadrature<1> & quad,const FiniteElement<dim> & fe_in,const unsigned int base_element_number)516 inline ShapeInfo<Number>::ShapeInfo(const Quadrature<1> & quad, 517 const FiniteElement<dim> &fe_in, 518 const unsigned int base_element_number) 519 : element_type(tensor_general) 520 , n_dimensions(0) 521 , n_components(0) 522 , n_q_points(0) 523 , dofs_per_component_on_cell(0) 524 , n_q_points_face(0) 525 , dofs_per_component_on_face(0) 526 { 527 reinit(quad, fe_in, base_element_number); 528 } 529 530 template <typename Number> 531 inline const UnivariateShapeData<Number> & get_shape_data(const unsigned int dimension,const unsigned int component)532 ShapeInfo<Number>::get_shape_data(const unsigned int dimension, 533 const unsigned int component) const 534 { 535 AssertDimension(n_dimensions, data_access.size(0)); 536 AssertDimension(n_components, data_access.size(1)); 537 AssertIndexRange(dimension, n_dimensions); 538 AssertIndexRange(component, n_components); 539 return *(data_access(dimension, component)); 540 } 541 542 } // end of namespace MatrixFreeFunctions 543 544 } // end of namespace internal 545 546 DEAL_II_NAMESPACE_CLOSE 547 548 #endif 549