1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 1998 - 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_update_flags_h 17 #define dealii_fe_update_flags_h 18 19 20 #include <deal.II/base/config.h> 21 22 #include <deal.II/base/derivative_form.h> 23 #include <deal.II/base/point.h> 24 #include <deal.II/base/table.h> 25 #include <deal.II/base/tensor.h> 26 27 #include <vector> 28 29 30 DEAL_II_NAMESPACE_OPEN 31 32 // Forward declaration 33 #ifndef DOXYGEN 34 template <int, int> 35 class FiniteElement; 36 #endif 37 38 /*!@addtogroup feaccess */ 39 /*@{*/ 40 41 /** 42 * The enum type given to the constructors of FEValues, FEFaceValues and 43 * FESubfaceValues, telling those objects which data will be needed on each 44 * mesh cell. 45 * 46 * Selecting these flags in a restrictive way is crucial for the efficiency of 47 * FEValues::reinit(), FEFaceValues::reinit() and FESubfaceValues::reinit(). 48 * Therefore, only the flags actually needed should be selected. It is the 49 * responsibility of the involved Mapping and FiniteElement to add additional 50 * flags according to their own requirements. For instance, most finite 51 * elements will add #update_covariant_transformation if #update_gradients is 52 * selected. By default, all flags are off, i.e. no reinitialization will be 53 * done. 54 * 55 * You can select more than one flag by concatenation using the bitwise or 56 * operator|(UpdateFlags,UpdateFlags). 57 * 58 * <h3>Use of these flags flags</h3> 59 * 60 * More information on the use of this type both in user code as well as 61 * internally can be found in the documentation modules on 62 * @ref UpdateFlags "The interplay of UpdateFlags, Mapping, and FiniteElement in FEValues" 63 * and 64 * @ref FE_vs_Mapping_vs_FEValues "How Mapping, FiniteElement, and FEValues work together". 65 */ 66 enum UpdateFlags 67 { 68 //! No update 69 update_default = 0, 70 //! Shape function values 71 /** 72 * Compute the values of the shape functions at the quadrature points on the 73 * real space cell. For the usual Lagrange elements, these values are equal 74 * to the values of the shape functions at the quadrature points on the unit 75 * cell, but they are different for more complicated elements, such as 76 * FE_RaviartThomas elements. 77 */ 78 update_values = 0x0001, 79 //! Shape function gradients 80 /** 81 * Compute the gradients of the shape functions in coordinates of the real 82 * cell. 83 */ 84 update_gradients = 0x0002, 85 //! Second derivatives of shape functions 86 /** 87 * Compute the second derivatives of the shape functions in coordinates of 88 * the real cell. 89 */ 90 update_hessians = 0x0004, 91 //! Third derivatives of shape functions 92 /** 93 * Compute the third derivatives of the shape functions in coordinates of 94 * the real cell 95 */ 96 update_3rd_derivatives = 0x0008, 97 //! Outer normal vector, not normalized 98 /** 99 * Vector product of tangential vectors, yielding a normal vector with a 100 * length corresponding to the surface element; may be more efficient than 101 * computing both. 102 */ 103 update_boundary_forms = 0x0010, 104 //! Transformed quadrature points 105 /** 106 * Compute the quadrature points location in real cell coordinates. 107 * 108 * FEValues objects take the quadrature point locations on the 109 * reference cell as an argument of the constructor (via the 110 * Quadrature object). For most finite elements, knowing the 111 * location of quadrature points on the reference cell is all that 112 * is necessary to evaluate shape functions, evaluate the mapping, 113 * and other things. On the other hand, if you want to evaluate a 114 * right hand side function $f(\mathbf x_q)$ at quadrature point 115 * locations $\mathbf x_q$ on the real cell, you need to pass this 116 * flag to the FEValues constructor to make sure you can later 117 * access them. 118 * 119 * In the context of DataPostprocessor, 120 * DataPostprocessorInputs::CommonInputs::evaluation_points will be updated. 121 */ 122 update_quadrature_points = 0x0020, 123 //! Transformed quadrature weights 124 /** 125 * Compute the quadrature weights on the real cell, i.e. the weights of the 126 * quadrature rule multiplied with the determinant of the Jacobian of the 127 * transformation from reference to real cell. 128 */ 129 update_JxW_values = 0x0040, 130 //! Normal vectors 131 /** 132 * Compute the normal vectors, either for a face or for a cell of 133 * codimension one. Setting this flag for any other object will raise an 134 * error. 135 */ 136 update_normal_vectors = 0x0080, 137 //! Volume element 138 /** 139 * Compute the Jacobian of the transformation from the reference cell to the 140 * real cell. 141 */ 142 update_jacobians = 0x0100, 143 //! Gradient of volume element 144 /** 145 * Compute the derivatives of the Jacobian of the transformation. 146 */ 147 update_jacobian_grads = 0x0200, 148 //! Volume element 149 /** 150 * Compute the inverse Jacobian of the transformation from the reference 151 * cell to the real cell. 152 */ 153 update_inverse_jacobians = 0x0400, 154 //! Covariant transformation 155 /** 156 * Compute all values the Mapping needs to perform a contravariant 157 * transformation of vectors. For special mappings like MappingCartesian 158 * this may be simpler than #update_inverse_jacobians. 159 */ 160 update_covariant_transformation = 0x0800, 161 //! Contravariant transformation 162 /** 163 * Compute all values the Mapping needs to perform a contravariant 164 * transformation of vectors. For special mappings like MappingCartesian 165 * this may be simpler than #update_jacobians. 166 */ 167 update_contravariant_transformation = 0x1000, 168 //! Shape function values of transformation 169 /** 170 * Compute the shape function values of the transformation defined by the 171 * Mapping. 172 */ 173 update_transformation_values = 0x2000, 174 //! Shape function gradients of transformation 175 /** 176 * Compute the shape function gradients of the transformation defined by the 177 * Mapping. 178 */ 179 update_transformation_gradients = 0x4000, 180 //! Determinant of the Jacobian 181 /** 182 * Compute the volume element in each quadrature point. 183 */ 184 update_volume_elements = 0x10000, 185 /** 186 * Compute the derivatives of the Jacobian of the transformation pushed 187 * forward to the real cell coordinates. 188 */ 189 update_jacobian_pushed_forward_grads = 0x100000, 190 /** 191 * Compute the second derivatives of the Jacobian of the transformation. 192 */ 193 update_jacobian_2nd_derivatives = 0x200000, 194 /** 195 * Compute the second derivatives of the Jacobian of the transformation 196 * pushed forward to the real cell coordinates. 197 */ 198 update_jacobian_pushed_forward_2nd_derivatives = 0x400000, 199 /** 200 * Compute the third derivatives of the Jacobian of the transformation. 201 */ 202 update_jacobian_3rd_derivatives = 0x800000, 203 /** 204 * Compute the third derivatives of the Jacobian of the transformation 205 * pushed forward to the real cell coordinates. 206 */ 207 update_jacobian_pushed_forward_3rd_derivatives = 0x1000000, 208 //! Values needed for Piola transform 209 /** 210 * Combination of the flags needed for Piola transform of Hdiv elements. 211 */ 212 update_piola = update_volume_elements | update_contravariant_transformation, 213 /** 214 * Combination of the flags that require a mapping calculation 215 */ 216 update_mapping = 217 // Direct data 218 update_quadrature_points | update_JxW_values | update_jacobians | 219 update_jacobian_grads | update_jacobian_pushed_forward_grads | 220 update_jacobian_2nd_derivatives | 221 update_jacobian_pushed_forward_2nd_derivatives | 222 update_jacobian_3rd_derivatives | 223 update_jacobian_pushed_forward_3rd_derivatives | update_inverse_jacobians | 224 update_boundary_forms | update_normal_vectors | 225 // Transformation dependence 226 update_covariant_transformation | update_contravariant_transformation | 227 update_transformation_values | update_transformation_gradients | 228 // Volume data 229 update_volume_elements 230 }; 231 232 233 /** 234 * Output operator which outputs update flags as a set of or'd text values. 235 * 236 * @ref UpdateFlags 237 */ 238 template <class StreamType> 239 inline StreamType & 240 operator<<(StreamType &s, const UpdateFlags u) 241 { 242 s << " UpdateFlags|"; 243 if (u & update_values) 244 s << "values|"; 245 if (u & update_gradients) 246 s << "gradients|"; 247 if (u & update_hessians) 248 s << "hessians|"; 249 if (u & update_3rd_derivatives) 250 s << "3rd_derivatives|"; 251 if (u & update_quadrature_points) 252 s << "quadrature_points|"; 253 if (u & update_JxW_values) 254 s << "JxW_values|"; 255 if (u & update_normal_vectors) 256 s << "normal_vectors|"; 257 if (u & update_jacobians) 258 s << "jacobians|"; 259 if (u & update_inverse_jacobians) 260 s << "inverse_jacobians|"; 261 if (u & update_jacobian_grads) 262 s << "jacobian_grads|"; 263 if (u & update_covariant_transformation) 264 s << "covariant_transformation|"; 265 if (u & update_contravariant_transformation) 266 s << "contravariant_transformation|"; 267 if (u & update_transformation_values) 268 s << "transformation_values|"; 269 if (u & update_transformation_gradients) 270 s << "transformation_gradients|"; 271 if (u & update_jacobian_pushed_forward_grads) 272 s << "jacobian_pushed_forward_grads|"; 273 if (u & update_jacobian_2nd_derivatives) 274 s << "jacobian_2nd_derivatives|"; 275 if (u & update_jacobian_pushed_forward_2nd_derivatives) 276 s << "jacobian_pushed_forward_2nd_derivatives|"; 277 if (u & update_jacobian_3rd_derivatives) 278 s << "jacobian_3rd_derivatives|"; 279 if (u & update_jacobian_pushed_forward_3rd_derivatives) 280 s << "jacobian_pushed_forward_3rd_derivatives|"; 281 282 // TODO: check that 'u' really only has the flags set that are handled above 283 return s; 284 } 285 286 287 /** 288 * Global operator which returns an object in which all bits are set which are 289 * either set in the first or the second argument. This operator exists since 290 * if it did not then the result of the bit-or <tt>operator |</tt> would be an 291 * integer which would in turn trigger a compiler warning when we tried to 292 * assign it to an object of type UpdateFlags. 293 * 294 * @ref UpdateFlags 295 */ 296 inline UpdateFlags 297 operator|(const UpdateFlags f1, const UpdateFlags f2) 298 { 299 return static_cast<UpdateFlags>(static_cast<unsigned int>(f1) | 300 static_cast<unsigned int>(f2)); 301 } 302 303 304 305 /** 306 * Global operator which sets the bits from the second argument also in the 307 * first one. 308 * 309 * @ref UpdateFlags 310 */ 311 inline UpdateFlags & 312 operator|=(UpdateFlags &f1, const UpdateFlags f2) 313 { 314 f1 = f1 | f2; 315 return f1; 316 } 317 318 319 /** 320 * Global operator which returns an object in which all bits are set which are 321 * set in the first as well as the second argument. This operator exists since 322 * if it did not then the result of the bit-and <tt>operator &</tt> would be 323 * an integer which would in turn trigger a compiler warning when we tried to 324 * assign it to an object of type UpdateFlags. 325 * 326 * @ref UpdateFlags 327 */ 328 inline UpdateFlags operator&(const UpdateFlags f1, const UpdateFlags f2) 329 { 330 return static_cast<UpdateFlags>(static_cast<unsigned int>(f1) & 331 static_cast<unsigned int>(f2)); 332 } 333 334 335 /** 336 * Global operator which clears all the bits in the first argument if they are 337 * not also set in the second argument. 338 * 339 * @ref UpdateFlags 340 */ 341 inline UpdateFlags & 342 operator&=(UpdateFlags &f1, const UpdateFlags f2) 343 { 344 f1 = f1 & f2; 345 return f1; 346 } 347 348 349 350 /** 351 * This enum definition is used for storing similarities of the current cell 352 * to the previously visited cell. This information is used for reusing data 353 * when calling the method FEValues::reinit() (like derivatives, which do not 354 * change if one cell is just a translation of the previous). Currently, this 355 * variable does only recognize a translation and an inverted translation (if 356 * dim<spacedim). However, this concept makes it easy to add additional states 357 * to be detected in FEValues/FEFaceValues for making use of these 358 * similarities as well. 359 */ 360 namespace CellSimilarity 361 { 362 enum Similarity 363 { 364 /** 365 * The cells differ by something besides a translation or inverted 366 * translations. 367 */ 368 none, 369 /** 370 * The cells differ by a translation. 371 */ 372 translation, 373 /** 374 * The cells differ by an inverted translation. 375 */ 376 inverted_translation, 377 /** 378 * The next cell is not valid. 379 */ 380 invalid_next_cell 381 }; 382 } 383 384 385 namespace internal 386 { 387 namespace FEValuesImplementation 388 { 389 /** 390 * A class that stores all of the mapping related data used in 391 * dealii::FEValues, dealii::FEFaceValues, and dealii::FESubfaceValues 392 * objects. Objects of this kind will be given as <i>output</i> argument 393 * when dealii::FEValues::reinit() calls Mapping::fill_fe_values() for a 394 * given cell, face, or subface. 395 * 396 * The data herein will then be provided as <i>input</i> argument in the 397 * following call to FiniteElement::fill_fe_values(). 398 * 399 * @ingroup feaccess 400 */ 401 template <int dim, int spacedim = dim> 402 class MappingRelatedData 403 { 404 public: 405 /** 406 * Initialize all vectors to correct size. 407 */ 408 void 409 initialize(const unsigned int n_quadrature_points, 410 const UpdateFlags flags); 411 412 /** 413 * Compute and return an estimate for the memory consumption (in bytes) 414 * of this object. 415 */ 416 std::size_t 417 memory_consumption() const; 418 419 /** 420 * Store an array of weights times the Jacobi determinant at the 421 * quadrature points. This function is reset each time reinit() is 422 * called. The Jacobi determinant is actually the reciprocal value of 423 * the Jacobi matrices stored in this class, see the general 424 * documentation of this class for more information. 425 * 426 * However, if this object refers to an FEFaceValues or FESubfaceValues 427 * object, then the JxW_values correspond to the Jacobian of the 428 * transformation of the face, not the cell, i.e. the dimensionality is 429 * that of a surface measure, not of a volume measure. In this case, it 430 * is computed from the boundary forms, rather than the Jacobian matrix. 431 */ 432 std::vector<double> JxW_values; 433 434 /** 435 * Array of the Jacobian matrices at the quadrature points. 436 */ 437 std::vector<DerivativeForm<1, dim, spacedim>> jacobians; 438 439 /** 440 * Array of the derivatives of the Jacobian matrices at the quadrature 441 * points. 442 */ 443 std::vector<DerivativeForm<2, dim, spacedim>> jacobian_grads; 444 445 /** 446 * Array of the inverse Jacobian matrices at the quadrature points. 447 */ 448 std::vector<DerivativeForm<1, spacedim, dim>> inverse_jacobians; 449 450 /** 451 * Array of the derivatives of the Jacobian matrices at the quadrature 452 * points, pushed forward to the real cell coordinates. 453 */ 454 std::vector<Tensor<3, spacedim>> jacobian_pushed_forward_grads; 455 456 /** 457 * Array of the second derivatives of the Jacobian matrices at the 458 * quadrature points. 459 */ 460 std::vector<DerivativeForm<3, dim, spacedim>> jacobian_2nd_derivatives; 461 462 /** 463 * Array of the second derivatives of the Jacobian matrices at the 464 * quadrature points, pushed forward to the real cell coordinates. 465 */ 466 std::vector<Tensor<4, spacedim>> jacobian_pushed_forward_2nd_derivatives; 467 468 /** 469 * Array of the third derivatives of the Jacobian matrices at the 470 * quadrature points. 471 */ 472 std::vector<DerivativeForm<4, dim, spacedim>> jacobian_3rd_derivatives; 473 474 /** 475 * Array of the third derivatives of the Jacobian matrices at the 476 * quadrature points, pushed forward to the real cell coordinates. 477 */ 478 std::vector<Tensor<5, spacedim>> jacobian_pushed_forward_3rd_derivatives; 479 480 /** 481 * Array of quadrature points. This array is set up upon calling 482 * reinit() and contains the quadrature points on the real element, 483 * rather than on the reference element. 484 */ 485 std::vector<Point<spacedim>> quadrature_points; 486 487 /** 488 * List of outward normal vectors at the quadrature points. 489 */ 490 std::vector<Tensor<1, spacedim>> normal_vectors; 491 492 /** 493 * List of boundary forms at the quadrature points. 494 */ 495 std::vector<Tensor<1, spacedim>> boundary_forms; 496 }; 497 498 499 /** 500 * A class that stores all of the shape function related data used in 501 * dealii::FEValues, dealii::FEFaceValues, and dealii::FESubfaceValues 502 * objects. Objects of this kind will be given as <i>output</i> argument 503 * when dealii::FEValues::reinit() calls FiniteElement::fill_fe_values(). 504 * 505 * @ingroup feaccess 506 */ 507 template <int dim, int spacedim = dim> 508 class FiniteElementRelatedData 509 { 510 public: 511 /** 512 * Initialize all vectors to correct size. 513 */ 514 void 515 initialize(const unsigned int n_quadrature_points, 516 const FiniteElement<dim, spacedim> &fe, 517 const UpdateFlags flags); 518 519 /** 520 * Compute and return an estimate for the memory consumption (in bytes) 521 * of this object. 522 */ 523 std::size_t 524 memory_consumption() const; 525 526 /** 527 * Storage type for shape values. Each row in the matrix denotes the 528 * values of a single shape function at the different points, columns 529 * are for a single point with the different shape functions. 530 * 531 * If a shape function has more than one non-zero component (in deal.II 532 * diction: it is non-primitive), then we allocate one row per non-zero 533 * component, and shift subsequent rows backward. Lookup of the correct 534 * row for a shape function is thus simple in case the entire finite 535 * element is primitive (i.e. all shape functions are primitive), since 536 * then the shape function number equals the row number. Otherwise, use 537 * the #shape_function_to_row_table array to get at the first row that 538 * belongs to this particular shape function, and navigate among all the 539 * rows for this shape function using the 540 * FiniteElement::get_nonzero_components() function which tells us which 541 * components are non-zero and thus have a row in the array presently 542 * under discussion. 543 */ 544 using ShapeVector = dealii::Table<2, double>; 545 546 /** 547 * Storage type for gradients. The layout of data is the same as for the 548 * #ShapeVector data type. 549 */ 550 using GradientVector = dealii::Table<2, Tensor<1, spacedim>>; 551 552 /** 553 * Likewise for second order derivatives. 554 */ 555 using HessianVector = dealii::Table<2, Tensor<2, spacedim>>; 556 557 /** 558 * And the same also applies to the third order derivatives. 559 */ 560 using ThirdDerivativeVector = dealii::Table<2, Tensor<3, spacedim>>; 561 562 /** 563 * Store the values of the shape functions at the quadrature points. See 564 * the description of the data type for the layout of the data in this 565 * field. 566 */ 567 ShapeVector shape_values; 568 569 /** 570 * Store the gradients of the shape functions at the quadrature points. 571 * See the description of the data type for the layout of the data in 572 * this field. 573 */ 574 GradientVector shape_gradients; 575 576 /** 577 * Store the 2nd derivatives of the shape functions at the quadrature 578 * points. See the description of the data type for the layout of the 579 * data in this field. 580 */ 581 HessianVector shape_hessians; 582 583 /** 584 * Store the 3nd derivatives of the shape functions at the quadrature 585 * points. See the description of the data type for the layout of the 586 * data in this field. 587 */ 588 ThirdDerivativeVector shape_3rd_derivatives; 589 590 /** 591 * When asked for the value (or gradient, or Hessian) of shape function 592 * i's c-th vector component, we need to look it up in the 593 * #shape_values, #shape_gradients and #shape_hessians arrays. The 594 * question is where in this array does the data for shape function i, 595 * component c reside. This is what this table answers. 596 * 597 * The format of the table is as follows: - It has dofs_per_cell times 598 * n_components entries. - The entry that corresponds to shape function 599 * i, component c is <code>i * n_components + c</code>. - The value 600 * stored at this position indicates the row in #shape_values and the 601 * other tables where the corresponding datum is stored for all the 602 * quadrature points. 603 * 604 * In the general, vector-valued context, the number of components is 605 * larger than one, but for a given shape function, not all vector 606 * components may be nonzero (e.g., if a shape function is primitive, 607 * then exactly one vector component is non-zero, while the others are 608 * all zero). For such zero components, #shape_values and friends do not 609 * have a row. Consequently, for vector components for which shape 610 * function i is zero, the entry in the current table is 611 * numbers::invalid_unsigned_int. 612 * 613 * On the other hand, the table is guaranteed to have at least one valid 614 * index for each shape function. In particular, for a primitive finite 615 * element, each shape function has exactly one nonzero component and so 616 * for each i, there is exactly one valid index within the range 617 * <code>[i*n_components, (i+1)*n_components)</code>. 618 */ 619 std::vector<unsigned int> shape_function_to_row_table; 620 }; 621 } // namespace FEValuesImplementation 622 } // namespace internal 623 624 625 /*@}*/ 626 627 628 629 DEAL_II_NAMESPACE_CLOSE 630 631 #endif 632