1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2001 - 2019 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 #ifndef dealii_fe_dgq_h 17 #define dealii_fe_dgq_h 18 19 #include <deal.II/base/config.h> 20 21 #include <deal.II/base/tensor_product_polynomials.h> 22 #include <deal.II/base/thread_management.h> 23 24 #include <deal.II/fe/fe_poly.h> 25 26 DEAL_II_NAMESPACE_OPEN 27 28 // Forward declarations 29 #ifndef DOXYGEN 30 template <int dim, int spacedim> 31 class MappingQ; 32 template <int dim> 33 class Quadrature; 34 #endif 35 36 /*!@addtogroup fe */ 37 /*@{*/ 38 39 /** 40 * Implementation of scalar, discontinuous tensor product elements based on 41 * equidistant support points. 42 * 43 * This is a discontinuous finite element based on tensor products of 44 * Lagrangian polynomials. The shape functions are Lagrangian interpolants of 45 * an equidistant grid of points on the unit cell. The points are numbered in 46 * lexicographical order, with <i>x</i> running fastest, then <i>y</i>, then 47 * <i>z</i> (if these coordinates are present for a given space dimension at 48 * all). For example, these are the node orderings for <tt>FE_DGQ(1)</tt> in 49 * 3d: 50 * @verbatim 51 * 6-------7 6-------7 52 * /| | / /| 53 * / | | / / | 54 * / | | / / | 55 * 4 | | 4-------5 | 56 * | 2-------3 | | 3 57 * | / / | | / 58 * | / / | | / 59 * |/ / | |/ 60 * 0-------1 0-------1 61 * @endverbatim 62 * and <tt>FE_DGQ(2)</tt>: 63 * @verbatim 64 * 24--25--26 24--25--26 65 * /| | / /| 66 * 21 | | 21 22 23 | 67 * / 15 16 17 / / 17 68 * 18 | | 18--19--20 | 69 * |12 6---7---8 | |14 8 70 * 9 / / 9 10 11 / 71 * | 3 4 5 | | 5 72 * |/ / | |/ 73 * 0---1---2 0---1---2 74 * @endverbatim 75 * with node 13 being placed in the interior of the hex. 76 * 77 * Note, however, that these are just the Lagrange interpolation points of the 78 * shape functions. Even though they may physically be on the boundary of the 79 * cell, they are logically in the interior since there are no continuity 80 * requirements for these shape functions across cell boundaries. While 81 * discontinuous, when restricted to a single cell the shape functions of this 82 * element are exactly the same as those of the FE_Q element where they are 83 * shown visually. 84 * 85 * <h3>Unit support point distribution and conditioning of interpolation</h3> 86 * 87 * When constructing an FE_DGQ element at polynomial degrees one or two, 88 * equidistant support points at 0 and 1 (linear case) or 0, 0.5, and 1 89 * (quadratic case) are used. The unit support or nodal points 90 * <i>x<sub>i</sub></i> are those points where the <i>j</i>th Lagrange 91 * polynomial satisfies the $\delta_{ij}$ property, i.e., where one polynomial 92 * is one and all the others are zero. For higher polynomial degrees, the 93 * support points are non-equidistant by default, and chosen to be the support 94 * points of the <tt>(degree+1)</tt>-order Gauss-Lobatto quadrature rule. This 95 * point distribution yields well-conditioned Lagrange interpolation at 96 * arbitrary polynomial degrees. By contrast, polynomials based on equidistant 97 * points get increasingly ill-conditioned as the polynomial degree 98 * increases. In interpolation, this effect is known as the Runge 99 * phenomenon. For Galerkin methods, the Runge phenomenon is typically not 100 * visible in the solution quality but rather in the condition number of the 101 * associated system matrices. For example, the elemental mass matrix of 102 * equidistant points at degree 10 has condition number 2.6e6, whereas the 103 * condition number for Gauss-Lobatto points is around 400. 104 * 105 * The Gauss-Lobatto points in 1D include the end points 0 and +1 of the unit 106 * interval. The interior points are shifted towards the end points, which 107 * gives a denser point distribution close to the element boundary. 108 */ 109 template <int dim, int spacedim = dim> 110 class FE_DGQ : public FE_Poly<dim, spacedim> 111 { 112 public: 113 /** 114 * Constructor for tensor product polynomials of degree <tt>p</tt>. The 115 * shape functions created using this constructor correspond to Lagrange 116 * interpolation polynomials for Gauss-Lobatto support (node) points in each 117 * coordinate direction. 118 */ 119 FE_DGQ(const unsigned int p); 120 121 /** 122 * Return a string that uniquely identifies a finite element. This class 123 * returns <tt>FE_DGQ<dim>(degree)</tt>, with <tt>dim</tt> and 124 * <tt>degree</tt> replaced by appropriate values. 125 */ 126 virtual std::string 127 get_name() const override; 128 129 /** 130 * Return the matrix interpolating from the given finite element to the 131 * present one. The size of the matrix is then @p dofs_per_cell times 132 * <tt>source.n_dofs_per_cell()</tt>. 133 * 134 * These matrices are only available if the source element is also a @p 135 * FE_DGQ element. Otherwise, an exception of type 136 * FiniteElement<dim>::ExcInterpolationNotImplemented is thrown. 137 */ 138 virtual void 139 get_interpolation_matrix(const FiniteElement<dim, spacedim> &source, 140 FullMatrix<double> &matrix) const override; 141 142 /** 143 * Return the matrix interpolating from a face of one element to the face 144 * of the neighboring element. The size of the matrix is then 145 * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. 146 * 147 * Derived elements will have to implement this function. They may only 148 * provide interpolation matrices for certain source finite elements, for 149 * example those from the same family. If they don't implement interpolation 150 * from a given element, then they must throw an exception of type 151 * FiniteElement<dim>::ExcInterpolationNotImplemented. 152 */ 153 virtual void 154 get_face_interpolation_matrix(const FiniteElement<dim, spacedim> &source, 155 FullMatrix<double> & matrix, 156 const unsigned int face_no = 0) const override; 157 158 /** 159 * Return the matrix interpolating from a face of one element to the face 160 * of the neighboring element. The size of the matrix is then 161 * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. 162 * 163 * Derived elements will have to implement this function. They may only 164 * provide interpolation matrices for certain source finite elements, for 165 * example those from the same family. If they don't implement interpolation 166 * from a given element, then they must throw an exception of type 167 * FiniteElement<dim>::ExcInterpolationNotImplemented. 168 */ 169 virtual void 170 get_subface_interpolation_matrix( 171 const FiniteElement<dim, spacedim> &source, 172 const unsigned int subface, 173 FullMatrix<double> & matrix, 174 const unsigned int face_no = 0) const override; 175 176 /** 177 * Projection from a fine grid space onto a coarse grid space. Overrides the 178 * respective method in FiniteElement, implementing lazy evaluation 179 * (initialize when requested). 180 * 181 * If this projection operator is associated with a matrix @p P, then the 182 * restriction of this matrix @p P_i to a single child cell is returned 183 * here. 184 * 185 * The matrix @p P is the concatenation or the sum of the cell matrices @p 186 * P_i, depending on the #restriction_is_additive_flags. This distinguishes 187 * interpolation (concatenation) and projection with respect to scalar 188 * products (summation). 189 * 190 * Row and column indices are related to coarse grid and fine grid spaces, 191 * respectively, consistent with the definition of the associated operator. 192 */ 193 virtual const FullMatrix<double> & 194 get_restriction_matrix( 195 const unsigned int child, 196 const RefinementCase<dim> &refinement_case = 197 RefinementCase<dim>::isotropic_refinement) const override; 198 199 /** 200 * Embedding matrix between grids. Overrides the respective method in 201 * FiniteElement, implementing lazy evaluation (initialize when queried). 202 * 203 * The identity operator from a coarse grid space into a fine grid space is 204 * associated with a matrix @p P. The restriction of this matrix @p P_i to a 205 * single child cell is returned here. 206 * 207 * The matrix @p P is the concatenation, not the sum of the cell matrices @p 208 * P_i. That is, if the same non-zero entry <tt>j,k</tt> exists in two 209 * different child matrices @p P_i, the value should be the same in both 210 * matrices and it is copied into the matrix @p P only once. 211 * 212 * Row and column indices are related to fine grid and coarse grid spaces, 213 * respectively, consistent with the definition of the associated operator. 214 * 215 * These matrices are used by routines assembling the prolongation matrix 216 * for multi-level methods. Upon assembling the transfer matrix between 217 * cells using this matrix array, zero elements in the prolongation matrix 218 * are discarded and will not fill up the transfer matrix. 219 */ 220 virtual const FullMatrix<double> & 221 get_prolongation_matrix( 222 const unsigned int child, 223 const RefinementCase<dim> &refinement_case = 224 RefinementCase<dim>::isotropic_refinement) const override; 225 226 /** 227 * @name Functions to support hp 228 * @{ 229 */ 230 231 /** 232 * If, on a vertex, several finite elements are active, the hp code first 233 * assigns the degrees of freedom of each of these FEs different global 234 * indices. It then calls this function to find out which of them should get 235 * identical values, and consequently can receive the same global DoF index. 236 * This function therefore returns a list of identities between DoFs of the 237 * present finite element object with the DoFs of @p fe_other, which is a 238 * reference to a finite element object representing one of the other finite 239 * elements active on this particular vertex. The function computes which of 240 * the degrees of freedom of the two finite element objects are equivalent, 241 * both numbered between zero and the corresponding value of 242 * n_dofs_per_vertex() of the two finite elements. The first index of each 243 * pair denotes one of the vertex dofs of the present element, whereas the 244 * second is the corresponding index of the other finite element. 245 * 246 * This being a discontinuous element, the set of such constraints is of 247 * course empty. 248 */ 249 virtual std::vector<std::pair<unsigned int, unsigned int>> 250 hp_vertex_dof_identities( 251 const FiniteElement<dim, spacedim> &fe_other) const override; 252 253 /** 254 * Same as hp_vertex_dof_indices(), except that the function treats degrees 255 * of freedom on lines. 256 * 257 * This being a discontinuous element, the set of such constraints is of 258 * course empty. 259 */ 260 virtual std::vector<std::pair<unsigned int, unsigned int>> 261 hp_line_dof_identities( 262 const FiniteElement<dim, spacedim> &fe_other) const override; 263 264 /** 265 * Same as hp_vertex_dof_indices(), except that the function treats degrees 266 * of freedom on quads. 267 * 268 * This being a discontinuous element, the set of such constraints is of 269 * course empty. 270 */ 271 virtual std::vector<std::pair<unsigned int, unsigned int>> 272 hp_quad_dof_identities(const FiniteElement<dim, spacedim> &fe_other, 273 const unsigned int face_no = 0) const override; 274 275 /** 276 * Return whether this element implements its hanging node constraints in 277 * the new way, which has to be used to make elements "hp compatible". 278 * 279 * For the FE_DGQ class the result is always true (independent of the degree 280 * of the element), as it has no hanging nodes (being a discontinuous 281 * element). 282 */ 283 virtual bool 284 hp_constraints_are_implemented() const override; 285 286 /** 287 * @copydoc FiniteElement::compare_for_domination() 288 */ 289 virtual FiniteElementDomination::Domination 290 compare_for_domination(const FiniteElement<dim, spacedim> &fe_other, 291 const unsigned int codim = 0) const override final; 292 293 /** 294 * @} 295 */ 296 297 /** 298 * This function returns @p true, if the shape function @p shape_index has 299 * non-zero function values somewhere on the face @p face_index. 300 */ 301 virtual bool 302 has_support_on_face(const unsigned int shape_index, 303 const unsigned int face_index) const override; 304 305 /** 306 * Return a list of constant modes of the element. For this element, it 307 * simply returns one row with all entries set to true. 308 */ 309 virtual std::pair<Table<2, bool>, std::vector<unsigned int>> 310 get_constant_modes() const override; 311 312 /** 313 * Implementation of the corresponding function in the FiniteElement 314 * class. Since the current element is interpolatory, the nodal 315 * values are exactly the support point values. Furthermore, since 316 * the current element is scalar, the support point values need to 317 * be vectors of length 1. 318 */ 319 virtual void 320 convert_generalized_support_point_values_to_dof_values( 321 const std::vector<Vector<double>> &support_point_values, 322 std::vector<double> & nodal_values) const override; 323 324 virtual std::unique_ptr<FiniteElement<dim, spacedim>> 325 clone() const override; 326 327 protected: 328 /** 329 * Constructor for tensor product polynomials based on an arbitrary vector 330 * of polynomials. This constructor is used in derived classes to construct 331 * e.g. elements with arbitrary nodes or elements based on Legendre 332 * polynomials. 333 * 334 * The degree of these polynomials is <tt>polynomials.size()-1</tt>. 335 */ 336 FE_DGQ(const std::vector<Polynomials::Polynomial<double>> &polynomials); 337 338 private: 339 /** 340 * Only for internal use. Its full name is @p get_dofs_per_object_vector 341 * function and it creates the @p dofs_per_object vector that is needed 342 * within the constructor to be passed to the constructor of @p 343 * FiniteElementData. 344 */ 345 static std::vector<unsigned int> 346 get_dpo_vector(const unsigned int degree); 347 348 /** 349 * Compute renumbering for rotation of degrees of freedom. 350 * 351 * This function rotates a tensor product numbering of degrees of 352 * freedom by 90 degrees. It is used to compute the transfer 353 * matrices of the children by using only the matrix for the first 354 * child. 355 * 356 * The direction parameter determines the type of rotation. It is one 357 * character of @p xXyYzZ. The character determines the axis of rotation, 358 * case determines the direction. Lower case is counter-clockwise seen in 359 * direction of the axis. 360 * 361 * Since rotation around the y-axis is not used, it is not implemented 362 * either. 363 */ 364 void 365 rotate_indices(std::vector<unsigned int> &indices, 366 const char direction) const; 367 368 /* 369 * Mutex for protecting initialization of restriction and embedding matrix. 370 */ 371 mutable Threads::Mutex mutex; 372 373 // Allow access from other dimensions. 374 template <int dim1, int spacedim1> 375 friend class FE_DGQ; 376 377 // Allow @p MappingQ class to access to build_renumbering function. 378 template <int dim1, int spacedim1> 379 friend class MappingQ; 380 }; 381 382 383 384 /** 385 * Implementation of scalar, discontinuous tensor product elements based on 386 * Lagrange polynomials with arbitrary nodes. The primary purpose of this 387 * class is to provide an element for which the mass matrix can be made 388 * diagonal by choosing basis functions that are not either zero or one at the 389 * vertices of the cell, but instead are zero or one at a given set of 390 * quadrature points. If this set of quadrature points is then also used in 391 * integrating the mass matrix, then it will be diagonal. The number of 392 * quadrature points automatically determines the polynomial degree chosen for 393 * this element. The typical applications are the Gauss quadrature or the 394 * Gauss-Lobatto quadrature (provided through the base class). 395 * 396 * See the base class documentation in FE_DGQ for details. 397 */ 398 template <int dim, int spacedim = dim> 399 class FE_DGQArbitraryNodes : public FE_DGQ<dim, spacedim> 400 { 401 public: 402 /** 403 * Constructor for tensor product polynomials based on Polynomials::Lagrange 404 * interpolation of the support points in the quadrature rule 405 * <tt>points</tt>. The degree of these polynomials is 406 * <tt>points.size()-1</tt>. 407 */ 408 FE_DGQArbitraryNodes(const Quadrature<1> &points); 409 410 /** 411 * Return a string that uniquely identifies a finite element. This class 412 * returns <tt>FE_DGQArbitraryNodes<dim>(degree)</tt>, with <tt>dim</tt> and 413 * <tt>degree</tt> replaced by appropriate values. 414 */ 415 virtual std::string 416 get_name() const override; 417 418 /** 419 * Implementation of the corresponding function in the FiniteElement 420 * class. Since the current element is interpolatory, the nodal 421 * values are exactly the support point values. Furthermore, since 422 * the current element is scalar, the support point values need to 423 * be vectors of length 1. 424 */ 425 virtual void 426 convert_generalized_support_point_values_to_dof_values( 427 const std::vector<Vector<double>> &support_point_values, 428 std::vector<double> & nodal_values) const override; 429 virtual std::unique_ptr<FiniteElement<dim, spacedim>> 430 clone() const override; 431 }; 432 433 434 435 /** 436 * Implementation of scalar, discontinuous tensor product elements based on 437 * Legendre polynomials, described by the tensor product of the polynomial 438 * space Polynomials::Legendre. As opposed to the basic FE_DGQ element, these 439 * elements are not interpolatory and no support points are defined. 440 * 441 * See the base class documentation in FE_DGQ for details. 442 */ 443 template <int dim, int spacedim = dim> 444 class FE_DGQLegendre : public FE_DGQ<dim, spacedim> 445 { 446 public: 447 /** 448 * Constructor for tensor product polynomials based on Polynomials::Legendre 449 * interpolation. 450 */ 451 FE_DGQLegendre(const unsigned int degree); 452 453 /** 454 * Return a list of constant modes of the element. For the Legendre basis, 455 * it returns one row where the first element (corresponding to the constant 456 * mode) is set to true and all other elements are set to false. 457 */ 458 virtual std::pair<Table<2, bool>, std::vector<unsigned int>> 459 get_constant_modes() const override; 460 461 /** 462 * Return a string that uniquely identifies a finite element. This class 463 * returns <tt>FE_DGQLegendre<dim>(degree)</tt> with <tt>dim</tt> and 464 * <tt>degree</tt> replaced by the values given by the template parameter 465 * and the argument passed to the constructor, respectively. 466 */ 467 virtual std::string 468 get_name() const override; 469 470 virtual std::unique_ptr<FiniteElement<dim, spacedim>> 471 clone() const override; 472 }; 473 474 475 476 /** 477 * Implementation of scalar, discontinuous tensor product elements based on 478 * Hermite-like polynomials, described by the polynomial space 479 * Polynomials::HermiteLikeInterpolation. As opposed to the basic FE_DGQ 480 * element, these elements are not interpolatory and no support points are 481 * defined. 482 * 483 * Note that Hermite polynomials are only available for degrees larger or 484 * equal to three, and thus the beneficial properties of 485 * Polynomials::HermiteLikeInterpolation with only two basis functions having 486 * a non-trivial value or derivative on a face per dimension is only present 487 * for higher degrees. To facilitate usage also for degrees zero to two, a 488 * usual Lagrange basis is constructed by this class. 489 * 490 * See the base class documentation in FE_DGQ for details. 491 */ 492 template <int dim, int spacedim = dim> 493 class FE_DGQHermite : public FE_DGQ<dim, spacedim> 494 { 495 public: 496 /** 497 * Constructor for tensor product polynomials based on 498 * Polynomials::HermiteLikeInterpolation. 499 */ 500 FE_DGQHermite(const unsigned int degree); 501 502 /** 503 * Return a string that uniquely identifies a finite element. This class 504 * returns <tt>FE_DGQHermite<dim>(degree)</tt>, with <tt>dim</tt> and 505 * <tt>degree</tt> replaced by the values given by the template parameter 506 * and the argument passed to the constructor, respectively. 507 */ 508 virtual std::string 509 get_name() const override; 510 511 virtual std::unique_ptr<FiniteElement<dim, spacedim>> 512 clone() const override; 513 }; 514 515 516 /*@}*/ 517 518 DEAL_II_NAMESPACE_CLOSE 519 520 #endif 521