1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 1998 - 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 #ifndef dealii_fe_values_h 17 #define dealii_fe_values_h 18 19 20 #include <deal.II/base/config.h> 21 22 #include <deal.II/base/derivative_form.h> 23 #include <deal.II/base/exceptions.h> 24 #include <deal.II/base/point.h> 25 #include <deal.II/base/quadrature.h> 26 #include <deal.II/base/std_cxx20/iota_view.h> 27 #include <deal.II/base/subscriptor.h> 28 #include <deal.II/base/symmetric_tensor.h> 29 30 #include <deal.II/dofs/dof_accessor.h> 31 #include <deal.II/dofs/dof_handler.h> 32 33 #include <deal.II/fe/fe.h> 34 #include <deal.II/fe/fe_update_flags.h> 35 #include <deal.II/fe/fe_values_extractors.h> 36 #include <deal.II/fe/mapping.h> 37 38 #include <deal.II/grid/tria.h> 39 #include <deal.II/grid/tria_iterator.h> 40 41 #include <algorithm> 42 #include <memory> 43 #include <type_traits> 44 45 46 // dummy include in order to have the 47 // definition of PetscScalar available 48 // without including other PETSc stuff 49 #ifdef DEAL_II_WITH_PETSC 50 # include <petsc.h> 51 #endif 52 53 DEAL_II_NAMESPACE_OPEN 54 55 // Forward declaration 56 #ifndef DOXYGEN 57 template <int dim, int spacedim = dim> 58 class FEValuesBase; 59 #endif 60 61 namespace internal 62 { 63 /** 64 * A class whose specialization is used to define what type the curl of a 65 * vector valued function corresponds to. 66 */ 67 template <int dim, class NumberType = double> 68 struct CurlType; 69 70 /** 71 * A class whose specialization is used to define what type the curl of a 72 * vector valued function corresponds to. 73 * 74 * In 1d, the curl is a scalar. 75 */ 76 template <class NumberType> 77 struct CurlType<1, NumberType> 78 { 79 using type = Tensor<1, 1, NumberType>; 80 }; 81 82 /** 83 * A class whose specialization is used to define what type the curl of a 84 * vector valued function corresponds to. 85 * 86 * In 2d, the curl is a scalar. 87 */ 88 template <class NumberType> 89 struct CurlType<2, NumberType> 90 { 91 using type = Tensor<1, 1, NumberType>; 92 }; 93 94 /** 95 * A class whose specialization is used to define what type the curl of a 96 * vector valued function corresponds to. 97 * 98 * In 3d, the curl is a vector. 99 */ 100 template <class NumberType> 101 struct CurlType<3, NumberType> 102 { 103 using type = Tensor<1, 3, NumberType>; 104 }; 105 } // namespace internal 106 107 108 109 /** 110 * A namespace for "views" on a FEValues, FEFaceValues, or FESubfaceValues 111 * object. A view represents only a certain part of the whole: whereas the 112 * FEValues object represents <i>all</i> values, gradients, or second 113 * derivatives of all components of a vector-valued element, views restrict 114 * the attention to only a single component or a subset of components. You 115 * typically get objects of classes defined in this namespace by applying 116 * FEValuesExtractors objects to a FEValues, FEFaceValues or FESubfaceValues 117 * objects using the square bracket operator. 118 * 119 * There are classes that present views for single scalar components, vector 120 * components consisting of <code>dim</code> elements, and symmetric second 121 * order tensor components consisting of <code>(dim*dim + dim)/2</code> 122 * elements 123 * 124 * See the description of the 125 * @ref vector_valued 126 * module for examples how to use the features of this namespace. 127 * 128 * @ingroup feaccess vector_valued 129 */ 130 namespace FEValuesViews 131 { 132 /** 133 * A class representing a view to a single scalar component of a possibly 134 * vector-valued finite element. Views are discussed in the 135 * @ref vector_valued 136 * module. 137 * 138 * You get an object of this type if you apply a FEValuesExtractors::Scalar 139 * to an FEValues, FEFaceValues or FESubfaceValues object. 140 * 141 * @ingroup feaccess vector_valued 142 */ 143 template <int dim, int spacedim = dim> 144 class Scalar 145 { 146 public: 147 /** 148 * An alias for the data type of values of the view this class 149 * represents. Since we deal with a single components, the value type is a 150 * scalar double. 151 */ 152 using value_type = double; 153 154 /** 155 * An alias for the type of gradients of the view this class represents. 156 * Here, for a scalar component of the finite element, the gradient is a 157 * <code>Tensor@<1,dim@></code>. 158 */ 159 using gradient_type = dealii::Tensor<1, spacedim>; 160 161 /** 162 * An alias for the type of second derivatives of the view this class 163 * represents. Here, for a scalar component of the finite element, the 164 * Hessian is a <code>Tensor@<2,dim@></code>. 165 */ 166 using hessian_type = dealii::Tensor<2, spacedim>; 167 168 /** 169 * An alias for the type of third derivatives of the view this class 170 * represents. Here, for a scalar component of the finite element, the 171 * Third derivative is a <code>Tensor@<3,dim@></code>. 172 */ 173 using third_derivative_type = dealii::Tensor<3, spacedim>; 174 175 /** 176 * A struct that provides the output type for the product of the value 177 * and derivatives of basis functions of the Scalar view and any @p Number type. 178 */ 179 template <typename Number> 180 struct OutputType 181 { 182 /** 183 * An alias for the data type of the product of a @p Number and the 184 * values of the view the Scalar class. 185 */ 186 using value_type = 187 typename ProductType<Number, 188 typename Scalar<dim, spacedim>::value_type>::type; 189 190 /** 191 * An alias for the data type of the product of a @p Number and the 192 * gradients of the view the Scalar class. 193 */ 194 using gradient_type = typename ProductType< 195 Number, 196 typename Scalar<dim, spacedim>::gradient_type>::type; 197 198 /** 199 * An alias for the data type of the product of a @p Number and the 200 * laplacians of the view the Scalar class. 201 */ 202 using laplacian_type = 203 typename ProductType<Number, 204 typename Scalar<dim, spacedim>::value_type>::type; 205 206 /** 207 * An alias for the data type of the product of a @p Number and the 208 * hessians of the view the Scalar class. 209 */ 210 using hessian_type = typename ProductType< 211 Number, 212 typename Scalar<dim, spacedim>::hessian_type>::type; 213 214 /** 215 * An alias for the data type of the product of a @p Number and the 216 * third derivatives of the view the Scalar class. 217 */ 218 using third_derivative_type = typename ProductType< 219 Number, 220 typename Scalar<dim, spacedim>::third_derivative_type>::type; 221 }; 222 223 /** 224 * A structure where for each shape function we pre-compute a bunch of 225 * data that will make later accesses much cheaper. 226 */ 227 struct ShapeFunctionData 228 { 229 /** 230 * For each shape function, store whether the selected vector component 231 * may be nonzero. For primitive shape functions we know for sure 232 * whether a certain scalar component of a given shape function is 233 * nonzero, whereas for non-primitive shape functions this may not be 234 * entirely clear (e.g. for RT elements it depends on the shape of a 235 * cell). 236 */ 237 bool is_nonzero_shape_function_component; 238 239 /** 240 * For each shape function, store the row index within the shape_values, 241 * shape_gradients, and shape_hessians tables (the column index is the 242 * quadrature point index). If the shape function is primitive, then we 243 * can get this information from the shape_function_to_row_table of the 244 * FEValues object; otherwise, we have to work a bit harder to compute 245 * this information. 246 */ 247 unsigned int row_index; 248 }; 249 250 /** 251 * Default constructor. Creates an invalid object. 252 */ 253 Scalar(); 254 255 /** 256 * Constructor for an object that represents a single scalar component of 257 * a FEValuesBase object (or of one of the classes derived from 258 * FEValuesBase). 259 */ 260 Scalar(const FEValuesBase<dim, spacedim> &fe_values_base, 261 const unsigned int component); 262 263 /** 264 * Copy operator. This is not a lightweight object so we don't allow 265 * copying and generate a compile-time error if this function is called. 266 */ 267 Scalar & 268 operator=(const Scalar<dim, spacedim> &) = delete; 269 270 /** 271 * Return the value of the vector component selected by this view, for the 272 * shape function and quadrature point selected by the arguments. 273 * 274 * @param shape_function Number of the shape function to be evaluated. 275 * Note that this number runs from zero to dofs_per_cell, even in the case 276 * of an FEFaceValues or FESubfaceValues object. 277 * 278 * @param q_point Number of the quadrature point at which function is to 279 * be evaluated. 280 * 281 * @dealiiRequiresUpdateFlags{update_values} 282 */ 283 value_type 284 value(const unsigned int shape_function, const unsigned int q_point) const; 285 286 /** 287 * Return the gradient (a tensor of rank 1) of the vector component 288 * selected by this view, for the shape function and quadrature point 289 * selected by the arguments. 290 * 291 * @note The meaning of the arguments is as documented for the value() 292 * function. 293 * 294 * @dealiiRequiresUpdateFlags{update_gradients} 295 */ 296 gradient_type 297 gradient(const unsigned int shape_function, 298 const unsigned int q_point) const; 299 300 /** 301 * Return the Hessian (the tensor of rank 2 of all second derivatives) of 302 * the vector component selected by this view, for the shape function and 303 * quadrature point selected by the arguments. 304 * 305 * @note The meaning of the arguments is as documented for the value() 306 * function. 307 * 308 * @dealiiRequiresUpdateFlags{update_hessians} 309 */ 310 hessian_type 311 hessian(const unsigned int shape_function, 312 const unsigned int q_point) const; 313 314 /** 315 * Return the tensor of rank 3 of all third derivatives of the vector 316 * component selected by this view, for the shape function and quadrature 317 * point selected by the arguments. 318 * 319 * @note The meaning of the arguments is as documented for the value() 320 * function. 321 * 322 * @dealiiRequiresUpdateFlags{update_third_derivatives} 323 */ 324 third_derivative_type 325 third_derivative(const unsigned int shape_function, 326 const unsigned int q_point) const; 327 328 /** 329 * Return the values of the selected scalar component of the finite 330 * element function characterized by <tt>fe_function</tt> at the 331 * quadrature points of the cell, face or subface selected the last time 332 * the <tt>reinit</tt> function of the FEValues object was called. 333 * 334 * This function is the equivalent of the 335 * FEValuesBase::get_function_values function but it only works on the 336 * selected scalar component. 337 * 338 * The data type stored by the output vector must be what you get when you 339 * multiply the values of shape functions (i.e., @p value_type) times the 340 * type used to store the values of the unknowns $U_j$ of your finite 341 * element vector $U$ (represented by the @p fe_function argument). 342 * 343 * @dealiiRequiresUpdateFlags{update_values} 344 */ 345 template <class InputVector> 346 void 347 get_function_values( 348 const InputVector &fe_function, 349 std::vector<typename ProductType<value_type, 350 typename InputVector::value_type>::type> 351 &values) const; 352 353 /** 354 * Same as above, but using a vector of local degree-of-freedom values. 355 * 356 * The @p dof_values vector must have a length equal to number of DoFs on 357 * a cell, and each entry @p dof_values[i] is the value of the local DoF 358 * @p i. The fundamental prerequisite for the @p InputVector is that it must 359 * be possible to create an ArrayView from it; this is satisfied by the 360 * @p std::vector class. 361 * 362 * The DoF values typically would be obtained in the following way: 363 * @code 364 * Vector<double> local_dof_values(cell->get_fe().n_dofs_per_cell()); 365 * cell->get_dof_values(solution, local_dof_values); 366 * @endcode 367 * or, for a generic @p Number type, 368 * @code 369 * std::vector<Number> local_dof_values(cell->get_fe().n_dofs_per_cell()); 370 * cell->get_dof_values(solution, 371 * local_dof_values.begin(), 372 * local_dof_values.end()); 373 * @endcode 374 */ 375 template <class InputVector> 376 void 377 get_function_values_from_local_dof_values( 378 const InputVector &dof_values, 379 std::vector< 380 typename OutputType<typename InputVector::value_type>::value_type> 381 &values) const; 382 383 /** 384 * Return the gradients of the selected scalar component of the finite 385 * element function characterized by <tt>fe_function</tt> at the 386 * quadrature points of the cell, face or subface selected the last time 387 * the <tt>reinit</tt> function of the FEValues object was called. 388 * 389 * This function is the equivalent of the 390 * FEValuesBase::get_function_gradients function but it only works on the 391 * selected scalar component. 392 * 393 * The data type stored by the output vector must be what you get when you 394 * multiply the gradients of shape functions (i.e., @p gradient_type) 395 * times the type used to store the values of the unknowns $U_j$ of your 396 * finite element vector $U$ (represented by the @p fe_function argument). 397 * 398 * @dealiiRequiresUpdateFlags{update_gradients} 399 */ 400 template <class InputVector> 401 void 402 get_function_gradients( 403 const InputVector &fe_function, 404 std::vector<typename ProductType<gradient_type, 405 typename InputVector::value_type>::type> 406 &gradients) const; 407 408 /** 409 * @copydoc FEValuesViews::Scalar::get_function_values_from_local_dof_values() 410 */ 411 template <class InputVector> 412 void 413 get_function_gradients_from_local_dof_values( 414 const InputVector &dof_values, 415 std::vector< 416 typename OutputType<typename InputVector::value_type>::gradient_type> 417 &gradients) const; 418 419 /** 420 * Return the Hessians of the selected scalar component of the finite 421 * element function characterized by <tt>fe_function</tt> at the 422 * quadrature points of the cell, face or subface selected the last time 423 * the <tt>reinit</tt> function of the FEValues object was called. 424 * 425 * This function is the equivalent of the 426 * FEValuesBase::get_function_hessians function but it only works on the 427 * selected scalar component. 428 * 429 * The data type stored by the output vector must be what you get when you 430 * multiply the Hessians of shape functions (i.e., @p hessian_type) times 431 * the type used to store the values of the unknowns $U_j$ of your finite 432 * element vector $U$ (represented by the @p fe_function argument). 433 * 434 * @dealiiRequiresUpdateFlags{update_hessians} 435 */ 436 template <class InputVector> 437 void 438 get_function_hessians( 439 const InputVector &fe_function, 440 std::vector<typename ProductType<hessian_type, 441 typename InputVector::value_type>::type> 442 &hessians) const; 443 444 /** 445 * @copydoc FEValuesViews::Scalar::get_function_values_from_local_dof_values() 446 */ 447 template <class InputVector> 448 void 449 get_function_hessians_from_local_dof_values( 450 const InputVector &dof_values, 451 std::vector< 452 typename OutputType<typename InputVector::value_type>::hessian_type> 453 &hessians) const; 454 455 456 /** 457 * Return the Laplacians of the selected scalar component of the finite 458 * element function characterized by <tt>fe_function</tt> at the 459 * quadrature points of the cell, face or subface selected the last time 460 * the <tt>reinit</tt> function of the FEValues object was called. The 461 * Laplacians are the trace of the Hessians. 462 * 463 * This function is the equivalent of the 464 * FEValuesBase::get_function_laplacians function but it only works on the 465 * selected scalar component. 466 * 467 * The data type stored by the output vector must be what you get when you 468 * multiply the Laplacians of shape functions (i.e., @p value_type) times 469 * the type used to store the values of the unknowns $U_j$ of your finite 470 * element vector $U$ (represented by the @p fe_function argument). 471 * 472 * @dealiiRequiresUpdateFlags{update_hessians} 473 */ 474 template <class InputVector> 475 void 476 get_function_laplacians( 477 const InputVector &fe_function, 478 std::vector<typename ProductType<value_type, 479 typename InputVector::value_type>::type> 480 &laplacians) const; 481 482 /** 483 * @copydoc FEValuesViews::Scalar::get_function_values_from_local_dof_values() 484 */ 485 template <class InputVector> 486 void 487 get_function_laplacians_from_local_dof_values( 488 const InputVector &dof_values, 489 std::vector< 490 typename OutputType<typename InputVector::value_type>::laplacian_type> 491 &laplacians) const; 492 493 494 /** 495 * Return the third derivatives of the selected scalar component of the 496 * finite element function characterized by <tt>fe_function</tt> at the 497 * quadrature points of the cell, face or subface selected the last time 498 * the <tt>reinit</tt> function of the FEValues object was called. 499 * 500 * This function is the equivalent of the 501 * FEValuesBase::get_function_third_derivatives function but it only works 502 * on the selected scalar component. 503 * 504 * The data type stored by the output vector must be what you get when you 505 * multiply the third derivatives of shape functions (i.e., @p 506 * third_derivative_type) times the type used to store the values of the 507 * unknowns $U_j$ of your finite element vector $U$ (represented by the @p 508 * fe_function argument). 509 * 510 * @dealiiRequiresUpdateFlags{update_third_derivatives} 511 */ 512 template <class InputVector> 513 void 514 get_function_third_derivatives( 515 const InputVector &fe_function, 516 std::vector<typename ProductType<third_derivative_type, 517 typename InputVector::value_type>::type> 518 &third_derivatives) const; 519 520 /** 521 * @copydoc FEValuesViews::Scalar::get_function_values_from_local_dof_values() 522 */ 523 template <class InputVector> 524 void 525 get_function_third_derivatives_from_local_dof_values( 526 const InputVector & dof_values, 527 std::vector<typename OutputType<typename InputVector::value_type>:: 528 third_derivative_type> &third_derivatives) const; 529 530 531 private: 532 /** 533 * A pointer to the FEValuesBase object we operate on. 534 */ 535 const SmartPointer<const FEValuesBase<dim, spacedim>> fe_values; 536 537 /** 538 * The single scalar component this view represents of the FEValuesBase 539 * object. 540 */ 541 const unsigned int component; 542 543 /** 544 * Store the data about shape functions. 545 */ 546 std::vector<ShapeFunctionData> shape_function_data; 547 }; 548 549 550 551 /** 552 * A class representing a view to a set of <code>spacedim</code> components 553 * forming a vector part of a vector-valued finite element. Views are 554 * discussed in the 555 * @ref vector_valued 556 * module. 557 * 558 * Note that in the current context, a vector is meant in the sense physics 559 * uses it: it has <code>spacedim</code> components that behave in specific 560 * ways under coordinate system transformations. Examples include velocity 561 * or displacement fields. This is opposed to how mathematics uses the word 562 * "vector" (and how we use this word in other contexts in the library, for 563 * example in the Vector class), where it really stands for a collection of 564 * numbers. An example of this latter use of the word could be the set of 565 * concentrations of chemical species in a flame; however, these are really 566 * just a collection of scalar variables, since they do not change if the 567 * coordinate system is rotated, unlike the components of a velocity vector, 568 * and consequently, this class should not be used for this context. 569 * 570 * This class allows to query the value, gradient and divergence of 571 * (components of) shape functions and solutions representing vectors. The 572 * gradient of a vector $d_{k}, 0\le k<\text{dim}$ is defined as $S_{ij} = 573 * \frac{\partial d_{i}}{\partial x_j}, 0\le i,j<\text{dim}$. 574 * 575 * You get an object of this type if you apply a FEValuesExtractors::Vector 576 * to an FEValues, FEFaceValues or FESubfaceValues object. 577 * 578 * @ingroup feaccess vector_valued 579 */ 580 template <int dim, int spacedim = dim> 581 class Vector 582 { 583 public: 584 /** 585 * An alias for the data type of values of the view this class 586 * represents. Since we deal with a set of <code>dim</code> components, 587 * the value type is a Tensor<1,spacedim>. 588 */ 589 using value_type = dealii::Tensor<1, spacedim>; 590 591 /** 592 * An alias for the type of gradients of the view this class represents. 593 * Here, for a set of <code>dim</code> components of the finite element, 594 * the gradient is a <code>Tensor@<2,spacedim@></code>. 595 * 596 * See the general documentation of this class for how exactly the 597 * gradient of a vector is defined. 598 */ 599 using gradient_type = dealii::Tensor<2, spacedim>; 600 601 /** 602 * An alias for the type of symmetrized gradients of the view this class 603 * represents. Here, for a set of <code>dim</code> components of the 604 * finite element, the symmetrized gradient is a 605 * <code>SymmetricTensor@<2,spacedim@></code>. 606 * 607 * The symmetric gradient of a vector field $\mathbf v$ is defined as 608 * $\varepsilon(\mathbf v)=\frac 12 (\nabla \mathbf v + \nabla \mathbf 609 * v^T)$. 610 */ 611 using symmetric_gradient_type = dealii::SymmetricTensor<2, spacedim>; 612 613 /** 614 * An alias for the type of the divergence of the view this class 615 * represents. Here, for a set of <code>dim</code> components of the 616 * finite element, the divergence of course is a scalar. 617 */ 618 using divergence_type = double; 619 620 /** 621 * An alias for the type of the curl of the view this class represents. 622 * Here, for a set of <code>spacedim=2</code> components of the finite 623 * element, the curl is a <code>Tensor@<1, 1@></code>. For 624 * <code>spacedim=3</code> it is a <code>Tensor@<1, dim@></code>. 625 */ 626 using curl_type = typename dealii::internal::CurlType<spacedim>::type; 627 628 /** 629 * An alias for the type of second derivatives of the view this class 630 * represents. Here, for a set of <code>dim</code> components of the 631 * finite element, the Hessian is a <code>Tensor@<3,dim@></code>. 632 */ 633 using hessian_type = dealii::Tensor<3, spacedim>; 634 635 /** 636 * An alias for the type of third derivatives of the view this class 637 * represents. Here, for a set of <code>dim</code> components of the 638 * finite element, the third derivative is a <code>Tensor@<4,dim@></code>. 639 */ 640 using third_derivative_type = dealii::Tensor<4, spacedim>; 641 642 /** 643 * A struct that provides the output type for the product of the value 644 * and derivatives of basis functions of the Vector view and any @p Number type. 645 */ 646 template <typename Number> 647 struct OutputType 648 { 649 /** 650 * An alias for the data type of the product of a @p Number and the 651 * values of the view the Vector class. 652 */ 653 using value_type = 654 typename ProductType<Number, 655 typename Vector<dim, spacedim>::value_type>::type; 656 657 /** 658 * An alias for the data type of the product of a @p Number and the 659 * gradients of the view the Vector class. 660 */ 661 using gradient_type = typename ProductType< 662 Number, 663 typename Vector<dim, spacedim>::gradient_type>::type; 664 665 /** 666 * An alias for the data type of the product of a @p Number and the 667 * symmetric gradients of the view the Vector class. 668 */ 669 using symmetric_gradient_type = typename ProductType< 670 Number, 671 typename Vector<dim, spacedim>::symmetric_gradient_type>::type; 672 673 /** 674 * An alias for the data type of the product of a @p Number and the 675 * divergences of the view the Vector class. 676 */ 677 using divergence_type = typename ProductType< 678 Number, 679 typename Vector<dim, spacedim>::divergence_type>::type; 680 681 /** 682 * An alias for the data type of the product of a @p Number and the 683 * laplacians of the view the Vector class. 684 */ 685 using laplacian_type = 686 typename ProductType<Number, 687 typename Vector<dim, spacedim>::value_type>::type; 688 689 /** 690 * An alias for the data type of the product of a @p Number and the 691 * curls of the view the Vector class. 692 */ 693 using curl_type = 694 typename ProductType<Number, 695 typename Vector<dim, spacedim>::curl_type>::type; 696 697 /** 698 * An alias for the data type of the product of a @p Number and the 699 * hessians of the view the Vector class. 700 */ 701 using hessian_type = typename ProductType< 702 Number, 703 typename Vector<dim, spacedim>::hessian_type>::type; 704 705 /** 706 * An alias for the data type of the product of a @p Number and the 707 * third derivatives of the view the Vector class. 708 */ 709 using third_derivative_type = typename ProductType< 710 Number, 711 typename Vector<dim, spacedim>::third_derivative_type>::type; 712 }; 713 714 /** 715 * A structure where for each shape function we pre-compute a bunch of 716 * data that will make later accesses much cheaper. 717 */ 718 struct ShapeFunctionData 719 { 720 /** 721 * For each pair (shape function,component within vector), store whether 722 * the selected vector component may be nonzero. For primitive shape 723 * functions we know for sure whether a certain scalar component of a 724 * given shape function is nonzero, whereas for non-primitive shape 725 * functions this may not be entirely clear (e.g. for RT elements it 726 * depends on the shape of a cell). 727 */ 728 bool is_nonzero_shape_function_component[spacedim]; 729 730 /** 731 * For each pair (shape function, component within vector), store the 732 * row index within the shape_values, shape_gradients, and 733 * shape_hessians tables (the column index is the quadrature point 734 * index). If the shape function is primitive, then we can get this 735 * information from the shape_function_to_row_table of the FEValues 736 * object; otherwise, we have to work a bit harder to compute this 737 * information. 738 */ 739 unsigned int row_index[spacedim]; 740 741 /** 742 * For each shape function say the following: if only a single entry in 743 * is_nonzero_shape_function_component for this shape function is 744 * nonzero, then store the corresponding value of row_index and 745 * single_nonzero_component_index represents the index between 0 and dim 746 * for which it is attained. If multiple components are nonzero, then 747 * store -1. If no components are nonzero then store -2. 748 */ 749 int single_nonzero_component; 750 unsigned int single_nonzero_component_index; 751 }; 752 753 /** 754 * Default constructor. Creates an invalid object. 755 */ 756 Vector(); 757 758 /** 759 * Constructor for an object that represents dim components of a 760 * FEValuesBase object (or of one of the classes derived from 761 * FEValuesBase), representing a vector-valued variable. 762 * 763 * The second argument denotes the index of the first component of the 764 * selected vector. 765 */ 766 Vector(const FEValuesBase<dim, spacedim> &fe_values_base, 767 const unsigned int first_vector_component); 768 769 /** 770 * Copy operator. This is not a lightweight object so we don't allow 771 * copying and generate a compile-time error if this function is called. 772 */ 773 Vector & 774 operator=(const Vector<dim, spacedim> &) = delete; 775 776 /** 777 * Return the value of the vector components selected by this view, for 778 * the shape function and quadrature point selected by the arguments. 779 * Here, since the view represents a vector-valued part of the FEValues 780 * object with <code>dim</code> components, the return type is a tensor of 781 * rank 1 with <code>dim</code> components. 782 * 783 * @param shape_function Number of the shape function to be evaluated. 784 * Note that this number runs from zero to dofs_per_cell, even in the case 785 * of an FEFaceValues or FESubfaceValues object. 786 * 787 * @param q_point Number of the quadrature point at which function is to 788 * be evaluated. 789 * 790 * @dealiiRequiresUpdateFlags{update_values} 791 */ 792 value_type 793 value(const unsigned int shape_function, const unsigned int q_point) const; 794 795 /** 796 * Return the gradient (a tensor of rank 2) of the vector component 797 * selected by this view, for the shape function and quadrature point 798 * selected by the arguments. 799 * 800 * See the general documentation of this class for how exactly the 801 * gradient of a vector is defined. 802 * 803 * @note The meaning of the arguments is as documented for the value() 804 * function. 805 * 806 * @dealiiRequiresUpdateFlags{update_gradients} 807 */ 808 gradient_type 809 gradient(const unsigned int shape_function, 810 const unsigned int q_point) const; 811 812 /** 813 * Return the symmetric gradient (a symmetric tensor of rank 2) of the 814 * vector component selected by this view, for the shape function and 815 * quadrature point selected by the arguments. 816 * 817 * The symmetric gradient is defined as $\frac 12 [(\nabla \phi_i(x_q)) + 818 * (\nabla \phi_i(x_q))^T]$, where $\phi_i$ represents the 819 * <code>dim</code> components selected from the FEValuesBase object, and 820 * $x_q$ is the location of the $q$-th quadrature point. 821 * 822 * @note The meaning of the arguments is as documented for the value() 823 * function. 824 * 825 * @dealiiRequiresUpdateFlags{update_gradients} 826 */ 827 symmetric_gradient_type 828 symmetric_gradient(const unsigned int shape_function, 829 const unsigned int q_point) const; 830 831 /** 832 * Return the scalar divergence of the vector components selected by this 833 * view, for the shape function and quadrature point selected by the 834 * arguments. 835 * 836 * @note The meaning of the arguments is as documented for the value() 837 * function. 838 * 839 * @dealiiRequiresUpdateFlags{update_gradients} 840 */ 841 divergence_type 842 divergence(const unsigned int shape_function, 843 const unsigned int q_point) const; 844 845 /** 846 * Return the vector curl of the vector components selected by this view, 847 * for the shape function and quadrature point selected by the arguments. 848 * For 1d this function does not make any sense. Thus it is not 849 * implemented for <code>spacedim=1</code>. In 2d the curl is defined as 850 * @f{equation*}{ 851 * \operatorname{curl}(u) \dealcoloneq \frac{du_2}{dx} -\frac{du_1}{dy}, 852 * @f} 853 * whereas in 3d it is given by 854 * @f{equation*}{ 855 * \operatorname{curl}(u) \dealcoloneq \left( \begin{array}{c} 856 * \frac{du_3}{dy}-\frac{du_2}{dz}\\ \frac{du_1}{dz}-\frac{du_3}{dx}\\ 857 * \frac{du_2}{dx}-\frac{du_1}{dy} \end{array} \right). 858 * @f} 859 * 860 * @note The meaning of the arguments is as documented for the value() 861 * function. 862 * 863 * @dealiiRequiresUpdateFlags{update_gradients} 864 */ 865 curl_type 866 curl(const unsigned int shape_function, const unsigned int q_point) const; 867 868 /** 869 * Return the Hessian (the tensor of rank 2 of all second derivatives) of 870 * the vector components selected by this view, for the shape function and 871 * quadrature point selected by the arguments. 872 * 873 * @note The meaning of the arguments is as documented for the value() 874 * function. 875 * 876 * @dealiiRequiresUpdateFlags{update_hessians} 877 */ 878 hessian_type 879 hessian(const unsigned int shape_function, 880 const unsigned int q_point) const; 881 882 /** 883 * Return the tensor of rank 3 of all third derivatives of the vector 884 * components selected by this view, for the shape function and quadrature 885 * point selected by the arguments. 886 * 887 * @note The meaning of the arguments is as documented for the value() 888 * function. 889 * 890 * @dealiiRequiresUpdateFlags{update_3rd_derivatives} 891 */ 892 third_derivative_type 893 third_derivative(const unsigned int shape_function, 894 const unsigned int q_point) const; 895 896 /** 897 * Return the values of the selected vector components of the finite 898 * element function characterized by <tt>fe_function</tt> at the 899 * quadrature points of the cell, face or subface selected the last time 900 * the <tt>reinit</tt> function of the FEValues object was called. 901 * 902 * This function is the equivalent of the 903 * FEValuesBase::get_function_values function but it only works on the 904 * selected vector components. 905 * 906 * The data type stored by the output vector must be what you get when you 907 * multiply the values of shape functions (i.e., @p value_type) times the 908 * type used to store the values of the unknowns $U_j$ of your finite 909 * element vector $U$ (represented by the @p fe_function argument). 910 * 911 * @dealiiRequiresUpdateFlags{update_values} 912 */ 913 template <class InputVector> 914 void 915 get_function_values( 916 const InputVector &fe_function, 917 std::vector<typename ProductType<value_type, 918 typename InputVector::value_type>::type> 919 &values) const; 920 921 /** 922 * Same as above, but using a vector of local degree-of-freedom values. 923 * 924 * The @p dof_values vector must have a length equal to number of DoFs on 925 * a cell, and each entry @p dof_values[i] is the value of the local DoF 926 * @p i. The fundamental prerequisite for the @p InputVector is that it must 927 * be possible to create an ArrayView from it; this is satisfied by the 928 * @p std::vector class. 929 * 930 * The DoF values typically would be obtained in the following way: 931 * @code 932 * Vector<double> local_dof_values(cell->get_fe().n_dofs_per_cell()); 933 * cell->get_dof_values(solution, local_dof_values); 934 * @endcode 935 * or, for a generic @p Number type, 936 * @code 937 * std::vector<Number> local_dof_values(cell->get_fe().n_dofs_per_cell()); 938 * cell->get_dof_values(solution, 939 * local_dof_values.begin(), 940 * local_dof_values.end()); 941 * @endcode 942 */ 943 template <class InputVector> 944 void 945 get_function_values_from_local_dof_values( 946 const InputVector &dof_values, 947 std::vector< 948 typename OutputType<typename InputVector::value_type>::value_type> 949 &values) const; 950 951 /** 952 * Return the gradients of the selected vector components of the finite 953 * element function characterized by <tt>fe_function</tt> at the 954 * quadrature points of the cell, face or subface selected the last time 955 * the <tt>reinit</tt> function of the FEValues object was called. 956 * 957 * This function is the equivalent of the 958 * FEValuesBase::get_function_gradients function but it only works on the 959 * selected vector components. 960 * 961 * The data type stored by the output vector must be what you get when you 962 * multiply the gradients of shape functions (i.e., @p gradient_type) 963 * times the type used to store the values of the unknowns $U_j$ of your 964 * finite element vector $U$ (represented by the @p fe_function argument). 965 * 966 * @dealiiRequiresUpdateFlags{update_gradients} 967 */ 968 template <class InputVector> 969 void 970 get_function_gradients( 971 const InputVector &fe_function, 972 std::vector<typename ProductType<gradient_type, 973 typename InputVector::value_type>::type> 974 &gradients) const; 975 976 /** 977 * @copydoc FEValuesViews::Vector::get_function_values_from_local_dof_values() 978 */ 979 template <class InputVector> 980 void 981 get_function_gradients_from_local_dof_values( 982 const InputVector &dof_values, 983 std::vector< 984 typename OutputType<typename InputVector::value_type>::gradient_type> 985 &gradients) const; 986 987 /** 988 * Return the symmetrized gradients of the selected vector components of 989 * the finite element function characterized by <tt>fe_function</tt> at 990 * the quadrature points of the cell, face or subface selected the last 991 * time the <tt>reinit</tt> function of the FEValues object was called. 992 * 993 * The symmetric gradient of a vector field $\mathbf v$ is defined as 994 * $\varepsilon(\mathbf v)=\frac 12 (\nabla \mathbf v + \nabla \mathbf 995 * v^T)$. 996 * 997 * @note There is no equivalent function such as 998 * FEValuesBase::get_function_symmetric_gradients in the FEValues classes 999 * but the information can be obtained from 1000 * FEValuesBase::get_function_gradients, of course. 1001 * 1002 * The data type stored by the output vector must be what you get when you 1003 * multiply the symmetric gradients of shape functions (i.e., @p 1004 * symmetric_gradient_type) times the type used to store the values of the 1005 * unknowns $U_j$ of your finite element vector $U$ (represented by the @p 1006 * fe_function argument). 1007 * 1008 * @dealiiRequiresUpdateFlags{update_gradients} 1009 */ 1010 template <class InputVector> 1011 void 1012 get_function_symmetric_gradients( 1013 const InputVector &fe_function, 1014 std::vector<typename ProductType<symmetric_gradient_type, 1015 typename InputVector::value_type>::type> 1016 &symmetric_gradients) const; 1017 1018 /** 1019 * @copydoc FEValuesViews::Vector::get_function_values_from_local_dof_values() 1020 */ 1021 template <class InputVector> 1022 void 1023 get_function_symmetric_gradients_from_local_dof_values( 1024 const InputVector & dof_values, 1025 std::vector<typename OutputType<typename InputVector::value_type>:: 1026 symmetric_gradient_type> &symmetric_gradients) const; 1027 1028 /** 1029 * Return the divergence of the selected vector components of the finite 1030 * element function characterized by <tt>fe_function</tt> at the 1031 * quadrature points of the cell, face or subface selected the last time 1032 * the <tt>reinit</tt> function of the FEValues object was called. 1033 * 1034 * There is no equivalent function such as 1035 * FEValuesBase::get_function_divergences in the FEValues classes but the 1036 * information can be obtained from FEValuesBase::get_function_gradients, 1037 * of course. 1038 * 1039 * The data type stored by the output vector must be what you get when you 1040 * multiply the divergences of shape functions (i.e., @p divergence_type) 1041 * times the type used to store the values of the unknowns $U_j$ of your 1042 * finite element vector $U$ (represented by the @p fe_function argument). 1043 * 1044 * @dealiiRequiresUpdateFlags{update_gradients} 1045 */ 1046 template <class InputVector> 1047 void 1048 get_function_divergences( 1049 const InputVector &fe_function, 1050 std::vector<typename ProductType<divergence_type, 1051 typename InputVector::value_type>::type> 1052 &divergences) const; 1053 1054 /** 1055 * @copydoc FEValuesViews::Vector::get_function_values_from_local_dof_values() 1056 */ 1057 template <class InputVector> 1058 void 1059 get_function_divergences_from_local_dof_values( 1060 const InputVector &dof_values, 1061 std::vector< 1062 typename OutputType<typename InputVector::value_type>::divergence_type> 1063 &divergences) const; 1064 1065 /** 1066 * Return the curl of the selected vector components of the finite element 1067 * function characterized by <tt>fe_function</tt> at the quadrature points 1068 * of the cell, face or subface selected the last time the <tt>reinit</tt> 1069 * function of the FEValues object was called. 1070 * 1071 * There is no equivalent function such as 1072 * FEValuesBase::get_function_curls in the FEValues classes but the 1073 * information can be obtained from FEValuesBase::get_function_gradients, 1074 * of course. 1075 * 1076 * The data type stored by the output vector must be what you get when you 1077 * multiply the curls of shape functions (i.e., @p curl_type) times the 1078 * type used to store the values of the unknowns $U_j$ of your finite 1079 * element vector $U$ (represented by the @p fe_function argument). 1080 * 1081 * @dealiiRequiresUpdateFlags{update_gradients} 1082 */ 1083 template <class InputVector> 1084 void 1085 get_function_curls( 1086 const InputVector &fe_function, 1087 std::vector< 1088 typename ProductType<curl_type, typename InputVector::value_type>::type> 1089 &curls) const; 1090 1091 /** 1092 * @copydoc FEValuesViews::Vector::get_function_values_from_local_dof_values() 1093 */ 1094 template <class InputVector> 1095 void 1096 get_function_curls_from_local_dof_values( 1097 const InputVector &dof_values, 1098 std::vector< 1099 typename OutputType<typename InputVector::value_type>::curl_type> 1100 &curls) const; 1101 1102 /** 1103 * Return the Hessians of the selected vector components of the finite 1104 * element function characterized by <tt>fe_function</tt> at the 1105 * quadrature points of the cell, face or subface selected the last time 1106 * the <tt>reinit</tt> function of the FEValues object was called. 1107 * 1108 * This function is the equivalent of the 1109 * FEValuesBase::get_function_hessians function but it only works on the 1110 * selected vector components. 1111 * 1112 * The data type stored by the output vector must be what you get when you 1113 * multiply the Hessians of shape functions (i.e., @p hessian_type) times 1114 * the type used to store the values of the unknowns $U_j$ of your finite 1115 * element vector $U$ (represented by the @p fe_function argument). 1116 * 1117 * @dealiiRequiresUpdateFlags{update_hessians} 1118 */ 1119 template <class InputVector> 1120 void 1121 get_function_hessians( 1122 const InputVector &fe_function, 1123 std::vector<typename ProductType<hessian_type, 1124 typename InputVector::value_type>::type> 1125 &hessians) const; 1126 1127 /** 1128 * @copydoc FEValuesViews::Vector::get_function_values_from_local_dof_values() 1129 */ 1130 template <class InputVector> 1131 void 1132 get_function_hessians_from_local_dof_values( 1133 const InputVector &dof_values, 1134 std::vector< 1135 typename OutputType<typename InputVector::value_type>::hessian_type> 1136 &hessians) const; 1137 1138 /** 1139 * Return the Laplacians of the selected vector components of the finite 1140 * element function characterized by <tt>fe_function</tt> at the 1141 * quadrature points of the cell, face or subface selected the last time 1142 * the <tt>reinit</tt> function of the FEValues object was called. The 1143 * Laplacians are the trace of the Hessians. 1144 * 1145 * This function is the equivalent of the 1146 * FEValuesBase::get_function_laplacians function but it only works on the 1147 * selected vector components. 1148 * 1149 * The data type stored by the output vector must be what you get when you 1150 * multiply the Laplacians of shape functions (i.e., @p laplacian_type) 1151 * times the type used to store the values of the unknowns $U_j$ of your 1152 * finite element vector $U$ (represented by the @p fe_function argument). 1153 * 1154 * @dealiiRequiresUpdateFlags{update_hessians} 1155 */ 1156 template <class InputVector> 1157 void 1158 get_function_laplacians( 1159 const InputVector &fe_function, 1160 std::vector<typename ProductType<value_type, 1161 typename InputVector::value_type>::type> 1162 &laplacians) const; 1163 1164 /** 1165 * @copydoc FEValuesViews::Vector::get_function_values_from_local_dof_values() 1166 */ 1167 template <class InputVector> 1168 void 1169 get_function_laplacians_from_local_dof_values( 1170 const InputVector &dof_values, 1171 std::vector< 1172 typename OutputType<typename InputVector::value_type>::laplacian_type> 1173 &laplacians) const; 1174 1175 /** 1176 * Return the third derivatives of the selected scalar component of the 1177 * finite element function characterized by <tt>fe_function</tt> at the 1178 * quadrature points of the cell, face or subface selected the last time 1179 * the <tt>reinit</tt> function of the FEValues object was called. 1180 * 1181 * This function is the equivalent of the 1182 * FEValuesBase::get_function_third_derivatives function but it only works 1183 * on the selected scalar component. 1184 * 1185 * The data type stored by the output vector must be what you get when you 1186 * multiply the third derivatives of shape functions (i.e., @p 1187 * third_derivative_type) times the type used to store the values of the 1188 * unknowns $U_j$ of your finite element vector $U$ (represented by the @p 1189 * fe_function argument). 1190 * 1191 * @dealiiRequiresUpdateFlags{update_third_derivatives} 1192 */ 1193 template <class InputVector> 1194 void 1195 get_function_third_derivatives( 1196 const InputVector &fe_function, 1197 std::vector<typename ProductType<third_derivative_type, 1198 typename InputVector::value_type>::type> 1199 &third_derivatives) const; 1200 1201 /** 1202 * @copydoc FEValuesViews::Vector::get_function_values_from_local_dof_values() 1203 */ 1204 template <class InputVector> 1205 void 1206 get_function_third_derivatives_from_local_dof_values( 1207 const InputVector & dof_values, 1208 std::vector<typename OutputType<typename InputVector::value_type>:: 1209 third_derivative_type> &third_derivatives) const; 1210 1211 private: 1212 /** 1213 * A pointer to the FEValuesBase object we operate on. 1214 */ 1215 const SmartPointer<const FEValuesBase<dim, spacedim>> fe_values; 1216 1217 /** 1218 * The first component of the vector this view represents of the 1219 * FEValuesBase object. 1220 */ 1221 const unsigned int first_vector_component; 1222 1223 /** 1224 * Store the data about shape functions. 1225 */ 1226 std::vector<ShapeFunctionData> shape_function_data; 1227 }; 1228 1229 1230 template <int rank, int dim, int spacedim = dim> 1231 class SymmetricTensor; 1232 1233 /** 1234 * A class representing a view to a set of <code>(dim*dim + dim)/2</code> 1235 * components forming a symmetric second-order tensor from a vector-valued 1236 * finite element. Views are discussed in the 1237 * @ref vector_valued 1238 * module. 1239 * 1240 * This class allows to query the value and divergence of (components of) 1241 * shape functions and solutions representing symmetric tensors. The 1242 * divergence of a symmetric tensor $S_{ij}, 0\le i,j<\text{dim}$ is defined 1243 * as $d_i = \sum_j \frac{\partial S_{ij}}{\partial x_j}, 0\le 1244 * i<\text{dim}$, which due to the symmetry of the tensor is also $d_i = 1245 * \sum_j \frac{\partial S_{ji}}{\partial x_j}$. In other words, it due to 1246 * the symmetry of $S$ it does not matter whether we apply the nabla 1247 * operator by row or by column to get the divergence. 1248 * 1249 * You get an object of this type if you apply a 1250 * FEValuesExtractors::SymmetricTensor to an FEValues, FEFaceValues or 1251 * FESubfaceValues object. 1252 * 1253 * @ingroup feaccess vector_valued 1254 */ 1255 template <int dim, int spacedim> 1256 class SymmetricTensor<2, dim, spacedim> 1257 { 1258 public: 1259 /** 1260 * An alias for the data type of values of the view this class 1261 * represents. Since we deal with a set of <code>(dim*dim + dim)/2</code> 1262 * components (i.e. the unique components of a symmetric second-order 1263 * tensor), the value type is a SymmetricTensor<2,spacedim>. 1264 */ 1265 using value_type = dealii::SymmetricTensor<2, spacedim>; 1266 1267 /** 1268 * An alias for the type of the divergence of the view this class 1269 * represents. Here, for a set of <code>(dim*dim + dim)/2</code> unique 1270 * components of the finite element representing a symmetric second-order 1271 * tensor, the divergence of course is a * <code>Tensor@<1,dim@></code>. 1272 * 1273 * See the general discussion of this class for a definition of the 1274 * divergence. 1275 */ 1276 using divergence_type = dealii::Tensor<1, spacedim>; 1277 1278 /** 1279 * A struct that provides the output type for the product of the value 1280 * and derivatives of basis functions of the SymmetricTensor view and any @p Number type. 1281 */ 1282 template <typename Number> 1283 struct OutputType 1284 { 1285 /** 1286 * An alias for the data type of the product of a @p Number and the 1287 * values of the view the SymmetricTensor class. 1288 */ 1289 using value_type = typename ProductType< 1290 Number, 1291 typename SymmetricTensor<2, dim, spacedim>::value_type>::type; 1292 1293 /** 1294 * An alias for the data type of the product of a @p Number and the 1295 * divergences of the view the SymmetricTensor class. 1296 */ 1297 using divergence_type = typename ProductType< 1298 Number, 1299 typename SymmetricTensor<2, dim, spacedim>::divergence_type>::type; 1300 }; 1301 1302 /** 1303 * A structure where for each shape function we pre-compute a bunch of 1304 * data that will make later accesses much cheaper. 1305 */ 1306 struct ShapeFunctionData 1307 { 1308 /** 1309 * For each pair (shape function,component within vector), store whether 1310 * the selected vector component may be nonzero. For primitive shape 1311 * functions we know for sure whether a certain scalar component of a 1312 * given shape function is nonzero, whereas for non-primitive shape 1313 * functions this may not be entirely clear (e.g. for RT elements it 1314 * depends on the shape of a cell). 1315 */ 1316 bool is_nonzero_shape_function_component 1317 [value_type::n_independent_components]; 1318 1319 /** 1320 * For each pair (shape function, component within vector), store the 1321 * row index within the shape_values, shape_gradients, and 1322 * shape_hessians tables (the column index is the quadrature point 1323 * index). If the shape function is primitive, then we can get this 1324 * information from the shape_function_to_row_table of the FEValues 1325 * object; otherwise, we have to work a bit harder to compute this 1326 * information. 1327 */ 1328 unsigned int row_index[value_type::n_independent_components]; 1329 1330 /** 1331 * For each shape function say the following: if only a single entry in 1332 * is_nonzero_shape_function_component for this shape function is 1333 * nonzero, then store the corresponding value of row_index and 1334 * single_nonzero_component_index represents the index between 0 and 1335 * (dim^2 + dim)/2 for which it is attained. If multiple components are 1336 * nonzero, then store -1. If no components are nonzero then store -2. 1337 */ 1338 int single_nonzero_component; 1339 1340 /** 1341 * Index of the @p single_nonzero_component . 1342 */ 1343 unsigned int single_nonzero_component_index; 1344 }; 1345 1346 /** 1347 * Default constructor. Creates an invalid object. 1348 */ 1349 SymmetricTensor(); 1350 1351 /** 1352 * Constructor for an object that represents <code>(dim*dim + 1353 * dim)/2</code> components of a FEValuesBase object (or of one of the 1354 * classes derived from FEValuesBase), representing the unique components 1355 * comprising a symmetric second- order tensor valued variable. 1356 * 1357 * The second argument denotes the index of the first component of the 1358 * selected symmetric second order tensor. 1359 */ 1360 SymmetricTensor(const FEValuesBase<dim, spacedim> &fe_values_base, 1361 const unsigned int first_tensor_component); 1362 1363 /** 1364 * Copy operator. This is not a lightweight object so we don't allow 1365 * copying and generate a compile-time error if this function is called. 1366 */ 1367 SymmetricTensor & 1368 operator=(const SymmetricTensor<2, dim, spacedim> &) = delete; 1369 1370 /** 1371 * Return the value of the vector components selected by this view, for 1372 * the shape function and quadrature point selected by the arguments. 1373 * Here, since the view represents a vector-valued part of the FEValues 1374 * object with <code>(dim*dim + dim)/2</code> components (the unique 1375 * components of a symmetric second-order tensor), the return type is a 1376 * symmetric tensor of rank 2. 1377 * 1378 * @param shape_function Number of the shape function to be evaluated. 1379 * Note that this number runs from zero to dofs_per_cell, even in the case 1380 * of an FEFaceValues or FESubfaceValues object. 1381 * 1382 * @param q_point Number of the quadrature point at which function is to 1383 * be evaluated. 1384 * 1385 * @dealiiRequiresUpdateFlags{update_values} 1386 */ 1387 value_type 1388 value(const unsigned int shape_function, const unsigned int q_point) const; 1389 1390 /** 1391 * Return the vector divergence of the vector components selected by this 1392 * view, for the shape function and quadrature point selected by the 1393 * arguments. 1394 * 1395 * See the general discussion of this class for a definition of the 1396 * divergence. 1397 * 1398 * @note The meaning of the arguments is as documented for the value() 1399 * function. 1400 * 1401 * @dealiiRequiresUpdateFlags{update_gradients} 1402 */ 1403 divergence_type 1404 divergence(const unsigned int shape_function, 1405 const unsigned int q_point) const; 1406 1407 /** 1408 * Return the values of the selected vector components of the finite 1409 * element function characterized by <tt>fe_function</tt> at the 1410 * quadrature points of the cell, face or subface selected the last time 1411 * the <tt>reinit</tt> function of the FEValues object was called. 1412 * 1413 * This function is the equivalent of the 1414 * FEValuesBase::get_function_values function but it only works on the 1415 * selected vector components. 1416 * 1417 * The data type stored by the output vector must be what you get when you 1418 * multiply the values of shape functions (i.e., @p value_type) times the 1419 * type used to store the values of the unknowns $U_j$ of your finite 1420 * element vector $U$ (represented by the @p fe_function argument). 1421 * 1422 * @dealiiRequiresUpdateFlags{update_values} 1423 */ 1424 template <class InputVector> 1425 void 1426 get_function_values( 1427 const InputVector &fe_function, 1428 std::vector<typename ProductType<value_type, 1429 typename InputVector::value_type>::type> 1430 &values) const; 1431 1432 /** 1433 * Same as above, but using a vector of local degree-of-freedom values. 1434 * 1435 * The @p dof_values vector must have a length equal to number of DoFs on 1436 * a cell, and each entry @p dof_values[i] is the value of the local DoF 1437 * @p i. The fundamental prerequisite for the @p InputVector is that it must 1438 * be possible to create an ArrayView from it; this is satisfied by the 1439 * @p std::vector class. 1440 * 1441 * The DoF values typically would be obtained in the following way: 1442 * @code 1443 * Vector<double> local_dof_values(cell->get_fe().n_dofs_per_cell()); 1444 * cell->get_dof_values(solution, local_dof_values); 1445 * @endcode 1446 * or, for a generic @p Number type, 1447 * @code 1448 * std::vector<Number> local_dof_values(cell->get_fe().n_dofs_per_cell()); 1449 * cell->get_dof_values(solution, 1450 * local_dof_values.begin(), 1451 * local_dof_values.end()); 1452 * @endcode 1453 */ 1454 template <class InputVector> 1455 void 1456 get_function_values_from_local_dof_values( 1457 const InputVector &dof_values, 1458 std::vector< 1459 typename OutputType<typename InputVector::value_type>::value_type> 1460 &values) const; 1461 1462 /** 1463 * Return the divergence of the selected vector components of the finite 1464 * element function characterized by <tt>fe_function</tt> at the 1465 * quadrature points of the cell, face or subface selected the last time 1466 * the <tt>reinit</tt> function of the FEValues object was called. 1467 * 1468 * There is no equivalent function such as 1469 * FEValuesBase::get_function_divergences in the FEValues classes but the 1470 * information can be obtained from FEValuesBase::get_function_gradients, 1471 * of course. 1472 * 1473 * See the general discussion of this class for a definition of the 1474 * divergence. 1475 * 1476 * The data type stored by the output vector must be what you get when you 1477 * multiply the divergences of shape functions (i.e., @p divergence_type) 1478 * times the type used to store the values of the unknowns $U_j$ of your 1479 * finite element vector $U$ (represented by the @p fe_function argument). 1480 * 1481 * @dealiiRequiresUpdateFlags{update_gradients} 1482 */ 1483 template <class InputVector> 1484 void 1485 get_function_divergences( 1486 const InputVector &fe_function, 1487 std::vector<typename ProductType<divergence_type, 1488 typename InputVector::value_type>::type> 1489 &divergences) const; 1490 1491 /** 1492 * @copydoc FEValuesViews::SymmetricTensor<2,dim,spacedim>::get_function_values_from_local_dof_values() 1493 */ 1494 template <class InputVector> 1495 void 1496 get_function_divergences_from_local_dof_values( 1497 const InputVector &dof_values, 1498 std::vector< 1499 typename OutputType<typename InputVector::value_type>::divergence_type> 1500 &divergences) const; 1501 1502 private: 1503 /** 1504 * A pointer to the FEValuesBase object we operate on. 1505 */ 1506 const SmartPointer<const FEValuesBase<dim, spacedim>> fe_values; 1507 1508 /** 1509 * The first component of the vector this view represents of the 1510 * FEValuesBase object. 1511 */ 1512 const unsigned int first_tensor_component; 1513 1514 /** 1515 * Store the data about shape functions. 1516 */ 1517 std::vector<ShapeFunctionData> shape_function_data; 1518 }; 1519 1520 1521 template <int rank, int dim, int spacedim = dim> 1522 class Tensor; 1523 1524 /** 1525 * A class representing a view to a set of <code>dim*dim</code> components 1526 * forming a second-order tensor from a vector-valued finite element. Views 1527 * are discussed in the 1528 * @ref vector_valued 1529 * module. 1530 * 1531 * This class allows to query the value, gradient and divergence of 1532 * (components of) shape functions and solutions representing tensors. The 1533 * divergence of a tensor $T_{ij},\, 0\le i,j<\text{dim}$ is defined as $d_i = 1534 * \sum_j \frac{\partial T_{ij}}{\partial x_j}, \, 0\le i<\text{dim}$, whereas 1535 * its gradient is $G_{ijk} = \frac{\partial T_{ij}}{\partial x_k}$. 1536 * 1537 * You get an object of this type if you apply a FEValuesExtractors::Tensor 1538 * to an FEValues, FEFaceValues or FESubfaceValues object. 1539 * 1540 * @ingroup feaccess vector_valued 1541 */ 1542 template <int dim, int spacedim> 1543 class Tensor<2, dim, spacedim> 1544 { 1545 public: 1546 /** 1547 * Data type for what you get when you apply an extractor of this kind to 1548 * a vector-valued finite element. 1549 */ 1550 using value_type = dealii::Tensor<2, spacedim>; 1551 1552 /** 1553 * Data type for taking the divergence of a tensor: a vector. 1554 */ 1555 using divergence_type = dealii::Tensor<1, spacedim>; 1556 1557 /** 1558 * Data type for taking the gradient of a second order tensor: a third order 1559 * tensor. 1560 */ 1561 using gradient_type = dealii::Tensor<3, spacedim>; 1562 1563 /** 1564 * A struct that provides the output type for the product of the value 1565 * and derivatives of basis functions of the Tensor view and any @p Number type. 1566 */ 1567 template <typename Number> 1568 struct OutputType 1569 { 1570 /** 1571 * An alias for the data type of the product of a @p Number and the 1572 * values of the view the Tensor class. 1573 */ 1574 using value_type = typename ProductType< 1575 Number, 1576 typename Tensor<2, dim, spacedim>::value_type>::type; 1577 1578 /** 1579 * An alias for the data type of the product of a @p Number and the 1580 * divergences of the view the Tensor class. 1581 */ 1582 using divergence_type = typename ProductType< 1583 Number, 1584 typename Tensor<2, dim, spacedim>::divergence_type>::type; 1585 1586 /** 1587 * An alias for the data type of the product of a @p Number and the 1588 * gradient of the view the Tensor class. 1589 */ 1590 using gradient_type = typename ProductType< 1591 Number, 1592 typename Tensor<2, dim, spacedim>::gradient_type>::type; 1593 }; 1594 1595 /** 1596 * A structure where for each shape function we pre-compute a bunch of 1597 * data that will make later accesses much cheaper. 1598 */ 1599 struct ShapeFunctionData 1600 { 1601 /** 1602 * For each pair (shape function,component within vector), store whether 1603 * the selected vector component may be nonzero. For primitive shape 1604 * functions we know for sure whether a certain scalar component of a 1605 * given shape function is nonzero, whereas for non-primitive shape 1606 * functions this may not be entirely clear (e.g. for RT elements it 1607 * depends on the shape of a cell). 1608 */ 1609 bool is_nonzero_shape_function_component 1610 [value_type::n_independent_components]; 1611 1612 /** 1613 * For each pair (shape function, component within vector), store the 1614 * row index within the shape_values, shape_gradients, and 1615 * shape_hessians tables (the column index is the quadrature point 1616 * index). If the shape function is primitive, then we can get this 1617 * information from the shape_function_to_row_table of the FEValues 1618 * object; otherwise, we have to work a bit harder to compute this 1619 * information. 1620 */ 1621 unsigned int row_index[value_type::n_independent_components]; 1622 1623 /** 1624 * For each shape function say the following: if only a single entry in 1625 * is_nonzero_shape_function_component for this shape function is 1626 * nonzero, then store the corresponding value of row_index and 1627 * single_nonzero_component_index represents the index between 0 and 1628 * (dim^2) for which it is attained. If multiple components are nonzero, 1629 * then store -1. If no components are nonzero then store -2. 1630 */ 1631 int single_nonzero_component; 1632 1633 /** 1634 * Index of the @p single_nonzero_component . 1635 */ 1636 unsigned int single_nonzero_component_index; 1637 }; 1638 1639 /** 1640 * Default constructor. Creates an invalid object. 1641 */ 1642 Tensor(); 1643 1644 /** 1645 * Constructor for an object that represents <code>(dim*dim)</code> 1646 * components of a FEValuesBase object (or of one of the classes derived 1647 * from FEValuesBase), representing the unique components comprising a 1648 * second-order tensor valued variable. 1649 * 1650 * The second argument denotes the index of the first component of the 1651 * selected symmetric second order tensor. 1652 */ 1653 Tensor(const FEValuesBase<dim, spacedim> &fe_values_base, 1654 const unsigned int first_tensor_component); 1655 1656 1657 /** 1658 * Copy operator. This is not a lightweight object so we don't allow 1659 * copying and generate a compile-time error if this function is called. 1660 */ 1661 Tensor & 1662 operator=(const Tensor<2, dim, spacedim> &) = delete; 1663 1664 /** 1665 * Return the value of the vector components selected by this view, for 1666 * the shape function and quadrature point selected by the arguments. 1667 * Here, since the view represents a vector-valued part of the FEValues 1668 * object with <code>(dim*dim)</code> components (the unique components of 1669 * a second-order tensor), the return type is a tensor of rank 2. 1670 * 1671 * @param shape_function Number of the shape function to be evaluated. 1672 * Note that this number runs from zero to dofs_per_cell, even in the case 1673 * of an FEFaceValues or FESubfaceValues object. 1674 * 1675 * @param q_point Number of the quadrature point at which function is to 1676 * be evaluated. 1677 * 1678 * @dealiiRequiresUpdateFlags{update_values} 1679 */ 1680 value_type 1681 value(const unsigned int shape_function, const unsigned int q_point) const; 1682 1683 /** 1684 * Return the vector divergence of the vector components selected by this 1685 * view, for the shape function and quadrature point selected by the 1686 * arguments. 1687 * 1688 * See the general discussion of this class for a definition of the 1689 * divergence. 1690 * 1691 * @note The meaning of the arguments is as documented for the value() 1692 * function. 1693 * 1694 * @dealiiRequiresUpdateFlags{update_gradients} 1695 */ 1696 divergence_type 1697 divergence(const unsigned int shape_function, 1698 const unsigned int q_point) const; 1699 1700 /** 1701 * Return the gradient (3-rd order tensor) of the vector components selected 1702 * by this view, for the shape function and quadrature point selected by the 1703 * arguments. 1704 * 1705 * See the general discussion of this class for a definition of the 1706 * gradient. 1707 * 1708 * @note The meaning of the arguments is as documented for the value() 1709 * function. 1710 * 1711 * @dealiiRequiresUpdateFlags{update_gradients} 1712 */ 1713 gradient_type 1714 gradient(const unsigned int shape_function, 1715 const unsigned int q_point) const; 1716 1717 /** 1718 * Return the values of the selected vector components of the finite 1719 * element function characterized by <tt>fe_function</tt> at the 1720 * quadrature points of the cell, face or subface selected the last time 1721 * the <tt>reinit</tt> function of the FEValues object was called. 1722 * 1723 * This function is the equivalent of the 1724 * FEValuesBase::get_function_values function but it only works on the 1725 * selected vector components. 1726 * 1727 * The data type stored by the output vector must be what you get when you 1728 * multiply the values of shape functions (i.e., @p value_type) times the 1729 * type used to store the values of the unknowns $U_j$ of your finite 1730 * element vector $U$ (represented by the @p fe_function argument). 1731 * 1732 * @dealiiRequiresUpdateFlags{update_values} 1733 */ 1734 template <class InputVector> 1735 void 1736 get_function_values( 1737 const InputVector &fe_function, 1738 std::vector<typename ProductType<value_type, 1739 typename InputVector::value_type>::type> 1740 &values) const; 1741 1742 /** 1743 * Same as above, but using a vector of local degree-of-freedom values. 1744 * 1745 * The @p dof_values vector must have a length equal to number of DoFs on 1746 * a cell, and each entry @p dof_values[i] is the value of the local DoF 1747 * @p i. The fundamental prerequisite for the @p InputVector is that it must 1748 * be possible to create an ArrayView from it; this is satisfied by the 1749 * @p std::vector class. 1750 * 1751 * The DoF values typically would be obtained in the following way: 1752 * @code 1753 * Vector<double> local_dof_values(cell->get_fe().n_dofs_per_cell()); 1754 * cell->get_dof_values(solution, local_dof_values); 1755 * @endcode 1756 * or, for a generic @p Number type, 1757 * @code 1758 * std::vector<Number> local_dof_values(cell->get_fe().n_dofs_per_cell()); 1759 * cell->get_dof_values(solution, 1760 * local_dof_values.begin(), 1761 * local_dof_values.end()); 1762 * @endcode 1763 */ 1764 template <class InputVector> 1765 void 1766 get_function_values_from_local_dof_values( 1767 const InputVector &dof_values, 1768 std::vector< 1769 typename OutputType<typename InputVector::value_type>::value_type> 1770 &values) const; 1771 1772 /** 1773 * Return the divergence of the selected vector components of the finite 1774 * element function characterized by <tt>fe_function</tt> at the 1775 * quadrature points of the cell, face or subface selected the last time 1776 * the <tt>reinit</tt> function of the FEValues object was called. 1777 * 1778 * There is no equivalent function such as 1779 * FEValuesBase::get_function_divergences in the FEValues classes but the 1780 * information can be obtained from FEValuesBase::get_function_gradients, 1781 * of course. 1782 * 1783 * See the general discussion of this class for a definition of the 1784 * divergence. 1785 * 1786 * The data type stored by the output vector must be what you get when you 1787 * multiply the divergences of shape functions (i.e., @p divergence_type) 1788 * times the type used to store the values of the unknowns $U_j$ of your 1789 * finite element vector $U$ (represented by the @p fe_function argument). 1790 * 1791 * @dealiiRequiresUpdateFlags{update_gradients} 1792 */ 1793 template <class InputVector> 1794 void 1795 get_function_divergences( 1796 const InputVector &fe_function, 1797 std::vector<typename ProductType<divergence_type, 1798 typename InputVector::value_type>::type> 1799 &divergences) const; 1800 1801 /** 1802 * @copydoc FEValuesViews::Tensor<2,dim,spacedim>::get_function_values_from_local_dof_values() 1803 */ 1804 template <class InputVector> 1805 void 1806 get_function_divergences_from_local_dof_values( 1807 const InputVector &dof_values, 1808 std::vector< 1809 typename OutputType<typename InputVector::value_type>::divergence_type> 1810 &divergences) const; 1811 1812 /** 1813 * Return the gradient of the selected vector components of the finite 1814 * element function characterized by <tt>fe_function</tt> at the 1815 * quadrature points of the cell, face or subface selected the last time 1816 * the <tt>reinit</tt> function of the FEValues object was called. 1817 * 1818 * See the general discussion of this class for a definition of the 1819 * gradient. 1820 * 1821 * The data type stored by the output vector must be what you get when you 1822 * multiply the gradients of shape functions (i.e., @p gradient_type) 1823 * times the type used to store the values of the unknowns $U_j$ of your 1824 * finite element vector $U$ (represented by the @p fe_function argument). 1825 * 1826 * @dealiiRequiresUpdateFlags{update_gradients} 1827 */ 1828 template <class InputVector> 1829 void 1830 get_function_gradients( 1831 const InputVector &fe_function, 1832 std::vector<typename ProductType<gradient_type, 1833 typename InputVector::value_type>::type> 1834 &gradients) const; 1835 1836 /** 1837 * @copydoc FEValuesViews::Tensor<2,dim,spacedim>::get_function_values_from_local_dof_values() 1838 */ 1839 template <class InputVector> 1840 void 1841 get_function_gradients_from_local_dof_values( 1842 const InputVector &dof_values, 1843 std::vector< 1844 typename OutputType<typename InputVector::value_type>::gradient_type> 1845 &gradients) const; 1846 1847 private: 1848 /** 1849 * A pointer to the FEValuesBase object we operate on. 1850 */ 1851 const SmartPointer<const FEValuesBase<dim, spacedim>> fe_values; 1852 1853 /** 1854 * The first component of the vector this view represents of the 1855 * FEValuesBase object. 1856 */ 1857 const unsigned int first_tensor_component; 1858 1859 /** 1860 * Store the data about shape functions. 1861 */ 1862 std::vector<ShapeFunctionData> shape_function_data; 1863 }; 1864 1865 } // namespace FEValuesViews 1866 1867 1868 namespace internal 1869 { 1870 namespace FEValuesViews 1871 { 1872 /** 1873 * A class whose specialization is used to define what FEValuesViews 1874 * object corresponds to the given FEValuesExtractors object. 1875 */ 1876 template <int dim, int spacedim, typename Extractor> 1877 struct ViewType 1878 {}; 1879 1880 /** 1881 * A class whose specialization is used to define what FEValuesViews 1882 * object corresponds to the given FEValuesExtractors object. 1883 * 1884 * When using FEValuesExtractors::Scalar, the corresponding view is an 1885 * FEValuesViews::Scalar<dim, spacedim>. 1886 */ 1887 template <int dim, int spacedim> 1888 struct ViewType<dim, spacedim, FEValuesExtractors::Scalar> 1889 { 1890 using type = typename dealii::FEValuesViews::Scalar<dim, spacedim>; 1891 }; 1892 1893 /** 1894 * A class whose specialization is used to define what FEValuesViews 1895 * object corresponds to the given FEValuesExtractors object. 1896 * 1897 * When using FEValuesExtractors::Vector, the corresponding view is an 1898 * FEValuesViews::Vector<dim, spacedim>. 1899 */ 1900 template <int dim, int spacedim> 1901 struct ViewType<dim, spacedim, FEValuesExtractors::Vector> 1902 { 1903 using type = typename dealii::FEValuesViews::Vector<dim, spacedim>; 1904 }; 1905 1906 /** 1907 * A class whose specialization is used to define what FEValuesViews 1908 * object corresponds to the given FEValuesExtractors object. 1909 * 1910 * When using FEValuesExtractors::Tensor<rank>, the corresponding view is an 1911 * FEValuesViews::Tensor<rank, dim, spacedim>. 1912 */ 1913 template <int dim, int spacedim, int rank> 1914 struct ViewType<dim, spacedim, FEValuesExtractors::Tensor<rank>> 1915 { 1916 using type = typename dealii::FEValuesViews::Tensor<rank, dim, spacedim>; 1917 }; 1918 1919 /** 1920 * A class whose specialization is used to define what FEValuesViews 1921 * object corresponds to the given FEValuesExtractors object. 1922 * 1923 * When using FEValuesExtractors::SymmetricTensor<rank>, the corresponding 1924 * view is an FEValuesViews::SymmetricTensor<rank, dim, spacedim>. 1925 */ 1926 template <int dim, int spacedim, int rank> 1927 struct ViewType<dim, spacedim, FEValuesExtractors::SymmetricTensor<rank>> 1928 { 1929 using type = 1930 typename dealii::FEValuesViews::SymmetricTensor<rank, dim, spacedim>; 1931 }; 1932 1933 /** 1934 * A class objects of which store a collection of FEValuesViews::Scalar, 1935 * FEValuesViews::Vector, etc object. The FEValuesBase class uses it to 1936 * generate all possible Views classes upon construction time; we do this 1937 * at construction time since the Views classes cache some information and 1938 * are therefore relatively expensive to create. 1939 */ 1940 template <int dim, int spacedim> 1941 struct Cache 1942 { 1943 /** 1944 * Caches for scalar and vector, and symmetric second-order tensor 1945 * valued views. 1946 */ 1947 std::vector<dealii::FEValuesViews::Scalar<dim, spacedim>> scalars; 1948 std::vector<dealii::FEValuesViews::Vector<dim, spacedim>> vectors; 1949 std::vector<dealii::FEValuesViews::SymmetricTensor<2, dim, spacedim>> 1950 symmetric_second_order_tensors; 1951 std::vector<dealii::FEValuesViews::Tensor<2, dim, spacedim>> 1952 second_order_tensors; 1953 1954 /** 1955 * Constructor. 1956 */ 1957 Cache(const FEValuesBase<dim, spacedim> &fe_values); 1958 }; 1959 } // namespace FEValuesViews 1960 } // namespace internal 1961 1962 namespace FEValuesViews 1963 { 1964 /** 1965 * A templated alias that associates to a given Extractor class 1966 * the corresponding view in FEValuesViews. 1967 */ 1968 template <int dim, int spacedim, typename Extractor> 1969 using View = typename dealii::internal::FEValuesViews:: 1970 ViewType<dim, spacedim, Extractor>::type; 1971 } // namespace FEValuesViews 1972 1973 1974 /** 1975 * FEValues, FEFaceValues and FESubfaceValues objects are interfaces to finite 1976 * element and mapping classes on the one hand side, to cells and quadrature 1977 * rules on the other side. They allow to evaluate values or derivatives of 1978 * shape functions at the quadrature points of a quadrature formula when 1979 * projected by a mapping from the unit cell onto a cell in real space. The 1980 * reason for this abstraction is possible optimization: Depending on the type 1981 * of finite element and mapping, some values can be computed once on the unit 1982 * cell. Others must be computed on each cell, but maybe computation of 1983 * several values at the same time offers ways for optimization. Since this 1984 * interplay may be complex and depends on the actual finite element, it 1985 * cannot be left to the applications programmer. 1986 * 1987 * FEValues, FEFaceValues and FESubfaceValues provide only data handling: 1988 * computations are left to objects of type Mapping and FiniteElement. These 1989 * provide functions <tt>get_*_data</tt> and <tt>fill_*_values</tt> which are 1990 * called by the constructor and <tt>reinit</tt> functions of 1991 * <tt>FEValues*</tt>, respectively. 1992 * 1993 * <h3>General usage</h3> 1994 * 1995 * Usually, an object of <tt>FEValues*</tt> is used in integration loops over 1996 * all cells of a triangulation (or faces of cells). To take full advantage of 1997 * the optimization features, it should be constructed before the loop so that 1998 * information that does not depend on the location and shape of cells can be 1999 * computed once and for all (this includes, for example, the values of shape 2000 * functions at quadrature points for the most common elements: we can 2001 * evaluate them on the unit cell and they will be the same when mapped to the 2002 * real cell). Then, in the loop over all cells, it must be re-initialized for 2003 * each grid cell to compute that part of the information that changes 2004 * depending on the actual cell (for example, the gradient of shape functions 2005 * equals the gradient on the unit cell -- which can be computed once and for 2006 * all -- times the Jacobian matrix of the mapping between unit and real cell, 2007 * which needs to be recomputed for each cell). 2008 * 2009 * A typical piece of code, adding up local contributions to the Laplace 2010 * matrix looks like this: 2011 * 2012 * @code 2013 * FEValues values (mapping, finite_element, quadrature, flags); 2014 * for (const auto &cell : dof_handler.active_cell_iterators()) 2015 * { 2016 * values.reinit(cell); 2017 * for (unsigned int q=0; q<quadrature.size(); ++q) 2018 * for (unsigned int i=0; i<finite_element.n_dofs_per_cell(); ++i) 2019 * for (unsigned int j=0; j<finite_element.n_dofs_per_cell(); ++j) 2020 * A(i,j) += fe_values.shape_value(i,q) * 2021 * fe_values.shape_value(j,q) * 2022 * fe_values.JxW(q); 2023 * ... 2024 * } 2025 * @endcode 2026 * 2027 * The individual functions used here are described below. Note that by 2028 * design, the order of quadrature points used inside the FEValues object is 2029 * the same as defined by the quadrature formula passed to the constructor of 2030 * the FEValues object above. 2031 * 2032 * <h3>Member functions</h3> 2033 * 2034 * The functions of this class fall into different categories: 2035 * <ul> 2036 * <li> shape_value(), shape_grad(), etc: return one of the values of this 2037 * object at a time. These functions are inlined, so this is the suggested 2038 * access to all finite element values. There should be no loss in performance 2039 * with an optimizing compiler. If the finite element is vector valued, then 2040 * these functions return the only non-zero component of the requested shape 2041 * function. However, some finite elements have shape functions that have more 2042 * than one non-zero component (we call them non-"primitive"), and in this 2043 * case this set of functions will throw an exception since they cannot 2044 * generate a useful result. Rather, use the next set of functions. 2045 * 2046 * <li> shape_value_component(), shape_grad_component(), etc: This is the same 2047 * set of functions as above, except that for vector valued finite elements 2048 * they return only one vector component. This is useful for elements of which 2049 * shape functions have more than one non-zero component, since then the above 2050 * functions cannot be used, and you have to walk over all (or only the non- 2051 * zero) components of the shape function using this set of functions. 2052 * 2053 * <li> get_function_values(), get_function_gradients(), etc.: Compute a 2054 * finite element function or its derivative in quadrature points. 2055 * 2056 * <li> reinit: initialize the FEValues object for a certain cell. This 2057 * function is not in the present class but only in the derived classes and 2058 * has a variable call syntax. See the docs for the derived classes for more 2059 * information. 2060 * </ul> 2061 * 2062 * 2063 * <h3>Internals about the implementation</h3> 2064 * 2065 * The mechanisms by which this class work are discussed on the page on 2066 * @ref UpdateFlags "Update flags" 2067 * and about the 2068 * @ref FE_vs_Mapping_vs_FEValues "How Mapping, FiniteElement, and FEValues work together". 2069 * 2070 * 2071 * @ingroup feaccess 2072 */ 2073 template <int dim, int spacedim> 2074 class FEValuesBase : public Subscriptor 2075 { 2076 public: 2077 /** 2078 * Dimension in which this object operates. 2079 */ 2080 static const unsigned int dimension = dim; 2081 2082 /** 2083 * Dimension of the space in which this object operates. 2084 */ 2085 static const unsigned int space_dimension = spacedim; 2086 2087 /** 2088 * Number of quadrature points. 2089 */ 2090 const unsigned int n_quadrature_points; 2091 2092 /** 2093 * Number of shape functions per cell. If we use this base class to evaluate 2094 * a finite element on faces of cells, this is still the number of degrees 2095 * of freedom per cell, not per face. 2096 */ 2097 const unsigned int dofs_per_cell; 2098 2099 2100 /** 2101 * Constructor. Set up the array sizes with <tt>n_q_points</tt> quadrature 2102 * points, <tt>dofs_per_cell</tt> trial functions per cell and with the 2103 * given pattern to update the fields when the <tt>reinit</tt> function of 2104 * the derived classes is called. The fields themselves are not set up, this 2105 * must happen in the constructor of the derived class. 2106 */ 2107 FEValuesBase(const unsigned int n_q_points, 2108 const unsigned int dofs_per_cell, 2109 const UpdateFlags update_flags, 2110 const Mapping<dim, spacedim> & mapping, 2111 const FiniteElement<dim, spacedim> &fe); 2112 2113 /** 2114 * The copy assignment is deleted since objects of this class are not 2115 * copyable. 2116 */ 2117 FEValuesBase & 2118 operator=(const FEValuesBase &) = delete; 2119 2120 /** 2121 * The copy constructor is deleted since objects of this class are not 2122 * copyable. 2123 */ 2124 FEValuesBase(const FEValuesBase &) = delete; 2125 2126 /** 2127 * Destructor. 2128 */ 2129 virtual ~FEValuesBase() override; 2130 2131 2132 /// @name Access to shape function values 2133 /// 2134 /// These fields are filled by the finite element. 2135 //@{ 2136 2137 /** 2138 * Value of a shape function at a quadrature point on the cell, face or 2139 * subface selected the last time the <tt>reinit</tt> function of the 2140 * derived class was called. 2141 * 2142 * If the shape function is vector-valued, then this returns the only non- 2143 * zero component. If the shape function has more than one non-zero 2144 * component (i.e. it is not primitive), then throw an exception of type 2145 * ExcShapeFunctionNotPrimitive. In that case, use the 2146 * shape_value_component() function. 2147 * 2148 * @param function_no Number of the shape function to be evaluated. Note 2149 * that this number runs from zero to dofs_per_cell, even in the case of an 2150 * FEFaceValues or FESubfaceValues object. 2151 * 2152 * @param point_no Number of the quadrature point at which function is to be 2153 * evaluated 2154 * 2155 * @dealiiRequiresUpdateFlags{update_values} 2156 */ 2157 const double & 2158 shape_value(const unsigned int function_no, 2159 const unsigned int point_no) const; 2160 2161 /** 2162 * Compute one vector component of the value of a shape function at a 2163 * quadrature point. If the finite element is scalar, then only component 2164 * zero is allowed and the return value equals that of the shape_value() 2165 * function. If the finite element is vector valued but all shape functions 2166 * are primitive (i.e. they are non-zero in only one component), then the 2167 * value returned by shape_value() equals that of this function for exactly 2168 * one component. This function is therefore only of greater interest if the 2169 * shape function is not primitive, but then it is necessary since the other 2170 * function cannot be used. 2171 * 2172 * @param function_no Number of the shape function to be evaluated. 2173 * 2174 * @param point_no Number of the quadrature point at which function is to be 2175 * evaluated. 2176 * 2177 * @param component vector component to be evaluated. 2178 * 2179 * @dealiiRequiresUpdateFlags{update_values} 2180 */ 2181 double 2182 shape_value_component(const unsigned int function_no, 2183 const unsigned int point_no, 2184 const unsigned int component) const; 2185 2186 /** 2187 * Compute the gradient of the <tt>function_no</tt>th shape function at the 2188 * <tt>quadrature_point</tt>th quadrature point with respect to real cell 2189 * coordinates. If you want to get the derivative in one of the coordinate 2190 * directions, use the appropriate function of the Tensor class to extract 2191 * one component of the Tensor returned by this function. Since only a 2192 * reference to the gradient's value is returned, there should be no major 2193 * performance drawback. 2194 * 2195 * If the shape function is vector-valued, then this returns the only non- 2196 * zero component. If the shape function has more than one non-zero 2197 * component (i.e. it is not primitive), then it will throw an exception of 2198 * type ExcShapeFunctionNotPrimitive. In that case, use the 2199 * shape_grad_component() function. 2200 * 2201 * The same holds for the arguments of this function as for the 2202 * shape_value() function. 2203 * 2204 * @param function_no Number of the shape function to be evaluated. 2205 * 2206 * @param quadrature_point Number of the quadrature point at which function 2207 * is to be evaluated. 2208 * 2209 * @dealiiRequiresUpdateFlags{update_gradients} 2210 */ 2211 const Tensor<1, spacedim> & 2212 shape_grad(const unsigned int function_no, 2213 const unsigned int quadrature_point) const; 2214 2215 /** 2216 * Return one vector component of the gradient of a shape function at a 2217 * quadrature point. If the finite element is scalar, then only component 2218 * zero is allowed and the return value equals that of the shape_grad() 2219 * function. If the finite element is vector valued but all shape functions 2220 * are primitive (i.e. they are non-zero in only one component), then the 2221 * value returned by shape_grad() equals that of this function for exactly 2222 * one component. This function is therefore only of greater interest if the 2223 * shape function is not primitive, but then it is necessary since the other 2224 * function cannot be used. 2225 * 2226 * The same holds for the arguments of this function as for the 2227 * shape_value_component() function. 2228 * 2229 * @dealiiRequiresUpdateFlags{update_gradients} 2230 */ 2231 Tensor<1, spacedim> 2232 shape_grad_component(const unsigned int function_no, 2233 const unsigned int point_no, 2234 const unsigned int component) const; 2235 2236 /** 2237 * Second derivatives of the <tt>function_no</tt>th shape function at the 2238 * <tt>point_no</tt>th quadrature point with respect to real cell 2239 * coordinates. If you want to get the derivatives in one of the coordinate 2240 * directions, use the appropriate function of the Tensor class to extract 2241 * one component. Since only a reference to the hessian values is returned, 2242 * there should be no major performance drawback. 2243 * 2244 * If the shape function is vector-valued, then this returns the only non- 2245 * zero component. If the shape function has more than one non-zero 2246 * component (i.e. it is not primitive), then throw an exception of type 2247 * ExcShapeFunctionNotPrimitive. In that case, use the 2248 * shape_hessian_component() function. 2249 * 2250 * The same holds for the arguments of this function as for the 2251 * shape_value() function. 2252 * 2253 * @dealiiRequiresUpdateFlags{update_hessians} 2254 */ 2255 const Tensor<2, spacedim> & 2256 shape_hessian(const unsigned int function_no, 2257 const unsigned int point_no) const; 2258 2259 /** 2260 * Return one vector component of the hessian of a shape function at a 2261 * quadrature point. If the finite element is scalar, then only component 2262 * zero is allowed and the return value equals that of the shape_hessian() 2263 * function. If the finite element is vector valued but all shape functions 2264 * are primitive (i.e. they are non-zero in only one component), then the 2265 * value returned by shape_hessian() equals that of this function for 2266 * exactly one component. This function is therefore only of greater 2267 * interest if the shape function is not primitive, but then it is necessary 2268 * since the other function cannot be used. 2269 * 2270 * The same holds for the arguments of this function as for the 2271 * shape_value_component() function. 2272 * 2273 * @dealiiRequiresUpdateFlags{update_hessians} 2274 */ 2275 Tensor<2, spacedim> 2276 shape_hessian_component(const unsigned int function_no, 2277 const unsigned int point_no, 2278 const unsigned int component) const; 2279 2280 /** 2281 * Third derivatives of the <tt>function_no</tt>th shape function at the 2282 * <tt>point_no</tt>th quadrature point with respect to real cell 2283 * coordinates. If you want to get the 3rd derivatives in one of the 2284 * coordinate directions, use the appropriate function of the Tensor class 2285 * to extract one component. Since only a reference to the 3rd derivative 2286 * values is returned, there should be no major performance drawback. 2287 * 2288 * If the shape function is vector-valued, then this returns the only non- 2289 * zero component. If the shape function has more than one non-zero 2290 * component (i.e. it is not primitive), then throw an exception of type 2291 * ExcShapeFunctionNotPrimitive. In that case, use the 2292 * shape_3rdderivative_component() function. 2293 * 2294 * The same holds for the arguments of this function as for the 2295 * shape_value() function. 2296 * 2297 * @dealiiRequiresUpdateFlags{update_3rd_derivatives} 2298 */ 2299 const Tensor<3, spacedim> & 2300 shape_3rd_derivative(const unsigned int function_no, 2301 const unsigned int point_no) const; 2302 2303 /** 2304 * Return one vector component of the third derivative of a shape function 2305 * at a quadrature point. If the finite element is scalar, then only 2306 * component zero is allowed and the return value equals that of the 2307 * shape_3rdderivative() function. If the finite element is vector valued 2308 * but all shape functions are primitive (i.e. they are non-zero in only one 2309 * component), then the value returned by shape_3rdderivative() equals that 2310 * of this function for exactly one component. This function is therefore 2311 * only of greater interest if the shape function is not primitive, but then 2312 * it is necessary since the other function cannot be used. 2313 * 2314 * The same holds for the arguments of this function as for the 2315 * shape_value_component() function. 2316 * 2317 * @dealiiRequiresUpdateFlags{update_3rd_derivatives} 2318 */ 2319 Tensor<3, spacedim> 2320 shape_3rd_derivative_component(const unsigned int function_no, 2321 const unsigned int point_no, 2322 const unsigned int component) const; 2323 2324 //@} 2325 /// @name Access to values of global finite element fields 2326 //@{ 2327 2328 /** 2329 * Return the values of a finite element function restricted to the current 2330 * cell, face or subface selected the last time the <tt>reinit</tt> function 2331 * of the derived class was called, at the quadrature points. 2332 * 2333 * If the present cell is not active then values are interpolated to the 2334 * current cell and point values are computed from that. 2335 * 2336 * This function may only be used if the finite element in use is a scalar 2337 * one, i.e. has only one vector component. To get values of multi- 2338 * component elements, there is another get_function_values() below, 2339 * returning a vector of vectors of results. 2340 * 2341 * @param[in] fe_function A vector of values that describes (globally) the 2342 * finite element function that this function should evaluate at the 2343 * quadrature points of the current cell. 2344 * 2345 * @param[out] values The values of the function specified by fe_function at 2346 * the quadrature points of the current cell. The object is assume to 2347 * already have the correct size. The data type stored by this output vector 2348 * must be what you get when you multiply the values of shape function times 2349 * the type used to store the values of the unknowns $U_j$ of your finite 2350 * element vector $U$ (represented by the @p fe_function argument). This 2351 * happens to be equal to the type of the elements of the solution vector. 2352 * 2353 * @post <code>values[q]</code> will contain the value of the field 2354 * described by fe_function at the $q$th quadrature point. 2355 * 2356 * @note The actual data type of the input vector may be either a 2357 * Vector<T>, BlockVector<T>, or one of the PETSc or Trilinos 2358 * vector wrapper classes. It represents a global vector of DoF values 2359 * associated with the DoFHandler object with which this FEValues object was 2360 * last initialized. 2361 * 2362 * @dealiiRequiresUpdateFlags{update_values} 2363 */ 2364 template <class InputVector> 2365 void 2366 get_function_values( 2367 const InputVector & fe_function, 2368 std::vector<typename InputVector::value_type> &values) const; 2369 2370 /** 2371 * This function does the same as the other get_function_values(), but 2372 * applied to multi-component (vector-valued) elements. The meaning of the 2373 * arguments is as explained there. 2374 * 2375 * @post <code>values[q]</code> is a vector of values of the field described 2376 * by fe_function at the $q$th quadrature point. The size of the vector 2377 * accessed by <code>values[q]</code> equals the number of components of the 2378 * finite element, i.e. <code>values[q](c)</code> returns the value of the 2379 * $c$th vector component at the $q$th quadrature point. 2380 * 2381 * @dealiiRequiresUpdateFlags{update_values} 2382 */ 2383 template <class InputVector> 2384 void 2385 get_function_values( 2386 const InputVector & fe_function, 2387 std::vector<Vector<typename InputVector::value_type>> &values) const; 2388 2389 /** 2390 * Generate function values from an arbitrary vector. 2391 * 2392 * This function offers the possibility to extract function values in 2393 * quadrature points from vectors not corresponding to a whole 2394 * discretization. 2395 * 2396 * The vector <tt>indices</tt> corresponds to the degrees of freedom on a 2397 * single cell. Its length may even be a multiple of the number of dofs per 2398 * cell. Then, the vectors in <tt>value</tt> should allow for the same 2399 * multiple of the components of the finite element. 2400 * 2401 * You may want to use this function, if you want to access just a single 2402 * block from a BlockVector, if you have a multi-level vector or if you 2403 * already have a local representation of your finite element data. 2404 * 2405 * @dealiiRequiresUpdateFlags{update_values} 2406 */ 2407 template <class InputVector> 2408 void 2409 get_function_values( 2410 const InputVector & fe_function, 2411 const ArrayView<const types::global_dof_index> &indices, 2412 std::vector<typename InputVector::value_type> & values) const; 2413 2414 /** 2415 * Generate vector function values from an arbitrary vector. 2416 * 2417 * This function offers the possibility to extract function values in 2418 * quadrature points from vectors not corresponding to a whole 2419 * discretization. 2420 * 2421 * The vector <tt>indices</tt> corresponds to the degrees of freedom on a 2422 * single cell. Its length may even be a multiple of the number of dofs per 2423 * cell. Then, the vectors in <tt>value</tt> should allow for the same 2424 * multiple of the components of the finite element. 2425 * 2426 * You may want to use this function, if you want to access just a single 2427 * block from a BlockVector, if you have a multi-level vector or if you 2428 * already have a local representation of your finite element data. 2429 * 2430 * Since this function allows for fairly general combinations of argument 2431 * sizes, be aware that the checks on the arguments may not detect errors. 2432 * 2433 * @dealiiRequiresUpdateFlags{update_values} 2434 */ 2435 template <class InputVector> 2436 void 2437 get_function_values( 2438 const InputVector & fe_function, 2439 const ArrayView<const types::global_dof_index> & indices, 2440 std::vector<Vector<typename InputVector::value_type>> &values) const; 2441 2442 2443 /** 2444 * Generate vector function values from an arbitrary vector. 2445 * 2446 * This function offers the possibility to extract function values in 2447 * quadrature points from vectors not corresponding to a whole 2448 * discretization. 2449 * 2450 * The vector <tt>indices</tt> corresponds to the degrees of freedom on a 2451 * single cell. Its length may even be a multiple of the number of dofs per 2452 * cell. Then, the vectors in <tt>value</tt> should allow for the same 2453 * multiple of the components of the finite element. 2454 * 2455 * Depending on the value of the last argument, the outer vector of 2456 * <tt>values</tt> has either the length of the quadrature rule 2457 * (<tt>quadrature_points_fastest == false</tt>) or the length of components 2458 * to be filled <tt>quadrature_points_fastest == true</tt>. If <tt>p</tt> is 2459 * the current quadrature point number and <tt>i</tt> is the vector 2460 * component of the solution desired, the access to <tt>values</tt> is 2461 * <tt>values[p][i]</tt> if <tt>quadrature_points_fastest == false</tt>, and 2462 * <tt>values[i][p]</tt> otherwise. 2463 * 2464 * You may want to use this function, if you want to access just a single 2465 * block from a BlockVector, if you have a multi-level vector or if you 2466 * already have a local representation of your finite element data. 2467 * 2468 * Since this function allows for fairly general combinations of argument 2469 * sizes, be aware that the checks on the arguments may not detect errors. 2470 * 2471 * @dealiiRequiresUpdateFlags{update_values} 2472 */ 2473 template <class InputVector> 2474 void 2475 get_function_values( 2476 const InputVector & fe_function, 2477 const ArrayView<const types::global_dof_index> & indices, 2478 ArrayView<std::vector<typename InputVector::value_type>> values, 2479 const bool quadrature_points_fastest) const; 2480 2481 //@} 2482 /// @name Access to derivatives of global finite element fields 2483 //@{ 2484 2485 /** 2486 * Compute the gradients of a finite element at the quadrature points of a 2487 * cell. This function is the equivalent of the corresponding 2488 * get_function_values() function (see there for more information) but 2489 * evaluates the finite element field's gradient instead of its value. 2490 * 2491 * This function may only be used if the finite element in use is a scalar 2492 * one, i.e. has only one vector component. There is a corresponding 2493 * function of the same name for vector-valued finite elements. 2494 * 2495 * @param[in] fe_function A vector of values that describes (globally) the 2496 * finite element function that this function should evaluate at the 2497 * quadrature points of the current cell. 2498 * 2499 * @param[out] gradients The gradients of the function specified by 2500 * fe_function at the quadrature points of the current cell. The gradients 2501 * are computed in real space (as opposed to on the unit cell). The object 2502 * is assume to already have the correct size. The data type stored by this 2503 * output vector must be what you get when you multiply the gradients of 2504 * shape function times the type used to store the values of the unknowns 2505 * $U_j$ of your finite element vector $U$ (represented by the @p 2506 * fe_function argument). 2507 * 2508 * @post <code>gradients[q]</code> will contain the gradient of the field 2509 * described by fe_function at the $q$th quadrature point. 2510 * <code>gradients[q][d]</code> represents the derivative in coordinate 2511 * direction $d$ at quadrature point $q$. 2512 * 2513 * @note The actual data type of the input vector may be either a 2514 * Vector<T>, BlockVector<T>, or one of the PETSc or Trilinos 2515 * vector wrapper classes. It represents a global vector of DoF values 2516 * associated with the DoFHandler object with which this FEValues object was 2517 * last initialized. 2518 * 2519 * @dealiiRequiresUpdateFlags{update_gradients} 2520 */ 2521 template <class InputVector> 2522 void 2523 get_function_gradients( 2524 const InputVector &fe_function, 2525 std::vector<Tensor<1, spacedim, typename InputVector::value_type>> 2526 &gradients) const; 2527 2528 /** 2529 * This function does the same as the other get_function_gradients(), but 2530 * applied to multi-component (vector-valued) elements. The meaning of the 2531 * arguments is as explained there. 2532 * 2533 * @post <code>gradients[q]</code> is a vector of gradients of the field 2534 * described by fe_function at the $q$th quadrature point. The size of the 2535 * vector accessed by <code>gradients[q]</code> equals the number of 2536 * components of the finite element, i.e. <code>gradients[q][c]</code> 2537 * returns the gradient of the $c$th vector component at the $q$th 2538 * quadrature point. Consequently, <code>gradients[q][c][d]</code> is the 2539 * derivative in coordinate direction $d$ of the $c$th vector component of 2540 * the vector field at quadrature point $q$ of the current cell. 2541 * 2542 * @dealiiRequiresUpdateFlags{update_gradients} 2543 */ 2544 template <class InputVector> 2545 void 2546 get_function_gradients( 2547 const InputVector &fe_function, 2548 std::vector< 2549 std::vector<Tensor<1, spacedim, typename InputVector::value_type>>> 2550 &gradients) const; 2551 2552 /** 2553 * Function gradient access with more flexibility. See get_function_values() 2554 * with corresponding arguments. 2555 * 2556 * @dealiiRequiresUpdateFlags{update_gradients} 2557 */ 2558 template <class InputVector> 2559 void 2560 get_function_gradients( 2561 const InputVector & fe_function, 2562 const ArrayView<const types::global_dof_index> &indices, 2563 std::vector<Tensor<1, spacedim, typename InputVector::value_type>> 2564 &gradients) const; 2565 2566 /** 2567 * Function gradient access with more flexibility. See get_function_values() 2568 * with corresponding arguments. 2569 * 2570 * @dealiiRequiresUpdateFlags{update_gradients} 2571 */ 2572 template <class InputVector> 2573 void 2574 get_function_gradients( 2575 const InputVector & fe_function, 2576 const ArrayView<const types::global_dof_index> &indices, 2577 ArrayView< 2578 std::vector<Tensor<1, spacedim, typename InputVector::value_type>>> 2579 gradients, 2580 bool quadrature_points_fastest = false) const; 2581 2582 //@} 2583 /// @name Access to second derivatives 2584 /// 2585 /// Hessian matrices and Laplacians of global finite element fields 2586 //@{ 2587 2588 /** 2589 * Compute the tensor of second derivatives of a finite element at the 2590 * quadrature points of a cell. This function is the equivalent of the 2591 * corresponding get_function_values() function (see there for more 2592 * information) but evaluates the finite element field's second derivatives 2593 * instead of its value. 2594 * 2595 * This function may only be used if the finite element in use is a scalar 2596 * one, i.e. has only one vector component. There is a corresponding 2597 * function of the same name for vector-valued finite elements. 2598 * 2599 * @param[in] fe_function A vector of values that describes (globally) the 2600 * finite element function that this function should evaluate at the 2601 * quadrature points of the current cell. 2602 * 2603 * @param[out] hessians The Hessians of the function specified by 2604 * fe_function at the quadrature points of the current cell. The Hessians 2605 * are computed in real space (as opposed to on the unit cell). The object 2606 * is assume to already have the correct size. The data type stored by this 2607 * output vector must be what you get when you multiply the Hessians of 2608 * shape function times the type used to store the values of the unknowns 2609 * $U_j$ of your finite element vector $U$ (represented by the @p 2610 * fe_function argument). 2611 * 2612 * @post <code>hessians[q]</code> will contain the Hessian of the field 2613 * described by fe_function at the $q$th quadrature point. 2614 * <code>hessians[q][i][j]</code> represents the $(i,j)$th component of the 2615 * matrix of second derivatives at quadrature point $q$. 2616 * 2617 * @note The actual data type of the input vector may be either a 2618 * Vector<T>, BlockVector<T>, or one of the PETSc or Trilinos 2619 * vector wrapper classes. It represents a global vector of DoF values 2620 * associated with the DoFHandler object with which this FEValues object was 2621 * last initialized. 2622 * 2623 * @dealiiRequiresUpdateFlags{update_hessians} 2624 */ 2625 template <class InputVector> 2626 void 2627 get_function_hessians( 2628 const InputVector &fe_function, 2629 std::vector<Tensor<2, spacedim, typename InputVector::value_type>> 2630 &hessians) const; 2631 2632 /** 2633 * This function does the same as the other get_function_hessians(), but 2634 * applied to multi-component (vector-valued) elements. The meaning of the 2635 * arguments is as explained there. 2636 * 2637 * @post <code>hessians[q]</code> is a vector of Hessians of the field 2638 * described by fe_function at the $q$th quadrature point. The size of the 2639 * vector accessed by <code>hessians[q]</code> equals the number of 2640 * components of the finite element, i.e. <code>hessians[q][c]</code> 2641 * returns the Hessian of the $c$th vector component at the $q$th quadrature 2642 * point. Consequently, <code>hessians[q][c][i][j]</code> is the $(i,j)$th 2643 * component of the matrix of second derivatives of the $c$th vector 2644 * component of the vector field at quadrature point $q$ of the current 2645 * cell. 2646 * 2647 * @dealiiRequiresUpdateFlags{update_hessians} 2648 */ 2649 template <class InputVector> 2650 void 2651 get_function_hessians( 2652 const InputVector &fe_function, 2653 std::vector< 2654 std::vector<Tensor<2, spacedim, typename InputVector::value_type>>> 2655 & hessians, 2656 bool quadrature_points_fastest = false) const; 2657 2658 /** 2659 * Access to the second derivatives of a function with more flexibility. See 2660 * get_function_values() with corresponding arguments. 2661 */ 2662 template <class InputVector> 2663 void 2664 get_function_hessians( 2665 const InputVector & fe_function, 2666 const ArrayView<const types::global_dof_index> &indices, 2667 std::vector<Tensor<2, spacedim, typename InputVector::value_type>> 2668 &hessians) const; 2669 2670 /** 2671 * Access to the second derivatives of a function with more flexibility. See 2672 * get_function_values() with corresponding arguments. 2673 * 2674 * @dealiiRequiresUpdateFlags{update_hessians} 2675 */ 2676 template <class InputVector> 2677 void 2678 get_function_hessians( 2679 const InputVector & fe_function, 2680 const ArrayView<const types::global_dof_index> &indices, 2681 ArrayView< 2682 std::vector<Tensor<2, spacedim, typename InputVector::value_type>>> 2683 hessians, 2684 bool quadrature_points_fastest = false) const; 2685 2686 /** 2687 * Compute the (scalar) Laplacian (i.e. the trace of the tensor of second 2688 * derivatives) of a finite element at the quadrature points of a cell. This 2689 * function is the equivalent of the corresponding get_function_values() 2690 * function (see there for more information) but evaluates the finite 2691 * element field's second derivatives instead of its value. 2692 * 2693 * This function may only be used if the finite element in use is a scalar 2694 * one, i.e. has only one vector component. There is a corresponding 2695 * function of the same name for vector-valued finite elements. 2696 * 2697 * @param[in] fe_function A vector of values that describes (globally) the 2698 * finite element function that this function should evaluate at the 2699 * quadrature points of the current cell. 2700 * 2701 * @param[out] laplacians The Laplacians of the function specified by 2702 * fe_function at the quadrature points of the current cell. The Laplacians 2703 * are computed in real space (as opposed to on the unit cell). The object 2704 * is assume to already have the correct size. The data type stored by this 2705 * output vector must be what you get when you multiply the Laplacians of 2706 * shape function times the type used to store the values of the unknowns 2707 * $U_j$ of your finite element vector $U$ (represented by the @p 2708 * fe_function argument). This happens to be equal to the type of the 2709 * elements of the input vector. 2710 * 2711 * @post <code>laplacians[q]</code> will contain the Laplacian of the field 2712 * described by fe_function at the $q$th quadrature point. 2713 * 2714 * @post For each component of the output vector, there holds 2715 * <code>laplacians[q]=trace(hessians[q])</code>, where <tt>hessians</tt> 2716 * would be the output of the get_function_hessians() function. 2717 * 2718 * @note The actual data type of the input vector may be either a 2719 * Vector<T>, BlockVector<T>, or one of the PETSc or Trilinos 2720 * vector wrapper classes. It represents a global vector of DoF values 2721 * associated with the DoFHandler object with which this FEValues object was 2722 * last initialized. 2723 * 2724 * @dealiiRequiresUpdateFlags{update_hessians} 2725 */ 2726 template <class InputVector> 2727 void 2728 get_function_laplacians( 2729 const InputVector & fe_function, 2730 std::vector<typename InputVector::value_type> &laplacians) const; 2731 2732 /** 2733 * This function does the same as the other get_function_laplacians(), but 2734 * applied to multi-component (vector-valued) elements. The meaning of the 2735 * arguments is as explained there. 2736 * 2737 * @post <code>laplacians[q]</code> is a vector of Laplacians of the field 2738 * described by fe_function at the $q$th quadrature point. The size of the 2739 * vector accessed by <code>laplacians[q]</code> equals the number of 2740 * components of the finite element, i.e. <code>laplacians[q][c]</code> 2741 * returns the Laplacian of the $c$th vector component at the $q$th 2742 * quadrature point. 2743 * 2744 * @post For each component of the output vector, there holds 2745 * <code>laplacians[q][c]=trace(hessians[q][c])</code>, where 2746 * <tt>hessians</tt> would be the output of the get_function_hessians() 2747 * function. 2748 * 2749 * @dealiiRequiresUpdateFlags{update_hessians} 2750 */ 2751 template <class InputVector> 2752 void 2753 get_function_laplacians( 2754 const InputVector & fe_function, 2755 std::vector<Vector<typename InputVector::value_type>> &laplacians) const; 2756 2757 /** 2758 * Access to the second derivatives of a function with more flexibility. See 2759 * get_function_values() with corresponding arguments. 2760 * 2761 * @dealiiRequiresUpdateFlags{update_hessians} 2762 */ 2763 template <class InputVector> 2764 void 2765 get_function_laplacians( 2766 const InputVector & fe_function, 2767 const ArrayView<const types::global_dof_index> &indices, 2768 std::vector<typename InputVector::value_type> & laplacians) const; 2769 2770 /** 2771 * Access to the second derivatives of a function with more flexibility. See 2772 * get_function_values() with corresponding arguments. 2773 * 2774 * @dealiiRequiresUpdateFlags{update_hessians} 2775 */ 2776 template <class InputVector> 2777 void 2778 get_function_laplacians( 2779 const InputVector & fe_function, 2780 const ArrayView<const types::global_dof_index> & indices, 2781 std::vector<Vector<typename InputVector::value_type>> &laplacians) const; 2782 2783 /** 2784 * Access to the second derivatives of a function with more flexibility. See 2785 * get_function_values() with corresponding arguments. 2786 * 2787 * @dealiiRequiresUpdateFlags{update_hessians} 2788 */ 2789 template <class InputVector> 2790 void 2791 get_function_laplacians( 2792 const InputVector & fe_function, 2793 const ArrayView<const types::global_dof_index> & indices, 2794 std::vector<std::vector<typename InputVector::value_type>> &laplacians, 2795 bool quadrature_points_fastest = false) const; 2796 2797 //@} 2798 /// @name Access to third derivatives of global finite element fields 2799 //@{ 2800 2801 /** 2802 * Compute the tensor of third derivatives of a finite element at the 2803 * quadrature points of a cell. This function is the equivalent of the 2804 * corresponding get_function_values() function (see there for more 2805 * information) but evaluates the finite element field's third derivatives 2806 * instead of its value. 2807 * 2808 * This function may only be used if the finite element in use is a scalar 2809 * one, i.e. has only one vector component. There is a corresponding 2810 * function of the same name for vector-valued finite elements. 2811 * 2812 * @param[in] fe_function A vector of values that describes (globally) the 2813 * finite element function that this function should evaluate at the 2814 * quadrature points of the current cell. 2815 * 2816 * @param[out] third_derivatives The third derivatives of the function 2817 * specified by fe_function at the quadrature points of the current cell. 2818 * The third derivatives are computed in real space (as opposed to on the 2819 * unit cell). The object is assumed to already have the correct size. The 2820 * data type stored by this output vector must be what you get when you 2821 * multiply the third derivatives of shape function times the type used to 2822 * store the values of the unknowns $U_j$ of your finite element vector $U$ 2823 * (represented by the @p fe_function argument). 2824 * 2825 * @post <code>third_derivatives[q]</code> will contain the third 2826 * derivatives of the field described by fe_function at the $q$th quadrature 2827 * point. <code>third_derivatives[q][i][j][k]</code> represents the 2828 * $(i,j,k)$th component of the 3rd order tensor of third derivatives at 2829 * quadrature point $q$. 2830 * 2831 * @note The actual data type of the input vector may be either a 2832 * Vector<T>, BlockVector<T>, or one of the PETSc or Trilinos 2833 * vector wrapper classes. It represents a global vector of DoF values 2834 * associated with the DoFHandler object with which this FEValues object was 2835 * last initialized. 2836 * 2837 * @dealiiRequiresUpdateFlags{update_3rd_derivatives} 2838 */ 2839 template <class InputVector> 2840 void 2841 get_function_third_derivatives( 2842 const InputVector &fe_function, 2843 std::vector<Tensor<3, spacedim, typename InputVector::value_type>> 2844 &third_derivatives) const; 2845 2846 /** 2847 * This function does the same as the other 2848 * get_function_third_derivatives(), but applied to multi-component (vector- 2849 * valued) elements. The meaning of the arguments is as explained there. 2850 * 2851 * @post <code>third_derivatives[q]</code> is a vector of third derivatives 2852 * of the field described by fe_function at the $q$th quadrature point. The 2853 * size of the vector accessed by <code>third_derivatives[q]</code> equals 2854 * the number of components of the finite element, i.e. 2855 * <code>third_derivatives[q][c]</code> returns the third derivative of the 2856 * $c$th vector component at the $q$th quadrature point. Consequently, 2857 * <code>third_derivatives[q][c][i][j][k]</code> is the $(i,j,k)$th 2858 * component of the tensor of third derivatives of the $c$th vector 2859 * component of the vector field at quadrature point $q$ of the current 2860 * cell. 2861 * 2862 * @dealiiRequiresUpdateFlags{update_3rd_derivatives} 2863 */ 2864 template <class InputVector> 2865 void 2866 get_function_third_derivatives( 2867 const InputVector &fe_function, 2868 std::vector< 2869 std::vector<Tensor<3, spacedim, typename InputVector::value_type>>> 2870 & third_derivatives, 2871 bool quadrature_points_fastest = false) const; 2872 2873 /** 2874 * Access to the third derivatives of a function with more flexibility. See 2875 * get_function_values() with corresponding arguments. 2876 */ 2877 template <class InputVector> 2878 void 2879 get_function_third_derivatives( 2880 const InputVector & fe_function, 2881 const ArrayView<const types::global_dof_index> &indices, 2882 std::vector<Tensor<3, spacedim, typename InputVector::value_type>> 2883 &third_derivatives) const; 2884 2885 /** 2886 * Access to the third derivatives of a function with more flexibility. See 2887 * get_function_values() with corresponding arguments. 2888 * 2889 * @dealiiRequiresUpdateFlags{update_3rd_derivatives} 2890 */ 2891 template <class InputVector> 2892 void 2893 get_function_third_derivatives( 2894 const InputVector & fe_function, 2895 const ArrayView<const types::global_dof_index> &indices, 2896 ArrayView< 2897 std::vector<Tensor<3, spacedim, typename InputVector::value_type>>> 2898 third_derivatives, 2899 bool quadrature_points_fastest = false) const; 2900 //@} 2901 2902 /// @name Cell degrees of freedom 2903 //@{ 2904 2905 /** 2906 * Return an object that can be thought of as an array containing all 2907 * indices from zero (inclusive) to `dofs_per_cell` (exclusive). This allows 2908 * one to write code using range-based `for` loops of the following kind: 2909 * @code 2910 * FEValues<dim> fe_values (...); 2911 * FullMatrix<double> cell_matrix (...); 2912 * 2913 * for (auto &cell : dof_handler.active_cell_iterators()) 2914 * { 2915 * cell_matrix = 0; 2916 * fe_values.reinit(cell); 2917 * for (const auto q : fe_values.quadrature_point_indices()) 2918 * for (const auto i : fe_values.dof_indices()) 2919 * for (const auto j : fe_values.dof_indices()) 2920 * cell_matrix(i,j) += ...; // Do something for DoF indices (i,j) 2921 * // at quadrature point q 2922 * } 2923 * @endcode 2924 * Here, we are looping over all degrees of freedom on all cells, with 2925 * `i` and `j` taking on all valid indices for cell degrees of freedom, as 2926 * defined by the finite element passed to `fe_values`. 2927 */ 2928 std_cxx20::ranges::iota_view<unsigned int, unsigned int> 2929 dof_indices() const; 2930 2931 /** 2932 * Return an object that can be thought of as an array containing all 2933 * indices from @p start_dof_index (inclusive) to `dofs_per_cell` (exclusive). 2934 * This allows one to write code using range-based `for` loops of the 2935 * following kind: 2936 * @code 2937 * FEValues<dim> fe_values (...); 2938 * FullMatrix<double> cell_matrix (...); 2939 * 2940 * for (auto &cell : dof_handler.active_cell_iterators()) 2941 * { 2942 * cell_matrix = 0; 2943 * fe_values.reinit(cell); 2944 * for (const auto q : fe_values.quadrature_point_indices()) 2945 * for (const auto i : fe_values.dof_indices()) 2946 * for (const auto j : fe_values.dof_indices_starting_at(i)) 2947 * cell_matrix(i,j) += ...; // Do something for DoF indices (i,j) 2948 * // at quadrature point q 2949 * } 2950 * @endcode 2951 * Here, we are looping over all local degrees of freedom on all cells, with 2952 * `i` taking on all valid indices for cell degrees of freedom, as 2953 * defined by the finite element passed to `fe_values`, and `j` taking 2954 * on a specified subset of `i`'s range, starting at `i` itself and ending at 2955 * the number of cell degrees of freedom. In this way, we can construct the 2956 * upper half and the diagonal of a stiffness matrix contribution (assuming it 2957 * is symmetric, and that only one half of it needs to be computed), for 2958 * example. 2959 * 2960 * @note If the @p start_dof_index is equal to the number of DoFs in the cell, 2961 * then the returned index range is empty. 2962 */ 2963 std_cxx20::ranges::iota_view<unsigned int, unsigned int> 2964 dof_indices_starting_at(const unsigned int start_dof_index) const; 2965 2966 /** 2967 * Return an object that can be thought of as an array containing all 2968 * indices from zero (inclusive) to @p end_dof_index (inclusive). This allows 2969 * one to write code using range-based `for` loops of the following kind: 2970 * @code 2971 * FEValues<dim> fe_values (...); 2972 * FullMatrix<double> cell_matrix (...); 2973 * 2974 * for (auto &cell : dof_handler.active_cell_iterators()) 2975 * { 2976 * cell_matrix = 0; 2977 * fe_values.reinit(cell); 2978 * for (const auto q : fe_values.quadrature_point_indices()) 2979 * for (const auto i : fe_values.dof_indices()) 2980 * for (const auto j : fe_values.dof_indices_ending_at(i)) 2981 * cell_matrix(i,j) += ...; // Do something for DoF indices (i,j) 2982 * // at quadrature point q 2983 * } 2984 * @endcode 2985 * Here, we are looping over all local degrees of freedom on all cells, with 2986 * `i` taking on all valid indices for cell degrees of freedom, as 2987 * defined by the finite element passed to `fe_values`, and `j` taking 2988 * on a specified subset of `i`'s range, starting at zero and ending at 2989 * `i` itself. In this way, we can construct the lower half and the 2990 * diagonal of a stiffness matrix contribution (assuming it is symmetric, and 2991 * that only one half of it needs to be computed), for example. 2992 * 2993 * @note If the @p end_dof_index is equal to zero, then the returned index 2994 * range is empty. 2995 */ 2996 std_cxx20::ranges::iota_view<unsigned int, unsigned int> 2997 dof_indices_ending_at(const unsigned int end_dof_index) const; 2998 2999 //@} 3000 3001 /// @name Geometry of the cell 3002 //@{ 3003 3004 /** 3005 * Return an object that can be thought of as an array containing all 3006 * indices from zero to `n_quadrature_points`. This allows to write code 3007 * using range-based `for` loops of the following kind: 3008 * @code 3009 * FEValues<dim> fe_values (...); 3010 * 3011 * for (auto &cell : dof_handler.active_cell_iterators()) 3012 * { 3013 * fe_values.reinit(cell); 3014 * for (const auto q_point : fe_values.quadrature_point_indices()) 3015 * ... do something at the quadrature point ... 3016 * } 3017 * @endcode 3018 * Here, we are looping over all quadrature points on all cells, with 3019 * `q_point` taking on all valid indices for quadrature points, as defined 3020 * by the quadrature rule passed to `fe_values`. 3021 * 3022 * @see CPP11 3023 */ 3024 std_cxx20::ranges::iota_view<unsigned int, unsigned int> 3025 quadrature_point_indices() const; 3026 3027 /** 3028 * Position of the <tt>q</tt>th quadrature point in real space. 3029 * 3030 * @dealiiRequiresUpdateFlags{update_quadrature_points} 3031 */ 3032 const Point<spacedim> & 3033 quadrature_point(const unsigned int q) const; 3034 3035 /** 3036 * Return a reference to the vector of quadrature points in real space. 3037 * 3038 * @dealiiRequiresUpdateFlags{update_quadrature_points} 3039 */ 3040 const std::vector<Point<spacedim>> & 3041 get_quadrature_points() const; 3042 3043 /** 3044 * Mapped quadrature weight. If this object refers to a volume evaluation 3045 * (i.e. the derived class is of type FEValues), then this is the Jacobi 3046 * determinant times the weight of the *<tt>i</tt>th unit quadrature point. 3047 * 3048 * For surface evaluations (i.e. classes FEFaceValues or FESubfaceValues), 3049 * it is the mapped surface element times the weight of the quadrature 3050 * point. 3051 * 3052 * You can think of the quantity returned by this function as the volume or 3053 * surface element $dx, ds$ in the integral that we implement here by 3054 * quadrature. 3055 * 3056 * @dealiiRequiresUpdateFlags{update_JxW_values} 3057 */ 3058 double 3059 JxW(const unsigned int quadrature_point) const; 3060 3061 /** 3062 * Return a reference to the array holding the values returned by JxW(). 3063 */ 3064 const std::vector<double> & 3065 get_JxW_values() const; 3066 3067 /** 3068 * Return the Jacobian of the transformation at the specified quadrature 3069 * point, i.e. $J_{ij}=dx_i/d\hat x_j$ 3070 * 3071 * @dealiiRequiresUpdateFlags{update_jacobians} 3072 */ 3073 const DerivativeForm<1, dim, spacedim> & 3074 jacobian(const unsigned int quadrature_point) const; 3075 3076 /** 3077 * Return a reference to the array holding the values returned by 3078 * jacobian(). 3079 * 3080 * @dealiiRequiresUpdateFlags{update_jacobians} 3081 */ 3082 const std::vector<DerivativeForm<1, dim, spacedim>> & 3083 get_jacobians() const; 3084 3085 /** 3086 * Return the second derivative of the transformation from unit to real 3087 * cell, i.e. the first derivative of the Jacobian, at the specified 3088 * quadrature point, i.e. $G_{ijk}=dJ_{jk}/d\hat x_i$. 3089 * 3090 * @dealiiRequiresUpdateFlags{update_jacobian_grads} 3091 */ 3092 const DerivativeForm<2, dim, spacedim> & 3093 jacobian_grad(const unsigned int quadrature_point) const; 3094 3095 /** 3096 * Return a reference to the array holding the values returned by 3097 * jacobian_grads(). 3098 * 3099 * @dealiiRequiresUpdateFlags{update_jacobian_grads} 3100 */ 3101 const std::vector<DerivativeForm<2, dim, spacedim>> & 3102 get_jacobian_grads() const; 3103 3104 /** 3105 * Return the second derivative of the transformation from unit to real 3106 * cell, i.e. the first derivative of the Jacobian, at the specified 3107 * quadrature point, pushed forward to the real cell coordinates, i.e. 3108 * $G_{ijk}=dJ_{iJ}/d\hat x_K (J_{jJ})^{-1} (J_{kK})^{-1}$. 3109 * 3110 * @dealiiRequiresUpdateFlags{update_jacobian_pushed_forward_grads} 3111 */ 3112 const Tensor<3, spacedim> & 3113 jacobian_pushed_forward_grad(const unsigned int quadrature_point) const; 3114 3115 /** 3116 * Return a reference to the array holding the values returned by 3117 * jacobian_pushed_forward_grads(). 3118 * 3119 * @dealiiRequiresUpdateFlags{update_jacobian_pushed_forward_grads} 3120 */ 3121 const std::vector<Tensor<3, spacedim>> & 3122 get_jacobian_pushed_forward_grads() const; 3123 3124 /** 3125 * Return the third derivative of the transformation from unit to real cell, 3126 * i.e. the second derivative of the Jacobian, at the specified quadrature 3127 * point, i.e. $G_{ijkl}=\frac{d^2J_{ij}}{d\hat x_k d\hat x_l}$. 3128 * 3129 * @dealiiRequiresUpdateFlags{update_jacobian_2nd_derivatives} 3130 */ 3131 const DerivativeForm<3, dim, spacedim> & 3132 jacobian_2nd_derivative(const unsigned int quadrature_point) const; 3133 3134 /** 3135 * Return a reference to the array holding the values returned by 3136 * jacobian_2nd_derivatives(). 3137 * 3138 * @dealiiRequiresUpdateFlags{update_jacobian_2nd_derivatives} 3139 */ 3140 const std::vector<DerivativeForm<3, dim, spacedim>> & 3141 get_jacobian_2nd_derivatives() const; 3142 3143 /** 3144 * Return the third derivative of the transformation from unit to real cell, 3145 * i.e. the second derivative of the Jacobian, at the specified quadrature 3146 * point, pushed forward to the real cell coordinates, i.e. 3147 * $G_{ijkl}=\frac{d^2J_{iJ}}{d\hat x_K d\hat x_L} (J_{jJ})^{-1} 3148 * (J_{kK})^{-1}(J_{lL})^{-1}$. 3149 * 3150 * @dealiiRequiresUpdateFlags{update_jacobian_pushed_forward_2nd_derivatives} 3151 */ 3152 const Tensor<4, spacedim> & 3153 jacobian_pushed_forward_2nd_derivative( 3154 const unsigned int quadrature_point) const; 3155 3156 /** 3157 * Return a reference to the array holding the values returned by 3158 * jacobian_pushed_forward_2nd_derivatives(). 3159 * 3160 * @dealiiRequiresUpdateFlags{update_jacobian_pushed_forward_2nd_derivatives} 3161 */ 3162 const std::vector<Tensor<4, spacedim>> & 3163 get_jacobian_pushed_forward_2nd_derivatives() const; 3164 3165 /** 3166 * Return the fourth derivative of the transformation from unit to real 3167 * cell, i.e. the third derivative of the Jacobian, at the specified 3168 * quadrature point, i.e. $G_{ijklm}=\frac{d^2J_{ij}}{d\hat x_k d\hat x_l 3169 * d\hat x_m}$. 3170 * 3171 * @dealiiRequiresUpdateFlags{update_jacobian_3rd_derivatives} 3172 */ 3173 const DerivativeForm<4, dim, spacedim> & 3174 jacobian_3rd_derivative(const unsigned int quadrature_point) const; 3175 3176 /** 3177 * Return a reference to the array holding the values returned by 3178 * jacobian_3rd_derivatives(). 3179 * 3180 * @dealiiRequiresUpdateFlags{update_jacobian_3rd_derivatives} 3181 */ 3182 const std::vector<DerivativeForm<4, dim, spacedim>> & 3183 get_jacobian_3rd_derivatives() const; 3184 3185 /** 3186 * Return the fourth derivative of the transformation from unit to real 3187 * cell, i.e. the third derivative of the Jacobian, at the specified 3188 * quadrature point, pushed forward to the real cell coordinates, i.e. 3189 * $G_{ijklm}=\frac{d^3J_{iJ}}{d\hat x_K d\hat x_L d\hat x_M} (J_{jJ})^{-1} 3190 * (J_{kK})^{-1} (J_{lL})^{-1} (J_{mM})^{-1}$. 3191 * 3192 * @dealiiRequiresUpdateFlags{update_jacobian_pushed_forward_3rd_derivatives} 3193 */ 3194 const Tensor<5, spacedim> & 3195 jacobian_pushed_forward_3rd_derivative( 3196 const unsigned int quadrature_point) const; 3197 3198 /** 3199 * Return a reference to the array holding the values returned by 3200 * jacobian_pushed_forward_3rd_derivatives(). 3201 * 3202 * @dealiiRequiresUpdateFlags{update_jacobian_pushed_forward_2nd_derivatives} 3203 */ 3204 const std::vector<Tensor<5, spacedim>> & 3205 get_jacobian_pushed_forward_3rd_derivatives() const; 3206 3207 /** 3208 * Return the inverse Jacobian of the transformation at the specified 3209 * quadrature point, i.e. $J_{ij}=d\hat x_i/dx_j$ 3210 * 3211 * @dealiiRequiresUpdateFlags{update_inverse_jacobians} 3212 */ 3213 const DerivativeForm<1, spacedim, dim> & 3214 inverse_jacobian(const unsigned int quadrature_point) const; 3215 3216 /** 3217 * Return a reference to the array holding the values returned by 3218 * inverse_jacobian(). 3219 * 3220 * @dealiiRequiresUpdateFlags{update_inverse_jacobians} 3221 */ 3222 const std::vector<DerivativeForm<1, spacedim, dim>> & 3223 get_inverse_jacobians() const; 3224 3225 /** 3226 * Return the normal vector at a quadrature point. If you call this 3227 * function for a face (i.e., when using a FEFaceValues or FESubfaceValues 3228 * object), then this function returns the outward normal vector to 3229 * the cell at the <tt>i</tt>th quadrature point of the face. 3230 * 3231 * In contrast, if you call this function for a cell of codimension one 3232 * (i.e., when using a `FEValues<dim,spacedim>` object with 3233 * `spacedim>dim`), then this function returns the normal vector to the 3234 * cell -- in other words, an approximation to the normal vector to the 3235 * manifold in which the triangulation is embedded. There are of 3236 * course two normal directions to a manifold in that case, and this 3237 * function returns the "up" direction as induced by the numbering of the 3238 * vertices. 3239 * 3240 * The length of the vector is normalized to one. 3241 * 3242 * @dealiiRequiresUpdateFlags{update_normal_vectors} 3243 */ 3244 const Tensor<1, spacedim> & 3245 normal_vector(const unsigned int i) const; 3246 3247 /** 3248 * Return the normal vectors at all quadrature points represented by 3249 * this object. See the normal_vector() function for what the normal 3250 * vectors represent. 3251 * 3252 * @dealiiRequiresUpdateFlags{update_normal_vectors} 3253 */ 3254 const std::vector<Tensor<1, spacedim>> & 3255 get_normal_vectors() const; 3256 3257 //@} 3258 3259 /// @name Extractors Methods to extract individual components 3260 //@{ 3261 3262 /** 3263 * Create a view of the current FEValues object that represents a particular 3264 * scalar component of the possibly vector-valued finite element. The 3265 * concept of views is explained in the documentation of the namespace 3266 * FEValuesViews and in particular in the 3267 * @ref vector_valued 3268 * module. 3269 */ 3270 const FEValuesViews::Scalar<dim, spacedim> & 3271 operator[](const FEValuesExtractors::Scalar &scalar) const; 3272 3273 /** 3274 * Create a view of the current FEValues object that represents a set of 3275 * <code>dim</code> scalar components (i.e. a vector) of the vector-valued 3276 * finite element. The concept of views is explained in the documentation of 3277 * the namespace FEValuesViews and in particular in the 3278 * @ref vector_valued 3279 * module. 3280 */ 3281 const FEValuesViews::Vector<dim, spacedim> & 3282 operator[](const FEValuesExtractors::Vector &vector) const; 3283 3284 /** 3285 * Create a view of the current FEValues object that represents a set of 3286 * <code>(dim*dim + dim)/2</code> scalar components (i.e. a symmetric 2nd 3287 * order tensor) of the vector-valued finite element. The concept of views 3288 * is explained in the documentation of the namespace FEValuesViews and in 3289 * particular in the 3290 * @ref vector_valued 3291 * module. 3292 */ 3293 const FEValuesViews::SymmetricTensor<2, dim, spacedim> & 3294 operator[](const FEValuesExtractors::SymmetricTensor<2> &tensor) const; 3295 3296 3297 /** 3298 * Create a view of the current FEValues object that represents a set of 3299 * <code>(dim*dim)</code> scalar components (i.e. a 2nd order tensor) of the 3300 * vector-valued finite element. The concept of views is explained in the 3301 * documentation of the namespace FEValuesViews and in particular in the 3302 * @ref vector_valued 3303 * module. 3304 */ 3305 const FEValuesViews::Tensor<2, dim, spacedim> & 3306 operator[](const FEValuesExtractors::Tensor<2> &tensor) const; 3307 3308 //@} 3309 3310 /// @name Access to the raw data 3311 //@{ 3312 3313 /** 3314 * Constant reference to the selected mapping object. 3315 */ 3316 const Mapping<dim, spacedim> & 3317 get_mapping() const; 3318 3319 /** 3320 * Constant reference to the selected finite element object. 3321 */ 3322 const FiniteElement<dim, spacedim> & 3323 get_fe() const; 3324 3325 /** 3326 * Return the update flags set for this object. 3327 */ 3328 UpdateFlags 3329 get_update_flags() const; 3330 3331 /** 3332 * Return a triangulation iterator to the current cell. 3333 */ 3334 const typename Triangulation<dim, spacedim>::cell_iterator 3335 get_cell() const; 3336 3337 /** 3338 * Return the relation of the current cell to the previous cell. This allows 3339 * re-use of some cell data (like local matrices for equations with constant 3340 * coefficients) if the result is <tt>CellSimilarity::translation</tt>. 3341 */ 3342 CellSimilarity::Similarity 3343 get_cell_similarity() const; 3344 3345 /** 3346 * Determine an estimate for the memory consumption (in bytes) of this 3347 * object. 3348 */ 3349 std::size_t 3350 memory_consumption() const; 3351 //@} 3352 3353 3354 /** 3355 * This exception is thrown if FEValuesBase is asked to return the value of 3356 * a field which was not required by the UpdateFlags for this FEValuesBase. 3357 * 3358 * @ingroup Exceptions 3359 */ 3360 DeclException1( 3361 ExcAccessToUninitializedField, 3362 std::string, 3363 << "You are requesting information from an FEValues/FEFaceValues/FESubfaceValues " 3364 << "object for which this kind of information has not been computed. What " 3365 << "information these objects compute is determined by the update_* flags you " 3366 << "pass to the constructor. Here, the operation you are attempting requires " 3367 << "the <" << arg1 3368 << "> flag to be set, but it was apparently not specified " 3369 << "upon construction."); 3370 3371 /** 3372 * Mismatch between the FEValues FiniteElement and 3373 * cell->get_dof_handler().get_fe() 3374 * 3375 * @ingroup Exceptions 3376 */ 3377 DeclExceptionMsg( 3378 ExcFEDontMatch, 3379 "The FiniteElement you provided to FEValues and the FiniteElement that belongs " 3380 "to the DoFHandler that provided the cell iterator do not match."); 3381 /** 3382 * A given shape function is not primitive, but it needs to be. 3383 * 3384 * @ingroup Exceptions 3385 */ 3386 DeclException1(ExcShapeFunctionNotPrimitive, 3387 int, 3388 << "The shape function with index " << arg1 3389 << " is not primitive, i.e. it is vector-valued and " 3390 << "has more than one non-zero vector component. This " 3391 << "function cannot be called for these shape functions. " 3392 << "Maybe you want to use the same function with the " 3393 << "_component suffix?"); 3394 3395 /** 3396 * The given FiniteElement is not a primitive element, see 3397 * FiniteElement::is_primitive(). 3398 * 3399 * @ingroup Exceptions 3400 */ 3401 DeclExceptionMsg( 3402 ExcFENotPrimitive, 3403 "The given FiniteElement is not a primitive element but the requested operation " 3404 "only works for those. See FiniteElement::is_primitive() for more information."); 3405 3406 protected: 3407 /** 3408 * Objects of the FEValues class need to store an iterator 3409 * to the present cell in order to be able to extract the values of the 3410 * degrees of freedom on this cell in the get_function_values() and assorted 3411 * functions. On the other hand, this class should also work for different 3412 * iterators, as long as they have the same interface to extract the DoF 3413 * values (i.e., for example, they need to have a @p 3414 * get_interpolated_dof_values function). 3415 * 3416 * This calls for a common base class of iterator classes, and making the 3417 * functions we need here @p virtual. On the other hand, this is the only 3418 * place in the library where we need this, and introducing a base class of 3419 * iterators and making a function virtual penalizes <em>all</em> users of 3420 * the iterators, which are basically intended as very fast accessor 3421 * functions. So we do not want to do this. Rather, what we do here is 3422 * making the functions we need virtual only for use with <em>this 3423 * class</em>. The idea is the following: have a common base class which 3424 * declares some pure virtual functions, and for each possible iterator 3425 * type, we have a derived class which stores the iterator to the cell and 3426 * implements these functions. Since the iterator classes have the same 3427 * interface, we can make the derived classes a template, templatized on the 3428 * iterator type. 3429 * 3430 * This way, the use of virtual functions is restricted to only this class, 3431 * and other users of iterators do not have to bear the negative effects. 3432 * 3433 * @note This class is an example of the 3434 * <a href="https://www.artima.com/cppsource/type_erasure.html">type 3435 * erasure</a> design pattern. 3436 */ 3437 class CellIteratorBase; 3438 3439 /** 3440 * Forward declaration of classes derived from CellIteratorBase. Their 3441 * definition and implementation is given in the .cc file. 3442 */ 3443 template <typename CI> 3444 class CellIterator; 3445 class TriaCellIterator; 3446 3447 /** 3448 * Store the cell selected last time the reinit() function was called. This 3449 * is necessary for the <tt>get_function_*</tt> functions as well as the 3450 * functions of same name in the extractor classes. 3451 */ 3452 std::unique_ptr<const CellIteratorBase> present_cell; 3453 3454 /** 3455 * A signal connection we use to ensure we get informed whenever the 3456 * triangulation changes by refinement. We need to know about that because 3457 * it invalidates all cell iterators and, as part of that, the 3458 * 'present_cell' iterator we keep around between subsequent calls to 3459 * reinit() in order to compute the cell similarity. 3460 */ 3461 boost::signals2::connection tria_listener_refinement; 3462 3463 /** 3464 * A signal connection we use to ensure we get informed whenever the 3465 * triangulation changes by mesh transformations. We need to know about that 3466 * because it invalidates all cell iterators and, as part of that, the 3467 * 'present_cell' iterator we keep around between subsequent calls to 3468 * reinit() in order to compute the cell similarity. 3469 */ 3470 boost::signals2::connection tria_listener_mesh_transform; 3471 3472 /** 3473 * A function that is connected to the triangulation in order to reset the 3474 * stored 'present_cell' iterator to an invalid one whenever the 3475 * triangulation is changed and the iterator consequently becomes invalid. 3476 */ 3477 void 3478 invalidate_present_cell(); 3479 3480 /** 3481 * This function is called by the various reinit() functions in derived 3482 * classes. Given the cell indicated by the argument, test whether we have 3483 * to throw away the previously stored present_cell argument because it 3484 * would require us to compare cells from different triangulations. In 3485 * checking all this, also make sure that we have tria_listener connected to 3486 * the triangulation to which we will set present_cell right after calling 3487 * this function. 3488 */ 3489 void 3490 maybe_invalidate_previous_present_cell( 3491 const typename Triangulation<dim, spacedim>::cell_iterator &cell); 3492 3493 /** 3494 * A pointer to the mapping object associated with this FEValues object. 3495 */ 3496 const SmartPointer<const Mapping<dim, spacedim>, FEValuesBase<dim, spacedim>> 3497 mapping; 3498 3499 /** 3500 * A pointer to the internal data object of mapping, obtained from 3501 * Mapping::get_data(), Mapping::get_face_data(), or 3502 * Mapping::get_subface_data(). 3503 */ 3504 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> 3505 mapping_data; 3506 3507 /** 3508 * An object into which the Mapping::fill_fe_values() and similar functions 3509 * place their output. 3510 */ 3511 dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> 3512 mapping_output; 3513 3514 3515 /** 3516 * A pointer to the finite element object associated with this FEValues 3517 * object. 3518 */ 3519 const SmartPointer<const FiniteElement<dim, spacedim>, 3520 FEValuesBase<dim, spacedim>> 3521 fe; 3522 3523 /** 3524 * A pointer to the internal data object of finite element, obtained from 3525 * FiniteElement::get_data(), Mapping::get_face_data(), or 3526 * FiniteElement::get_subface_data(). 3527 */ 3528 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase> 3529 fe_data; 3530 3531 /** 3532 * An object into which the FiniteElement::fill_fe_values() and similar 3533 * functions place their output. 3534 */ 3535 dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, 3536 spacedim> 3537 finite_element_output; 3538 3539 3540 /** 3541 * Original update flags handed to the constructor of FEValues. 3542 */ 3543 UpdateFlags update_flags; 3544 3545 /** 3546 * Initialize some update flags. Called from the @p initialize functions of 3547 * derived classes, which are in turn called from their constructors. 3548 * 3549 * Basically, this function finds out using the finite element and mapping 3550 * object already stored which flags need to be set to compute everything 3551 * the user wants, as expressed through the flags passed as argument. 3552 */ 3553 UpdateFlags 3554 compute_update_flags(const UpdateFlags update_flags) const; 3555 3556 /** 3557 * An enum variable that can store different states of the current cell in 3558 * comparison to the previously visited cell. If wanted, additional states 3559 * can be checked here and used in one of the methods used during reinit. 3560 */ 3561 CellSimilarity::Similarity cell_similarity; 3562 3563 /** 3564 * A function that checks whether the new cell is similar to the one 3565 * previously used. Then, a significant amount of the data can be reused, 3566 * e.g. the derivatives of the basis functions in real space, shape_grad. 3567 */ 3568 void 3569 check_cell_similarity( 3570 const typename Triangulation<dim, spacedim>::cell_iterator &cell); 3571 3572 private: 3573 /** 3574 * A cache for all possible FEValuesViews objects. 3575 */ 3576 dealii::internal::FEValuesViews::Cache<dim, spacedim> fe_values_views_cache; 3577 3578 // Make the view classes friends of this class, since they access internal 3579 // data. 3580 template <int, int> 3581 friend class FEValuesViews::Scalar; 3582 template <int, int> 3583 friend class FEValuesViews::Vector; 3584 template <int, int, int> 3585 friend class FEValuesViews::SymmetricTensor; 3586 template <int, int, int> 3587 friend class FEValuesViews::Tensor; 3588 }; 3589 3590 3591 3592 /** 3593 * Finite element evaluated in quadrature points of a cell. 3594 * 3595 * This function implements the initialization routines for FEValuesBase, if 3596 * values in quadrature points of a cell are needed. For further documentation 3597 * see this class. 3598 * 3599 * @ingroup feaccess 3600 */ 3601 template <int dim, int spacedim = dim> 3602 class FEValues : public FEValuesBase<dim, spacedim> 3603 { 3604 public: 3605 /** 3606 * Dimension of the object over which we integrate. For the present class, 3607 * this is equal to <code>dim</code>. 3608 */ 3609 static const unsigned int integral_dimension = dim; 3610 3611 /** 3612 * Constructor. Gets cell independent data from mapping and finite element 3613 * objects, matching the quadrature rule and update flags. 3614 */ 3615 FEValues(const Mapping<dim, spacedim> & mapping, 3616 const FiniteElement<dim, spacedim> &fe, 3617 const Quadrature<dim> & quadrature, 3618 const UpdateFlags update_flags); 3619 3620 /** 3621 * Constructor. This constructor is equivalent to the other one except that 3622 * it makes the object use a $Q_1$ mapping (i.e., an object of type 3623 * MappingQGeneric(1)) implicitly. 3624 */ 3625 FEValues(const FiniteElement<dim, spacedim> &fe, 3626 const Quadrature<dim> & quadrature, 3627 const UpdateFlags update_flags); 3628 3629 /** 3630 * Reinitialize the gradients, Jacobi determinants, etc for the given cell 3631 * of type "iterator into a DoFHandler object", and the finite element 3632 * associated with this object. It is assumed that the finite element used 3633 * by the given cell is also the one used by this FEValues object. 3634 */ 3635 template <bool level_dof_access> 3636 void 3637 reinit( 3638 const TriaIterator<DoFCellAccessor<dim, spacedim, level_dof_access>> &cell); 3639 3640 /** 3641 * Reinitialize the gradients, Jacobi determinants, etc for the given cell 3642 * of type "iterator into a Triangulation object", and the given finite 3643 * element. Since iterators into triangulation alone only convey information 3644 * about the geometry of a cell, but not about degrees of freedom possibly 3645 * associated with this cell, you will not be able to call some functions of 3646 * this class if they need information about degrees of freedom. These 3647 * functions are, above all, the 3648 * <tt>get_function_value/gradients/hessians/laplacians/third_derivatives</tt> 3649 * functions. If you want to call these functions, you have to call the @p 3650 * reinit variants that take iterators into DoFHandler or other DoF handler 3651 * type objects. 3652 */ 3653 void 3654 reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell); 3655 3656 /** 3657 * Return a reference to the copy of the quadrature formula stored by this 3658 * object. 3659 */ 3660 const Quadrature<dim> & 3661 get_quadrature() const; 3662 3663 /** 3664 * Determine an estimate for the memory consumption (in bytes) of this 3665 * object. 3666 */ 3667 std::size_t 3668 memory_consumption() const; 3669 3670 /** 3671 * Return a reference to this very object. 3672 * 3673 * Though it seems that it is not very useful, this function is there to 3674 * provide capability to the hp::FEValues class, in which case it provides 3675 * the FEValues object for the present cell (remember that for hp finite 3676 * elements, the actual FE object used may change from cell to cell, so we 3677 * also need different FEValues objects for different cells; once you 3678 * reinitialize the hp::FEValues object for a specific cell, it retrieves 3679 * the FEValues object for the FE on that cell and returns it through a 3680 * function of the same name as this one; this function here therefore only 3681 * provides the same interface so that one can templatize on FEValues and 3682 * hp::FEValues). 3683 */ 3684 const FEValues<dim, spacedim> & 3685 get_present_fe_values() const; 3686 3687 private: 3688 /** 3689 * Store a copy of the quadrature formula here. 3690 */ 3691 const Quadrature<dim> quadrature; 3692 3693 /** 3694 * Do work common to the two constructors. 3695 */ 3696 void 3697 initialize(const UpdateFlags update_flags); 3698 3699 /** 3700 * The reinit() functions do only that part of the work that requires 3701 * knowledge of the type of iterator. After setting present_cell(), they 3702 * pass on to this function, which does the real work, and which is 3703 * independent of the actual type of the cell iterator. 3704 */ 3705 void 3706 do_reinit(); 3707 }; 3708 3709 3710 /** 3711 * Extend the interface of FEValuesBase to values that only make sense when 3712 * evaluating something on the surface of a cell. All the data that is 3713 * available in the interior of cells is also available here. 3714 * 3715 * See FEValuesBase 3716 * 3717 * @ingroup feaccess 3718 */ 3719 template <int dim, int spacedim = dim> 3720 class FEFaceValuesBase : public FEValuesBase<dim, spacedim> 3721 { 3722 public: 3723 /** 3724 * Dimension of the object over which we integrate. For the present class, 3725 * this is equal to <code>dim-1</code>. 3726 */ 3727 static const unsigned int integral_dimension = dim - 1; 3728 3729 /** 3730 * Constructor. Call the constructor of the base class and set up the arrays 3731 * of this class with the right sizes. Actually filling these arrays is a 3732 * duty of the derived class's constructors. 3733 * 3734 * @p n_faces_or_subfaces is the number of faces or subfaces that this 3735 * object is to store. The actual number depends on the derived class, for 3736 * FEFaceValues it is <tt>2*dim</tt>, while for the FESubfaceValues class it 3737 * is <tt>2*dim*(1<<(dim-1))</tt>, i.e. the number of faces times the number 3738 * of subfaces per face. 3739 */ 3740 FEFaceValuesBase(const unsigned int n_q_points, 3741 const unsigned int dofs_per_cell, 3742 const UpdateFlags update_flags, 3743 const Mapping<dim, spacedim> & mapping, 3744 const FiniteElement<dim, spacedim> &fe, 3745 const Quadrature<dim - 1> & quadrature); 3746 3747 /** 3748 * Boundary form of the transformation of the cell at the <tt>i</tt>th 3749 * quadrature point. See 3750 * @ref GlossBoundaryForm. 3751 * 3752 * @dealiiRequiresUpdateFlags{update_boundary_forms} 3753 */ 3754 const Tensor<1, spacedim> & 3755 boundary_form(const unsigned int i) const; 3756 3757 /** 3758 * Return the list of outward normal vectors times the Jacobian of the 3759 * surface mapping. 3760 * 3761 * @dealiiRequiresUpdateFlags{update_boundary_forms} 3762 */ 3763 const std::vector<Tensor<1, spacedim>> & 3764 get_boundary_forms() const; 3765 3766 /** 3767 * Return the index of the face selected the last time the reinit() function 3768 * was called. 3769 */ 3770 unsigned int 3771 get_face_index() const; 3772 3773 /** 3774 * Return a reference to the copy of the quadrature formula stored by this 3775 * object. 3776 */ 3777 const Quadrature<dim - 1> & 3778 get_quadrature() const; 3779 3780 /** 3781 * Determine an estimate for the memory consumption (in bytes) of this 3782 * object. 3783 */ 3784 std::size_t 3785 memory_consumption() const; 3786 3787 protected: 3788 /** 3789 * Index of the face selected the last time the reinit() function was 3790 * called. 3791 */ 3792 unsigned int present_face_index; 3793 3794 /** 3795 * Store a copy of the quadrature formula here. 3796 */ 3797 const Quadrature<dim - 1> quadrature; 3798 }; 3799 3800 3801 3802 /** 3803 * Finite element evaluated in quadrature points on a face. 3804 * 3805 * This class adds the functionality of FEFaceValuesBase to FEValues; see 3806 * there for more documentation. 3807 * 3808 * Since finite element functions and their derivatives may be discontinuous 3809 * at cell boundaries, there is no restriction of this function to a mesh 3810 * face. But, there are limits of these values approaching the face from 3811 * either of the neighboring cells. 3812 * 3813 * @ingroup feaccess 3814 */ 3815 template <int dim, int spacedim = dim> 3816 class FEFaceValues : public FEFaceValuesBase<dim, spacedim> 3817 { 3818 public: 3819 /** 3820 * Dimension in which this object operates. 3821 */ 3822 3823 static const unsigned int dimension = dim; 3824 3825 static const unsigned int space_dimension = spacedim; 3826 3827 /** 3828 * Dimension of the object over which we integrate. For the present class, 3829 * this is equal to <code>dim-1</code>. 3830 */ 3831 static const unsigned int integral_dimension = dim - 1; 3832 3833 /** 3834 * Constructor. Gets cell independent data from mapping and finite element 3835 * objects, matching the quadrature rule and update flags. 3836 */ 3837 FEFaceValues(const Mapping<dim, spacedim> & mapping, 3838 const FiniteElement<dim, spacedim> &fe, 3839 const Quadrature<dim - 1> & quadrature, 3840 const UpdateFlags update_flags); 3841 3842 /** 3843 * Constructor. This constructor is equivalent to the other one except that 3844 * it makes the object use a $Q_1$ mapping (i.e., an object of type 3845 * MappingQGeneric(1)) implicitly. 3846 */ 3847 FEFaceValues(const FiniteElement<dim, spacedim> &fe, 3848 const Quadrature<dim - 1> & quadrature, 3849 const UpdateFlags update_flags); 3850 3851 /** 3852 * Reinitialize the gradients, Jacobi determinants, etc for the face with 3853 * number @p face_no of @p cell and the given finite element. 3854 */ 3855 template <bool level_dof_access> 3856 void 3857 reinit( 3858 const TriaIterator<DoFCellAccessor<dim, spacedim, level_dof_access>> &cell, 3859 const unsigned int face_no); 3860 3861 /** 3862 * Reinitialize the gradients, Jacobi determinants, etc for face @p face 3863 * and cell @p cell. 3864 * 3865 * @note @p face must be one of @p cell's face iterators. 3866 */ 3867 template <bool level_dof_access> 3868 void 3869 reinit( 3870 const TriaIterator<DoFCellAccessor<dim, spacedim, level_dof_access>> &cell, 3871 const typename Triangulation<dim, spacedim>::face_iterator & face); 3872 3873 /** 3874 * Reinitialize the gradients, Jacobi determinants, etc for the given face 3875 * on a given cell of type "iterator into a Triangulation object", and the 3876 * given finite element. Since iterators into a triangulation alone only 3877 * convey information about the geometry of a cell, but not about degrees of 3878 * freedom possibly associated with this cell, you will not be able to call 3879 * some functions of this class if they need information about degrees of 3880 * freedom. These functions are, above all, the 3881 * <tt>get_function_value/gradients/hessians/third_derivatives</tt> 3882 * functions. If you want to call these functions, you have to call the @p 3883 * reinit variants that take iterators into DoFHandler or other DoF handler 3884 * type objects. 3885 */ 3886 void 3887 reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell, 3888 const unsigned int face_no); 3889 3890 /* 3891 * Reinitialize the gradients, Jacobi determinants, etc for the given face 3892 * on a given cell of type "iterator into a Triangulation object", and the 3893 * given finite element. Since iterators into a triangulation alone only 3894 * convey information about the geometry of a cell, but not about degrees of 3895 * freedom possibly associated with this cell, you will not be able to call 3896 * some functions of this class if they need information about degrees of 3897 * freedom. These functions are, above all, the 3898 * <tt>get_function_value/gradients/hessians/third_derivatives</tt> 3899 * functions. If you want to call these functions, you have to call the @p 3900 * reinit variants that take iterators into DoFHandler or other DoF handler 3901 * type objects. 3902 * 3903 * @note @p face must be one of @p cell's face iterators. 3904 */ 3905 void 3906 reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell, 3907 const typename Triangulation<dim, spacedim>::face_iterator &face); 3908 3909 /** 3910 * Return a reference to this very object. 3911 * 3912 * Though it seems that it is not very useful, this function is there to 3913 * provide capability to the hp::FEValues class, in which case it provides 3914 * the FEValues object for the present cell (remember that for hp finite 3915 * elements, the actual FE object used may change from cell to cell, so we 3916 * also need different FEValues objects for different cells; once you 3917 * reinitialize the hp::FEValues object for a specific cell, it retrieves 3918 * the FEValues object for the FE on that cell and returns it through a 3919 * function of the same name as this one; this function here therefore only 3920 * provides the same interface so that one can templatize on FEValues and 3921 * hp::FEValues). 3922 */ 3923 const FEFaceValues<dim, spacedim> & 3924 get_present_fe_values() const; 3925 3926 private: 3927 /** 3928 * Do work common to the two constructors. 3929 */ 3930 void 3931 initialize(const UpdateFlags update_flags); 3932 3933 /** 3934 * The reinit() functions do only that part of the work that requires 3935 * knowledge of the type of iterator. After setting present_cell(), they 3936 * pass on to this function, which does the real work, and which is 3937 * independent of the actual type of the cell iterator. 3938 */ 3939 void 3940 do_reinit(const unsigned int face_no); 3941 }; 3942 3943 3944 /** 3945 * Finite element evaluated in quadrature points on a face. 3946 * 3947 * This class adds the functionality of FEFaceValuesBase to FEValues; see 3948 * there for more documentation. 3949 * 3950 * This class is used for faces lying on a refinement edge. In this case, the 3951 * neighboring cell is refined. To be able to compute differences between 3952 * interior and exterior function values, the refinement of the neighboring 3953 * cell must be simulated on this cell. This is achieved by applying a 3954 * quadrature rule that simulates the refinement. The resulting data fields 3955 * are split up to reflect the refinement structure of the neighbor: a subface 3956 * number corresponds to the number of the child of the neighboring face. 3957 * 3958 * @ingroup feaccess 3959 */ 3960 template <int dim, int spacedim = dim> 3961 class FESubfaceValues : public FEFaceValuesBase<dim, spacedim> 3962 { 3963 public: 3964 /** 3965 * Dimension in which this object operates. 3966 */ 3967 static const unsigned int dimension = dim; 3968 3969 /** 3970 * Dimension of the space in which this object operates. 3971 */ 3972 static const unsigned int space_dimension = spacedim; 3973 3974 /** 3975 * Dimension of the object over which we integrate. For the present class, 3976 * this is equal to <code>dim-1</code>. 3977 */ 3978 static const unsigned int integral_dimension = dim - 1; 3979 3980 /** 3981 * Constructor. Gets cell independent data from mapping and finite element 3982 * objects, matching the quadrature rule and update flags. 3983 */ 3984 FESubfaceValues(const Mapping<dim, spacedim> & mapping, 3985 const FiniteElement<dim, spacedim> &fe, 3986 const Quadrature<dim - 1> & face_quadrature, 3987 const UpdateFlags update_flags); 3988 3989 /** 3990 * Constructor. This constructor is equivalent to the other one except that 3991 * it makes the object use a $Q_1$ mapping (i.e., an object of type 3992 * MappingQGeneric(1)) implicitly. 3993 */ 3994 FESubfaceValues(const FiniteElement<dim, spacedim> &fe, 3995 const Quadrature<dim - 1> & face_quadrature, 3996 const UpdateFlags update_flags); 3997 3998 /** 3999 * Reinitialize the gradients, Jacobi determinants, etc for the given cell 4000 * of type "iterator into a DoFHandler object", and the finite element 4001 * associated with this object. It is assumed that the finite element used 4002 * by the given cell is also the one used by this FESubfaceValues object. 4003 */ 4004 template <bool level_dof_access> 4005 void 4006 reinit( 4007 const TriaIterator<DoFCellAccessor<dim, spacedim, level_dof_access>> &cell, 4008 const unsigned int face_no, 4009 const unsigned int subface_no); 4010 4011 /** 4012 * Alternative reinitialization function that takes, as arguments, iterators 4013 * to the face and subface instead of their numbers. 4014 */ 4015 template <bool level_dof_access> 4016 void 4017 reinit( 4018 const TriaIterator<DoFCellAccessor<dim, spacedim, level_dof_access>> &cell, 4019 const typename Triangulation<dim, spacedim>::face_iterator & face, 4020 const typename Triangulation<dim, spacedim>::face_iterator &subface); 4021 4022 /** 4023 * Reinitialize the gradients, Jacobi determinants, etc for the given 4024 * subface on a given cell of type "iterator into a Triangulation object", and 4025 * the given finite element. Since iterators into a triangulation alone only 4026 * convey information about the geometry of a cell, but not about degrees of 4027 * freedom possibly associated with this cell, you will not be able to call 4028 * some functions of this class if they need information about degrees of 4029 * freedom. These functions are, above all, the 4030 * <tt>get_function_value/gradients/hessians/third_derivatives</tt> 4031 * functions. If you want to call these functions, you have to call the @p 4032 * reinit variants that take iterators into DoFHandler or other DoF handler 4033 * type objects. 4034 */ 4035 void 4036 reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell, 4037 const unsigned int face_no, 4038 const unsigned int subface_no); 4039 4040 /** 4041 * Reinitialize the gradients, Jacobi determinants, etc for the given 4042 * subface on a given cell of type "iterator into a Triangulation object", and 4043 * the given finite element. Since iterators into a triangulation alone only 4044 * convey information about the geometry of a cell, but not about degrees of 4045 * freedom possibly associated with this cell, you will not be able to call 4046 * some functions of this class if they need information about degrees of 4047 * freedom. These functions are, above all, the 4048 * <tt>get_function_value/gradients/hessians/third_derivatives</tt> 4049 * functions. If you want to call these functions, you have to call the @p 4050 * reinit variants that take iterators into DoFHandler or other DoF handler 4051 * type objects. 4052 * 4053 * This does the same thing as the previous function but takes iterators 4054 * instead of numbers as arguments. 4055 * 4056 * @note @p face and @p subface must correspond to a face (and a subface of 4057 * that face) of @p cell. 4058 */ 4059 void 4060 reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell, 4061 const typename Triangulation<dim, spacedim>::face_iterator &face, 4062 const typename Triangulation<dim, spacedim>::face_iterator &subface); 4063 4064 /** 4065 * Return a reference to this very object. 4066 * 4067 * Though it seems that it is not very useful, this function is there to 4068 * provide capability to the hp::FEValues class, in which case it provides 4069 * the FEValues object for the present cell (remember that for hp finite 4070 * elements, the actual FE object used may change from cell to cell, so we 4071 * also need different FEValues objects for different cells; once you 4072 * reinitialize the hp::FEValues object for a specific cell, it retrieves 4073 * the FEValues object for the FE on that cell and returns it through a 4074 * function of the same name as this one; this function here therefore only 4075 * provides the same interface so that one can templatize on FEValues and 4076 * hp::FEValues). 4077 */ 4078 const FESubfaceValues<dim, spacedim> & 4079 get_present_fe_values() const; 4080 4081 /** 4082 * @todo Document this 4083 * 4084 * @ingroup Exceptions 4085 */ 4086 DeclException0(ExcReinitCalledWithBoundaryFace); 4087 4088 /** 4089 * @todo Document this 4090 * 4091 * @ingroup Exceptions 4092 */ 4093 DeclException0(ExcFaceHasNoSubfaces); 4094 4095 private: 4096 /** 4097 * Do work common to the two constructors. 4098 */ 4099 void 4100 initialize(const UpdateFlags update_flags); 4101 4102 /** 4103 * The reinit() functions do only that part of the work that requires 4104 * knowledge of the type of iterator. After setting present_cell(), they 4105 * pass on to this function, which does the real work, and which is 4106 * independent of the actual type of the cell iterator. 4107 */ 4108 void 4109 do_reinit(const unsigned int face_no, const unsigned int subface_no); 4110 }; 4111 4112 4113 #ifndef DOXYGEN 4114 4115 4116 /*------------------------ Inline functions: namespace FEValuesViews --------*/ 4117 4118 namespace FEValuesViews 4119 { 4120 template <int dim, int spacedim> 4121 inline typename Scalar<dim, spacedim>::value_type 4122 Scalar<dim, spacedim>::value(const unsigned int shape_function, 4123 const unsigned int q_point) const 4124 { 4125 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4126 Assert( 4127 fe_values->update_flags & update_values, 4128 ((typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4129 "update_values")))); 4130 4131 // an adaptation of the FEValuesBase::shape_value_component function 4132 // except that here we know the component as fixed and we have 4133 // pre-computed and cached a bunch of information. See the comments there. 4134 if (shape_function_data[shape_function].is_nonzero_shape_function_component) 4135 return fe_values->finite_element_output.shape_values( 4136 shape_function_data[shape_function].row_index, q_point); 4137 else 4138 return 0; 4139 } 4140 4141 4142 4143 template <int dim, int spacedim> 4144 inline typename Scalar<dim, spacedim>::gradient_type 4145 Scalar<dim, spacedim>::gradient(const unsigned int shape_function, 4146 const unsigned int q_point) const 4147 { 4148 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4149 Assert(fe_values->update_flags & update_gradients, 4150 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4151 "update_gradients"))); 4152 4153 // an adaptation of the FEValuesBase::shape_grad_component 4154 // function except that here we know the component as fixed and we have 4155 // pre-computed and cached a bunch of information. See the comments there. 4156 if (shape_function_data[shape_function].is_nonzero_shape_function_component) 4157 return fe_values->finite_element_output 4158 .shape_gradients[shape_function_data[shape_function].row_index] 4159 [q_point]; 4160 else 4161 return gradient_type(); 4162 } 4163 4164 4165 4166 template <int dim, int spacedim> 4167 inline typename Scalar<dim, spacedim>::hessian_type 4168 Scalar<dim, spacedim>::hessian(const unsigned int shape_function, 4169 const unsigned int q_point) const 4170 { 4171 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4172 Assert(fe_values->update_flags & update_hessians, 4173 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4174 "update_hessians"))); 4175 4176 // an adaptation of the FEValuesBase::shape_hessian_component 4177 // function except that here we know the component as fixed and we have 4178 // pre-computed and cached a bunch of information. See the comments there. 4179 if (shape_function_data[shape_function].is_nonzero_shape_function_component) 4180 return fe_values->finite_element_output 4181 .shape_hessians[shape_function_data[shape_function].row_index][q_point]; 4182 else 4183 return hessian_type(); 4184 } 4185 4186 4187 4188 template <int dim, int spacedim> 4189 inline typename Scalar<dim, spacedim>::third_derivative_type 4190 Scalar<dim, spacedim>::third_derivative(const unsigned int shape_function, 4191 const unsigned int q_point) const 4192 { 4193 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4194 Assert(fe_values->update_flags & update_3rd_derivatives, 4195 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4196 "update_3rd_derivatives"))); 4197 4198 // an adaptation of the FEValuesBase::shape_3rdderivative_component 4199 // function except that here we know the component as fixed and we have 4200 // pre-computed and cached a bunch of information. See the comments there. 4201 if (shape_function_data[shape_function].is_nonzero_shape_function_component) 4202 return fe_values->finite_element_output 4203 .shape_3rd_derivatives[shape_function_data[shape_function].row_index] 4204 [q_point]; 4205 else 4206 return third_derivative_type(); 4207 } 4208 4209 4210 4211 template <int dim, int spacedim> 4212 inline typename Vector<dim, spacedim>::value_type 4213 Vector<dim, spacedim>::value(const unsigned int shape_function, 4214 const unsigned int q_point) const 4215 { 4216 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4217 Assert(fe_values->update_flags & update_values, 4218 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4219 "update_values"))); 4220 4221 // same as for the scalar case except that we have one more index 4222 const int snc = 4223 shape_function_data[shape_function].single_nonzero_component; 4224 if (snc == -2) 4225 return value_type(); 4226 else if (snc != -1) 4227 { 4228 value_type return_value; 4229 return_value[shape_function_data[shape_function] 4230 .single_nonzero_component_index] = 4231 fe_values->finite_element_output.shape_values(snc, q_point); 4232 return return_value; 4233 } 4234 else 4235 { 4236 value_type return_value; 4237 for (unsigned int d = 0; d < dim; ++d) 4238 if (shape_function_data[shape_function] 4239 .is_nonzero_shape_function_component[d]) 4240 return_value[d] = fe_values->finite_element_output.shape_values( 4241 shape_function_data[shape_function].row_index[d], q_point); 4242 4243 return return_value; 4244 } 4245 } 4246 4247 4248 4249 template <int dim, int spacedim> 4250 inline typename Vector<dim, spacedim>::gradient_type 4251 Vector<dim, spacedim>::gradient(const unsigned int shape_function, 4252 const unsigned int q_point) const 4253 { 4254 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4255 Assert(fe_values->update_flags & update_gradients, 4256 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4257 "update_gradients"))); 4258 4259 // same as for the scalar case except that we have one more index 4260 const int snc = 4261 shape_function_data[shape_function].single_nonzero_component; 4262 if (snc == -2) 4263 return gradient_type(); 4264 else if (snc != -1) 4265 { 4266 gradient_type return_value; 4267 return_value[shape_function_data[shape_function] 4268 .single_nonzero_component_index] = 4269 fe_values->finite_element_output.shape_gradients[snc][q_point]; 4270 return return_value; 4271 } 4272 else 4273 { 4274 gradient_type return_value; 4275 for (unsigned int d = 0; d < dim; ++d) 4276 if (shape_function_data[shape_function] 4277 .is_nonzero_shape_function_component[d]) 4278 return_value[d] = 4279 fe_values->finite_element_output.shape_gradients 4280 [shape_function_data[shape_function].row_index[d]][q_point]; 4281 4282 return return_value; 4283 } 4284 } 4285 4286 4287 4288 template <int dim, int spacedim> 4289 inline typename Vector<dim, spacedim>::divergence_type 4290 Vector<dim, spacedim>::divergence(const unsigned int shape_function, 4291 const unsigned int q_point) const 4292 { 4293 // this function works like in the case above 4294 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4295 Assert(fe_values->update_flags & update_gradients, 4296 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4297 "update_gradients"))); 4298 4299 // same as for the scalar case except that we have one more index 4300 const int snc = 4301 shape_function_data[shape_function].single_nonzero_component; 4302 if (snc == -2) 4303 return divergence_type(); 4304 else if (snc != -1) 4305 return fe_values->finite_element_output 4306 .shape_gradients[snc][q_point][shape_function_data[shape_function] 4307 .single_nonzero_component_index]; 4308 else 4309 { 4310 divergence_type return_value = 0; 4311 for (unsigned int d = 0; d < dim; ++d) 4312 if (shape_function_data[shape_function] 4313 .is_nonzero_shape_function_component[d]) 4314 return_value += 4315 fe_values->finite_element_output.shape_gradients 4316 [shape_function_data[shape_function].row_index[d]][q_point][d]; 4317 4318 return return_value; 4319 } 4320 } 4321 4322 4323 4324 template <int dim, int spacedim> 4325 inline typename Vector<dim, spacedim>::curl_type 4326 Vector<dim, spacedim>::curl(const unsigned int shape_function, 4327 const unsigned int q_point) const 4328 { 4329 // this function works like in the case above 4330 4331 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4332 Assert(fe_values->update_flags & update_gradients, 4333 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4334 "update_gradients"))); 4335 // same as for the scalar case except that we have one more index 4336 const int snc = 4337 shape_function_data[shape_function].single_nonzero_component; 4338 4339 if (snc == -2) 4340 return curl_type(); 4341 4342 else 4343 switch (dim) 4344 { 4345 case 1: 4346 { 4347 Assert(false, 4348 ExcMessage( 4349 "Computing the curl in 1d is not a useful operation")); 4350 return curl_type(); 4351 } 4352 4353 case 2: 4354 { 4355 if (snc != -1) 4356 { 4357 curl_type return_value; 4358 4359 // the single nonzero component can only be zero or one in 2d 4360 if (shape_function_data[shape_function] 4361 .single_nonzero_component_index == 0) 4362 return_value[0] = 4363 -1.0 * fe_values->finite_element_output 4364 .shape_gradients[snc][q_point][1]; 4365 else 4366 return_value[0] = fe_values->finite_element_output 4367 .shape_gradients[snc][q_point][0]; 4368 4369 return return_value; 4370 } 4371 4372 else 4373 { 4374 curl_type return_value; 4375 4376 return_value[0] = 0.0; 4377 4378 if (shape_function_data[shape_function] 4379 .is_nonzero_shape_function_component[0]) 4380 return_value[0] -= 4381 fe_values->finite_element_output 4382 .shape_gradients[shape_function_data[shape_function] 4383 .row_index[0]][q_point][1]; 4384 4385 if (shape_function_data[shape_function] 4386 .is_nonzero_shape_function_component[1]) 4387 return_value[0] += 4388 fe_values->finite_element_output 4389 .shape_gradients[shape_function_data[shape_function] 4390 .row_index[1]][q_point][0]; 4391 4392 return return_value; 4393 } 4394 } 4395 4396 case 3: 4397 { 4398 if (snc != -1) 4399 { 4400 curl_type return_value; 4401 4402 switch (shape_function_data[shape_function] 4403 .single_nonzero_component_index) 4404 { 4405 case 0: 4406 { 4407 return_value[0] = 0; 4408 return_value[1] = fe_values->finite_element_output 4409 .shape_gradients[snc][q_point][2]; 4410 return_value[2] = 4411 -1.0 * fe_values->finite_element_output 4412 .shape_gradients[snc][q_point][1]; 4413 return return_value; 4414 } 4415 4416 case 1: 4417 { 4418 return_value[0] = 4419 -1.0 * fe_values->finite_element_output 4420 .shape_gradients[snc][q_point][2]; 4421 return_value[1] = 0; 4422 return_value[2] = fe_values->finite_element_output 4423 .shape_gradients[snc][q_point][0]; 4424 return return_value; 4425 } 4426 4427 default: 4428 { 4429 return_value[0] = fe_values->finite_element_output 4430 .shape_gradients[snc][q_point][1]; 4431 return_value[1] = 4432 -1.0 * fe_values->finite_element_output 4433 .shape_gradients[snc][q_point][0]; 4434 return_value[2] = 0; 4435 return return_value; 4436 } 4437 } 4438 } 4439 4440 else 4441 { 4442 curl_type return_value; 4443 4444 for (unsigned int i = 0; i < dim; ++i) 4445 return_value[i] = 0.0; 4446 4447 if (shape_function_data[shape_function] 4448 .is_nonzero_shape_function_component[0]) 4449 { 4450 return_value[1] += 4451 fe_values->finite_element_output 4452 .shape_gradients[shape_function_data[shape_function] 4453 .row_index[0]][q_point][2]; 4454 return_value[2] -= 4455 fe_values->finite_element_output 4456 .shape_gradients[shape_function_data[shape_function] 4457 .row_index[0]][q_point][1]; 4458 } 4459 4460 if (shape_function_data[shape_function] 4461 .is_nonzero_shape_function_component[1]) 4462 { 4463 return_value[0] -= 4464 fe_values->finite_element_output 4465 .shape_gradients[shape_function_data[shape_function] 4466 .row_index[1]][q_point][2]; 4467 return_value[2] += 4468 fe_values->finite_element_output 4469 .shape_gradients[shape_function_data[shape_function] 4470 .row_index[1]][q_point][0]; 4471 } 4472 4473 if (shape_function_data[shape_function] 4474 .is_nonzero_shape_function_component[2]) 4475 { 4476 return_value[0] += 4477 fe_values->finite_element_output 4478 .shape_gradients[shape_function_data[shape_function] 4479 .row_index[2]][q_point][1]; 4480 return_value[1] -= 4481 fe_values->finite_element_output 4482 .shape_gradients[shape_function_data[shape_function] 4483 .row_index[2]][q_point][0]; 4484 } 4485 4486 return return_value; 4487 } 4488 } 4489 } 4490 // should not end up here 4491 Assert(false, ExcInternalError()); 4492 return curl_type(); 4493 } 4494 4495 4496 4497 template <int dim, int spacedim> 4498 inline typename Vector<dim, spacedim>::hessian_type 4499 Vector<dim, spacedim>::hessian(const unsigned int shape_function, 4500 const unsigned int q_point) const 4501 { 4502 // this function works like in the case above 4503 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4504 Assert(fe_values->update_flags & update_hessians, 4505 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4506 "update_hessians"))); 4507 4508 // same as for the scalar case except that we have one more index 4509 const int snc = 4510 shape_function_data[shape_function].single_nonzero_component; 4511 if (snc == -2) 4512 return hessian_type(); 4513 else if (snc != -1) 4514 { 4515 hessian_type return_value; 4516 return_value[shape_function_data[shape_function] 4517 .single_nonzero_component_index] = 4518 fe_values->finite_element_output.shape_hessians[snc][q_point]; 4519 return return_value; 4520 } 4521 else 4522 { 4523 hessian_type return_value; 4524 for (unsigned int d = 0; d < dim; ++d) 4525 if (shape_function_data[shape_function] 4526 .is_nonzero_shape_function_component[d]) 4527 return_value[d] = 4528 fe_values->finite_element_output.shape_hessians 4529 [shape_function_data[shape_function].row_index[d]][q_point]; 4530 4531 return return_value; 4532 } 4533 } 4534 4535 4536 4537 template <int dim, int spacedim> 4538 inline typename Vector<dim, spacedim>::third_derivative_type 4539 Vector<dim, spacedim>::third_derivative(const unsigned int shape_function, 4540 const unsigned int q_point) const 4541 { 4542 // this function works like in the case above 4543 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4544 Assert(fe_values->update_flags & update_3rd_derivatives, 4545 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4546 "update_3rd_derivatives"))); 4547 4548 // same as for the scalar case except that we have one more index 4549 const int snc = 4550 shape_function_data[shape_function].single_nonzero_component; 4551 if (snc == -2) 4552 return third_derivative_type(); 4553 else if (snc != -1) 4554 { 4555 third_derivative_type return_value; 4556 return_value[shape_function_data[shape_function] 4557 .single_nonzero_component_index] = 4558 fe_values->finite_element_output.shape_3rd_derivatives[snc][q_point]; 4559 return return_value; 4560 } 4561 else 4562 { 4563 third_derivative_type return_value; 4564 for (unsigned int d = 0; d < dim; ++d) 4565 if (shape_function_data[shape_function] 4566 .is_nonzero_shape_function_component[d]) 4567 return_value[d] = 4568 fe_values->finite_element_output.shape_3rd_derivatives 4569 [shape_function_data[shape_function].row_index[d]][q_point]; 4570 4571 return return_value; 4572 } 4573 } 4574 4575 4576 4577 namespace internal 4578 { 4579 /** 4580 * Return the symmetrized version of a tensor whose n'th row equals the 4581 * second argument, with all other rows equal to zero. 4582 */ 4583 inline dealii::SymmetricTensor<2, 1> 4584 symmetrize_single_row(const unsigned int n, const dealii::Tensor<1, 1> &t) 4585 { 4586 AssertIndexRange(n, 1); 4587 (void)n; 4588 4589 return {{t[0]}}; 4590 } 4591 4592 4593 4594 inline dealii::SymmetricTensor<2, 2> 4595 symmetrize_single_row(const unsigned int n, const dealii::Tensor<1, 2> &t) 4596 { 4597 switch (n) 4598 { 4599 case 0: 4600 { 4601 return {{t[0], 0, t[1] / 2}}; 4602 } 4603 case 1: 4604 { 4605 return {{0, t[1], t[0] / 2}}; 4606 } 4607 default: 4608 { 4609 AssertIndexRange(n, 2); 4610 return {}; 4611 } 4612 } 4613 } 4614 4615 4616 4617 inline dealii::SymmetricTensor<2, 3> 4618 symmetrize_single_row(const unsigned int n, const dealii::Tensor<1, 3> &t) 4619 { 4620 switch (n) 4621 { 4622 case 0: 4623 { 4624 return {{t[0], 0, 0, t[1] / 2, t[2] / 2, 0}}; 4625 } 4626 case 1: 4627 { 4628 return {{0, t[1], 0, t[0] / 2, 0, t[2] / 2}}; 4629 } 4630 case 2: 4631 { 4632 return {{0, 0, t[2], 0, t[0] / 2, t[1] / 2}}; 4633 } 4634 default: 4635 { 4636 AssertIndexRange(n, 3); 4637 return {}; 4638 } 4639 } 4640 } 4641 } // namespace internal 4642 4643 4644 4645 template <int dim, int spacedim> 4646 inline typename Vector<dim, spacedim>::symmetric_gradient_type 4647 Vector<dim, spacedim>::symmetric_gradient(const unsigned int shape_function, 4648 const unsigned int q_point) const 4649 { 4650 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4651 Assert(fe_values->update_flags & update_gradients, 4652 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4653 "update_gradients"))); 4654 4655 // same as for the scalar case except that we have one more index 4656 const int snc = 4657 shape_function_data[shape_function].single_nonzero_component; 4658 if (snc == -2) 4659 return symmetric_gradient_type(); 4660 else if (snc != -1) 4661 return internal::symmetrize_single_row( 4662 shape_function_data[shape_function].single_nonzero_component_index, 4663 fe_values->finite_element_output.shape_gradients[snc][q_point]); 4664 else 4665 { 4666 gradient_type return_value; 4667 for (unsigned int d = 0; d < dim; ++d) 4668 if (shape_function_data[shape_function] 4669 .is_nonzero_shape_function_component[d]) 4670 return_value[d] = 4671 fe_values->finite_element_output.shape_gradients 4672 [shape_function_data[shape_function].row_index[d]][q_point]; 4673 4674 return symmetrize(return_value); 4675 } 4676 } 4677 4678 4679 4680 template <int dim, int spacedim> 4681 inline typename SymmetricTensor<2, dim, spacedim>::value_type 4682 SymmetricTensor<2, dim, spacedim>::value(const unsigned int shape_function, 4683 const unsigned int q_point) const 4684 { 4685 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4686 Assert(fe_values->update_flags & update_values, 4687 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4688 "update_values"))); 4689 4690 // similar to the vector case where we have more then one index and we need 4691 // to convert between unrolled and component indexing for tensors 4692 const int snc = 4693 shape_function_data[shape_function].single_nonzero_component; 4694 4695 if (snc == -2) 4696 { 4697 // shape function is zero for the selected components 4698 return value_type(); 4699 } 4700 else if (snc != -1) 4701 { 4702 value_type return_value; 4703 const unsigned int comp = 4704 shape_function_data[shape_function].single_nonzero_component_index; 4705 return_value[value_type::unrolled_to_component_indices(comp)] = 4706 fe_values->finite_element_output.shape_values(snc, q_point); 4707 return return_value; 4708 } 4709 else 4710 { 4711 value_type return_value; 4712 for (unsigned int d = 0; d < value_type::n_independent_components; ++d) 4713 if (shape_function_data[shape_function] 4714 .is_nonzero_shape_function_component[d]) 4715 return_value[value_type::unrolled_to_component_indices(d)] = 4716 fe_values->finite_element_output.shape_values( 4717 shape_function_data[shape_function].row_index[d], q_point); 4718 return return_value; 4719 } 4720 } 4721 4722 4723 4724 template <int dim, int spacedim> 4725 inline typename SymmetricTensor<2, dim, spacedim>::divergence_type 4726 SymmetricTensor<2, dim, spacedim>::divergence( 4727 const unsigned int shape_function, 4728 const unsigned int q_point) const 4729 { 4730 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4731 Assert(fe_values->update_flags & update_gradients, 4732 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4733 "update_gradients"))); 4734 4735 const int snc = 4736 shape_function_data[shape_function].single_nonzero_component; 4737 4738 if (snc == -2) 4739 { 4740 // shape function is zero for the selected components 4741 return divergence_type(); 4742 } 4743 else if (snc != -1) 4744 { 4745 // we have a single non-zero component when the symmetric tensor is 4746 // represented in unrolled form. this implies we potentially have 4747 // two non-zero components when represented in component form! we 4748 // will only have one non-zero entry if the non-zero component lies on 4749 // the diagonal of the tensor. 4750 // 4751 // the divergence of a second-order tensor is a first order tensor. 4752 // 4753 // assume the second-order tensor is A with components A_{ij}. then 4754 // A_{ij} = A_{ji} and there is only one (if diagonal) or two non-zero 4755 // entries in the tensorial representation. define the 4756 // divergence as: 4757 // b_i \dealcoloneq \dfrac{\partial phi_{ij}}{\partial x_j}. 4758 // (which is incidentally also 4759 // b_j \dealcoloneq \dfrac{\partial phi_{ij}}{\partial x_i}). 4760 // In both cases, a sum is implied. 4761 // 4762 // Now, we know the nonzero component in unrolled form: it is indicated 4763 // by 'snc'. we can figure out which tensor components belong to this: 4764 const unsigned int comp = 4765 shape_function_data[shape_function].single_nonzero_component_index; 4766 const unsigned int ii = 4767 value_type::unrolled_to_component_indices(comp)[0]; 4768 const unsigned int jj = 4769 value_type::unrolled_to_component_indices(comp)[1]; 4770 4771 // given the form of the divergence above, if ii=jj there is only a 4772 // single nonzero component of the full tensor and the gradient 4773 // equals 4774 // b_ii \dealcoloneq \dfrac{\partial phi_{ii,ii}}{\partial x_ii}. 4775 // all other entries of 'b' are zero 4776 // 4777 // on the other hand, if ii!=jj, then there are two nonzero entries in 4778 // the full tensor and 4779 // b_ii \dealcoloneq \dfrac{\partial phi_{ii,jj}}{\partial x_ii}. 4780 // b_jj \dealcoloneq \dfrac{\partial phi_{ii,jj}}{\partial x_jj}. 4781 // again, all other entries of 'b' are zero 4782 const dealii::Tensor<1, spacedim> &phi_grad = 4783 fe_values->finite_element_output.shape_gradients[snc][q_point]; 4784 4785 divergence_type return_value; 4786 return_value[ii] = phi_grad[jj]; 4787 4788 if (ii != jj) 4789 return_value[jj] = phi_grad[ii]; 4790 4791 return return_value; 4792 } 4793 else 4794 { 4795 Assert(false, ExcNotImplemented()); 4796 divergence_type return_value; 4797 return return_value; 4798 } 4799 } 4800 4801 4802 4803 template <int dim, int spacedim> 4804 inline typename Tensor<2, dim, spacedim>::value_type 4805 Tensor<2, dim, spacedim>::value(const unsigned int shape_function, 4806 const unsigned int q_point) const 4807 { 4808 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4809 Assert(fe_values->update_flags & update_values, 4810 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4811 "update_values"))); 4812 4813 // similar to the vector case where we have more then one index and we need 4814 // to convert between unrolled and component indexing for tensors 4815 const int snc = 4816 shape_function_data[shape_function].single_nonzero_component; 4817 4818 if (snc == -2) 4819 { 4820 // shape function is zero for the selected components 4821 return value_type(); 4822 } 4823 else if (snc != -1) 4824 { 4825 value_type return_value; 4826 const unsigned int comp = 4827 shape_function_data[shape_function].single_nonzero_component_index; 4828 const TableIndices<2> indices = 4829 dealii::Tensor<2, spacedim>::unrolled_to_component_indices(comp); 4830 return_value[indices] = 4831 fe_values->finite_element_output.shape_values(snc, q_point); 4832 return return_value; 4833 } 4834 else 4835 { 4836 value_type return_value; 4837 for (unsigned int d = 0; d < dim * dim; ++d) 4838 if (shape_function_data[shape_function] 4839 .is_nonzero_shape_function_component[d]) 4840 { 4841 const TableIndices<2> indices = 4842 dealii::Tensor<2, spacedim>::unrolled_to_component_indices(d); 4843 return_value[indices] = 4844 fe_values->finite_element_output.shape_values( 4845 shape_function_data[shape_function].row_index[d], q_point); 4846 } 4847 return return_value; 4848 } 4849 } 4850 4851 4852 4853 template <int dim, int spacedim> 4854 inline typename Tensor<2, dim, spacedim>::divergence_type 4855 Tensor<2, dim, spacedim>::divergence(const unsigned int shape_function, 4856 const unsigned int q_point) const 4857 { 4858 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4859 Assert(fe_values->update_flags & update_gradients, 4860 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4861 "update_gradients"))); 4862 4863 const int snc = 4864 shape_function_data[shape_function].single_nonzero_component; 4865 4866 if (snc == -2) 4867 { 4868 // shape function is zero for the selected components 4869 return divergence_type(); 4870 } 4871 else if (snc != -1) 4872 { 4873 // we have a single non-zero component when the tensor is 4874 // represented in unrolled form. 4875 // 4876 // the divergence of a second-order tensor is a first order tensor. 4877 // 4878 // assume the second-order tensor is A with components A_{ij}, 4879 // then divergence is d_i := \frac{\partial A_{ij}}{\partial x_j} 4880 // 4881 // Now, we know the nonzero component in unrolled form: it is indicated 4882 // by 'snc'. we can figure out which tensor components belong to this: 4883 const unsigned int comp = 4884 shape_function_data[shape_function].single_nonzero_component_index; 4885 const TableIndices<2> indices = 4886 dealii::Tensor<2, spacedim>::unrolled_to_component_indices(comp); 4887 const unsigned int ii = indices[0]; 4888 const unsigned int jj = indices[1]; 4889 4890 const dealii::Tensor<1, spacedim> &phi_grad = 4891 fe_values->finite_element_output.shape_gradients[snc][q_point]; 4892 4893 divergence_type return_value; 4894 // note that we contract \nabla from the right 4895 return_value[ii] = phi_grad[jj]; 4896 4897 return return_value; 4898 } 4899 else 4900 { 4901 Assert(false, ExcNotImplemented()); 4902 divergence_type return_value; 4903 return return_value; 4904 } 4905 } 4906 4907 4908 4909 template <int dim, int spacedim> 4910 inline typename Tensor<2, dim, spacedim>::gradient_type 4911 Tensor<2, dim, spacedim>::gradient(const unsigned int shape_function, 4912 const unsigned int q_point) const 4913 { 4914 AssertIndexRange(shape_function, fe_values->fe->n_dofs_per_cell()); 4915 Assert(fe_values->update_flags & update_gradients, 4916 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 4917 "update_gradients"))); 4918 4919 const int snc = 4920 shape_function_data[shape_function].single_nonzero_component; 4921 4922 if (snc == -2) 4923 { 4924 // shape function is zero for the selected components 4925 return gradient_type(); 4926 } 4927 else if (snc != -1) 4928 { 4929 // we have a single non-zero component when the tensor is 4930 // represented in unrolled form. 4931 // 4932 // the gradient of a second-order tensor is a third order tensor. 4933 // 4934 // assume the second-order tensor is A with components A_{ij}, 4935 // then gradient is B_{ijk} := \frac{\partial A_{ij}}{\partial x_k} 4936 // 4937 // Now, we know the nonzero component in unrolled form: it is indicated 4938 // by 'snc'. we can figure out which tensor components belong to this: 4939 const unsigned int comp = 4940 shape_function_data[shape_function].single_nonzero_component_index; 4941 const TableIndices<2> indices = 4942 dealii::Tensor<2, spacedim>::unrolled_to_component_indices(comp); 4943 const unsigned int ii = indices[0]; 4944 const unsigned int jj = indices[1]; 4945 4946 const dealii::Tensor<1, spacedim> &phi_grad = 4947 fe_values->finite_element_output.shape_gradients[snc][q_point]; 4948 4949 gradient_type return_value; 4950 return_value[ii][jj] = phi_grad; 4951 4952 return return_value; 4953 } 4954 else 4955 { 4956 Assert(false, ExcNotImplemented()); 4957 gradient_type return_value; 4958 return return_value; 4959 } 4960 } 4961 4962 } // namespace FEValuesViews 4963 4964 4965 4966 /*---------------------- Inline functions: FEValuesBase ---------------------*/ 4967 4968 4969 4970 template <int dim, int spacedim> 4971 inline const FEValuesViews::Scalar<dim, spacedim> &FEValuesBase<dim, spacedim>:: 4972 operator[](const FEValuesExtractors::Scalar &scalar) const 4973 { 4974 AssertIndexRange(scalar.component, fe_values_views_cache.scalars.size()); 4975 4976 return fe_values_views_cache.scalars[scalar.component]; 4977 } 4978 4979 4980 4981 template <int dim, int spacedim> 4982 inline const FEValuesViews::Vector<dim, spacedim> &FEValuesBase<dim, spacedim>:: 4983 operator[](const FEValuesExtractors::Vector &vector) const 4984 { 4985 AssertIndexRange(vector.first_vector_component, 4986 fe_values_views_cache.vectors.size()); 4987 4988 return fe_values_views_cache.vectors[vector.first_vector_component]; 4989 } 4990 4991 4992 4993 template <int dim, int spacedim> 4994 inline const FEValuesViews::SymmetricTensor<2, dim, spacedim> & 4995 FEValuesBase<dim, spacedim>:: 4996 operator[](const FEValuesExtractors::SymmetricTensor<2> &tensor) const 4997 { 4998 Assert( 4999 tensor.first_tensor_component < 5000 fe_values_views_cache.symmetric_second_order_tensors.size(), 5001 ExcIndexRange(tensor.first_tensor_component, 5002 0, 5003 fe_values_views_cache.symmetric_second_order_tensors.size())); 5004 5005 return fe_values_views_cache 5006 .symmetric_second_order_tensors[tensor.first_tensor_component]; 5007 } 5008 5009 5010 5011 template <int dim, int spacedim> 5012 inline const FEValuesViews::Tensor<2, dim, spacedim> & 5013 FEValuesBase<dim, spacedim>:: 5014 operator[](const FEValuesExtractors::Tensor<2> &tensor) const 5015 { 5016 AssertIndexRange(tensor.first_tensor_component, 5017 fe_values_views_cache.second_order_tensors.size()); 5018 5019 return fe_values_views_cache 5020 .second_order_tensors[tensor.first_tensor_component]; 5021 } 5022 5023 5024 5025 template <int dim, int spacedim> 5026 inline const double & 5027 FEValuesBase<dim, spacedim>::shape_value(const unsigned int i, 5028 const unsigned int j) const 5029 { 5030 AssertIndexRange(i, fe->n_dofs_per_cell()); 5031 Assert(this->update_flags & update_values, 5032 ExcAccessToUninitializedField("update_values")); 5033 Assert(fe->is_primitive(i), ExcShapeFunctionNotPrimitive(i)); 5034 Assert(present_cell.get() != nullptr, 5035 ExcMessage("FEValues object is not reinit'ed to any cell")); 5036 // if the entire FE is primitive, 5037 // then we can take a short-cut: 5038 if (fe->is_primitive()) 5039 return this->finite_element_output.shape_values(i, j); 5040 else 5041 { 5042 // otherwise, use the mapping 5043 // between shape function 5044 // numbers and rows. note that 5045 // by the assertions above, we 5046 // know that this particular 5047 // shape function is primitive, 5048 // so we can call 5049 // system_to_component_index 5050 const unsigned int row = 5051 this->finite_element_output 5052 .shape_function_to_row_table[i * fe->n_components() + 5053 fe->system_to_component_index(i).first]; 5054 return this->finite_element_output.shape_values(row, j); 5055 } 5056 } 5057 5058 5059 5060 template <int dim, int spacedim> 5061 inline double 5062 FEValuesBase<dim, spacedim>::shape_value_component( 5063 const unsigned int i, 5064 const unsigned int j, 5065 const unsigned int component) const 5066 { 5067 AssertIndexRange(i, fe->n_dofs_per_cell()); 5068 Assert(this->update_flags & update_values, 5069 ExcAccessToUninitializedField("update_values")); 5070 AssertIndexRange(component, fe->n_components()); 5071 Assert(present_cell.get() != nullptr, 5072 ExcMessage("FEValues object is not reinit'ed to any cell")); 5073 5074 // check whether the shape function 5075 // is non-zero at all within 5076 // this component: 5077 if (fe->get_nonzero_components(i)[component] == false) 5078 return 0; 5079 5080 // look up the right row in the 5081 // table and take the data from 5082 // there 5083 const unsigned int row = 5084 this->finite_element_output 5085 .shape_function_to_row_table[i * fe->n_components() + component]; 5086 return this->finite_element_output.shape_values(row, j); 5087 } 5088 5089 5090 5091 template <int dim, int spacedim> 5092 inline const Tensor<1, spacedim> & 5093 FEValuesBase<dim, spacedim>::shape_grad(const unsigned int i, 5094 const unsigned int j) const 5095 { 5096 AssertIndexRange(i, fe->n_dofs_per_cell()); 5097 Assert(this->update_flags & update_gradients, 5098 ExcAccessToUninitializedField("update_gradients")); 5099 Assert(fe->is_primitive(i), ExcShapeFunctionNotPrimitive(i)); 5100 Assert(present_cell.get() != nullptr, 5101 ExcMessage("FEValues object is not reinit'ed to any cell")); 5102 // if the entire FE is primitive, 5103 // then we can take a short-cut: 5104 if (fe->is_primitive()) 5105 return this->finite_element_output.shape_gradients[i][j]; 5106 else 5107 { 5108 // otherwise, use the mapping 5109 // between shape function 5110 // numbers and rows. note that 5111 // by the assertions above, we 5112 // know that this particular 5113 // shape function is primitive, 5114 // so we can call 5115 // system_to_component_index 5116 const unsigned int row = 5117 this->finite_element_output 5118 .shape_function_to_row_table[i * fe->n_components() + 5119 fe->system_to_component_index(i).first]; 5120 return this->finite_element_output.shape_gradients[row][j]; 5121 } 5122 } 5123 5124 5125 5126 template <int dim, int spacedim> 5127 inline Tensor<1, spacedim> 5128 FEValuesBase<dim, spacedim>::shape_grad_component( 5129 const unsigned int i, 5130 const unsigned int j, 5131 const unsigned int component) const 5132 { 5133 AssertIndexRange(i, fe->n_dofs_per_cell()); 5134 Assert(this->update_flags & update_gradients, 5135 ExcAccessToUninitializedField("update_gradients")); 5136 AssertIndexRange(component, fe->n_components()); 5137 Assert(present_cell.get() != nullptr, 5138 ExcMessage("FEValues object is not reinit'ed to any cell")); 5139 // check whether the shape function 5140 // is non-zero at all within 5141 // this component: 5142 if (fe->get_nonzero_components(i)[component] == false) 5143 return Tensor<1, spacedim>(); 5144 5145 // look up the right row in the 5146 // table and take the data from 5147 // there 5148 const unsigned int row = 5149 this->finite_element_output 5150 .shape_function_to_row_table[i * fe->n_components() + component]; 5151 return this->finite_element_output.shape_gradients[row][j]; 5152 } 5153 5154 5155 5156 template <int dim, int spacedim> 5157 inline const Tensor<2, spacedim> & 5158 FEValuesBase<dim, spacedim>::shape_hessian(const unsigned int i, 5159 const unsigned int j) const 5160 { 5161 AssertIndexRange(i, fe->n_dofs_per_cell()); 5162 Assert(this->update_flags & update_hessians, 5163 ExcAccessToUninitializedField("update_hessians")); 5164 Assert(fe->is_primitive(i), ExcShapeFunctionNotPrimitive(i)); 5165 Assert(present_cell.get() != nullptr, 5166 ExcMessage("FEValues object is not reinit'ed to any cell")); 5167 // if the entire FE is primitive, 5168 // then we can take a short-cut: 5169 if (fe->is_primitive()) 5170 return this->finite_element_output.shape_hessians[i][j]; 5171 else 5172 { 5173 // otherwise, use the mapping 5174 // between shape function 5175 // numbers and rows. note that 5176 // by the assertions above, we 5177 // know that this particular 5178 // shape function is primitive, 5179 // so we can call 5180 // system_to_component_index 5181 const unsigned int row = 5182 this->finite_element_output 5183 .shape_function_to_row_table[i * fe->n_components() + 5184 fe->system_to_component_index(i).first]; 5185 return this->finite_element_output.shape_hessians[row][j]; 5186 } 5187 } 5188 5189 5190 5191 template <int dim, int spacedim> 5192 inline Tensor<2, spacedim> 5193 FEValuesBase<dim, spacedim>::shape_hessian_component( 5194 const unsigned int i, 5195 const unsigned int j, 5196 const unsigned int component) const 5197 { 5198 AssertIndexRange(i, fe->n_dofs_per_cell()); 5199 Assert(this->update_flags & update_hessians, 5200 ExcAccessToUninitializedField("update_hessians")); 5201 AssertIndexRange(component, fe->n_components()); 5202 Assert(present_cell.get() != nullptr, 5203 ExcMessage("FEValues object is not reinit'ed to any cell")); 5204 // check whether the shape function 5205 // is non-zero at all within 5206 // this component: 5207 if (fe->get_nonzero_components(i)[component] == false) 5208 return Tensor<2, spacedim>(); 5209 5210 // look up the right row in the 5211 // table and take the data from 5212 // there 5213 const unsigned int row = 5214 this->finite_element_output 5215 .shape_function_to_row_table[i * fe->n_components() + component]; 5216 return this->finite_element_output.shape_hessians[row][j]; 5217 } 5218 5219 5220 5221 template <int dim, int spacedim> 5222 inline const Tensor<3, spacedim> & 5223 FEValuesBase<dim, spacedim>::shape_3rd_derivative(const unsigned int i, 5224 const unsigned int j) const 5225 { 5226 AssertIndexRange(i, fe->n_dofs_per_cell()); 5227 Assert(this->update_flags & update_3rd_derivatives, 5228 ExcAccessToUninitializedField("update_3rd_derivatives")); 5229 Assert(fe->is_primitive(i), ExcShapeFunctionNotPrimitive(i)); 5230 Assert(present_cell.get() != nullptr, 5231 ExcMessage("FEValues object is not reinit'ed to any cell")); 5232 // if the entire FE is primitive, 5233 // then we can take a short-cut: 5234 if (fe->is_primitive()) 5235 return this->finite_element_output.shape_3rd_derivatives[i][j]; 5236 else 5237 { 5238 // otherwise, use the mapping 5239 // between shape function 5240 // numbers and rows. note that 5241 // by the assertions above, we 5242 // know that this particular 5243 // shape function is primitive, 5244 // so we can call 5245 // system_to_component_index 5246 const unsigned int row = 5247 this->finite_element_output 5248 .shape_function_to_row_table[i * fe->n_components() + 5249 fe->system_to_component_index(i).first]; 5250 return this->finite_element_output.shape_3rd_derivatives[row][j]; 5251 } 5252 } 5253 5254 5255 5256 template <int dim, int spacedim> 5257 inline Tensor<3, spacedim> 5258 FEValuesBase<dim, spacedim>::shape_3rd_derivative_component( 5259 const unsigned int i, 5260 const unsigned int j, 5261 const unsigned int component) const 5262 { 5263 AssertIndexRange(i, fe->n_dofs_per_cell()); 5264 Assert(this->update_flags & update_3rd_derivatives, 5265 ExcAccessToUninitializedField("update_3rd_derivatives")); 5266 AssertIndexRange(component, fe->n_components()); 5267 Assert(present_cell.get() != nullptr, 5268 ExcMessage("FEValues object is not reinit'ed to any cell")); 5269 // check whether the shape function 5270 // is non-zero at all within 5271 // this component: 5272 if (fe->get_nonzero_components(i)[component] == false) 5273 return Tensor<3, spacedim>(); 5274 5275 // look up the right row in the 5276 // table and take the data from 5277 // there 5278 const unsigned int row = 5279 this->finite_element_output 5280 .shape_function_to_row_table[i * fe->n_components() + component]; 5281 return this->finite_element_output.shape_3rd_derivatives[row][j]; 5282 } 5283 5284 5285 5286 template <int dim, int spacedim> 5287 inline const FiniteElement<dim, spacedim> & 5288 FEValuesBase<dim, spacedim>::get_fe() const 5289 { 5290 return *fe; 5291 } 5292 5293 5294 5295 template <int dim, int spacedim> 5296 inline const Mapping<dim, spacedim> & 5297 FEValuesBase<dim, spacedim>::get_mapping() const 5298 { 5299 return *mapping; 5300 } 5301 5302 5303 5304 template <int dim, int spacedim> 5305 inline UpdateFlags 5306 FEValuesBase<dim, spacedim>::get_update_flags() const 5307 { 5308 return this->update_flags; 5309 } 5310 5311 5312 5313 template <int dim, int spacedim> 5314 inline const std::vector<Point<spacedim>> & 5315 FEValuesBase<dim, spacedim>::get_quadrature_points() const 5316 { 5317 Assert(this->update_flags & update_quadrature_points, 5318 ExcAccessToUninitializedField("update_quadrature_points")); 5319 Assert(present_cell.get() != nullptr, 5320 ExcMessage("FEValues object is not reinit'ed to any cell")); 5321 return this->mapping_output.quadrature_points; 5322 } 5323 5324 5325 5326 template <int dim, int spacedim> 5327 inline const std::vector<double> & 5328 FEValuesBase<dim, spacedim>::get_JxW_values() const 5329 { 5330 Assert(this->update_flags & update_JxW_values, 5331 ExcAccessToUninitializedField("update_JxW_values")); 5332 Assert(present_cell.get() != nullptr, 5333 ExcMessage("FEValues object is not reinit'ed to any cell")); 5334 return this->mapping_output.JxW_values; 5335 } 5336 5337 5338 5339 template <int dim, int spacedim> 5340 inline const std::vector<DerivativeForm<1, dim, spacedim>> & 5341 FEValuesBase<dim, spacedim>::get_jacobians() const 5342 { 5343 Assert(this->update_flags & update_jacobians, 5344 ExcAccessToUninitializedField("update_jacobians")); 5345 Assert(present_cell.get() != nullptr, 5346 ExcMessage("FEValues object is not reinit'ed to any cell")); 5347 return this->mapping_output.jacobians; 5348 } 5349 5350 5351 5352 template <int dim, int spacedim> 5353 inline const std::vector<DerivativeForm<2, dim, spacedim>> & 5354 FEValuesBase<dim, spacedim>::get_jacobian_grads() const 5355 { 5356 Assert(this->update_flags & update_jacobian_grads, 5357 ExcAccessToUninitializedField("update_jacobians_grads")); 5358 Assert(present_cell.get() != nullptr, 5359 ExcMessage("FEValues object is not reinit'ed to any cell")); 5360 return this->mapping_output.jacobian_grads; 5361 } 5362 5363 5364 5365 template <int dim, int spacedim> 5366 inline const Tensor<3, spacedim> & 5367 FEValuesBase<dim, spacedim>::jacobian_pushed_forward_grad( 5368 const unsigned int i) const 5369 { 5370 Assert(this->update_flags & update_jacobian_pushed_forward_grads, 5371 ExcAccessToUninitializedField("update_jacobian_pushed_forward_grads")); 5372 Assert(present_cell.get() != nullptr, 5373 ExcMessage("FEValues object is not reinit'ed to any cell")); 5374 return this->mapping_output.jacobian_pushed_forward_grads[i]; 5375 } 5376 5377 5378 5379 template <int dim, int spacedim> 5380 inline const std::vector<Tensor<3, spacedim>> & 5381 FEValuesBase<dim, spacedim>::get_jacobian_pushed_forward_grads() const 5382 { 5383 Assert(this->update_flags & update_jacobian_pushed_forward_grads, 5384 ExcAccessToUninitializedField("update_jacobian_pushed_forward_grads")); 5385 Assert(present_cell.get() != nullptr, 5386 ExcMessage("FEValues object is not reinit'ed to any cell")); 5387 return this->mapping_output.jacobian_pushed_forward_grads; 5388 } 5389 5390 5391 5392 template <int dim, int spacedim> 5393 inline const DerivativeForm<3, dim, spacedim> & 5394 FEValuesBase<dim, spacedim>::jacobian_2nd_derivative(const unsigned int i) const 5395 { 5396 Assert(this->update_flags & update_jacobian_2nd_derivatives, 5397 ExcAccessToUninitializedField("update_jacobian_2nd_derivatives")); 5398 Assert(present_cell.get() != nullptr, 5399 ExcMessage("FEValues object is not reinit'ed to any cell")); 5400 return this->mapping_output.jacobian_2nd_derivatives[i]; 5401 } 5402 5403 5404 5405 template <int dim, int spacedim> 5406 inline const std::vector<DerivativeForm<3, dim, spacedim>> & 5407 FEValuesBase<dim, spacedim>::get_jacobian_2nd_derivatives() const 5408 { 5409 Assert(this->update_flags & update_jacobian_2nd_derivatives, 5410 ExcAccessToUninitializedField("update_jacobian_2nd_derivatives")); 5411 Assert(present_cell.get() != nullptr, 5412 ExcMessage("FEValues object is not reinit'ed to any cell")); 5413 return this->mapping_output.jacobian_2nd_derivatives; 5414 } 5415 5416 5417 5418 template <int dim, int spacedim> 5419 inline const Tensor<4, spacedim> & 5420 FEValuesBase<dim, spacedim>::jacobian_pushed_forward_2nd_derivative( 5421 const unsigned int i) const 5422 { 5423 Assert(this->update_flags & update_jacobian_pushed_forward_2nd_derivatives, 5424 ExcAccessToUninitializedField( 5425 "update_jacobian_pushed_forward_2nd_derivatives")); 5426 Assert(present_cell.get() != nullptr, 5427 ExcMessage("FEValues object is not reinit'ed to any cell")); 5428 return this->mapping_output.jacobian_pushed_forward_2nd_derivatives[i]; 5429 } 5430 5431 5432 5433 template <int dim, int spacedim> 5434 inline const std::vector<Tensor<4, spacedim>> & 5435 FEValuesBase<dim, spacedim>::get_jacobian_pushed_forward_2nd_derivatives() const 5436 { 5437 Assert(this->update_flags & update_jacobian_pushed_forward_2nd_derivatives, 5438 ExcAccessToUninitializedField( 5439 "update_jacobian_pushed_forward_2nd_derivatives")); 5440 Assert(present_cell.get() != nullptr, 5441 ExcMessage("FEValues object is not reinit'ed to any cell")); 5442 return this->mapping_output.jacobian_pushed_forward_2nd_derivatives; 5443 } 5444 5445 5446 5447 template <int dim, int spacedim> 5448 inline const DerivativeForm<4, dim, spacedim> & 5449 FEValuesBase<dim, spacedim>::jacobian_3rd_derivative(const unsigned int i) const 5450 { 5451 Assert(this->update_flags & update_jacobian_3rd_derivatives, 5452 ExcAccessToUninitializedField("update_jacobian_3rd_derivatives")); 5453 Assert(present_cell.get() != nullptr, 5454 ExcMessage("FEValues object is not reinit'ed to any cell")); 5455 return this->mapping_output.jacobian_3rd_derivatives[i]; 5456 } 5457 5458 5459 5460 template <int dim, int spacedim> 5461 inline const std::vector<DerivativeForm<4, dim, spacedim>> & 5462 FEValuesBase<dim, spacedim>::get_jacobian_3rd_derivatives() const 5463 { 5464 Assert(this->update_flags & update_jacobian_3rd_derivatives, 5465 ExcAccessToUninitializedField("update_jacobian_3rd_derivatives")); 5466 Assert(present_cell.get() != nullptr, 5467 ExcMessage("FEValues object is not reinit'ed to any cell")); 5468 return this->mapping_output.jacobian_3rd_derivatives; 5469 } 5470 5471 5472 5473 template <int dim, int spacedim> 5474 inline const Tensor<5, spacedim> & 5475 FEValuesBase<dim, spacedim>::jacobian_pushed_forward_3rd_derivative( 5476 const unsigned int i) const 5477 { 5478 Assert(this->update_flags & update_jacobian_pushed_forward_3rd_derivatives, 5479 ExcAccessToUninitializedField( 5480 "update_jacobian_pushed_forward_3rd_derivatives")); 5481 Assert(present_cell.get() != nullptr, 5482 ExcMessage("FEValues object is not reinit'ed to any cell")); 5483 return this->mapping_output.jacobian_pushed_forward_3rd_derivatives[i]; 5484 } 5485 5486 5487 5488 template <int dim, int spacedim> 5489 inline const std::vector<Tensor<5, spacedim>> & 5490 FEValuesBase<dim, spacedim>::get_jacobian_pushed_forward_3rd_derivatives() const 5491 { 5492 Assert(this->update_flags & update_jacobian_pushed_forward_3rd_derivatives, 5493 ExcAccessToUninitializedField( 5494 "update_jacobian_pushed_forward_3rd_derivatives")); 5495 Assert(present_cell.get() != nullptr, 5496 ExcMessage("FEValues object is not reinit'ed to any cell")); 5497 return this->mapping_output.jacobian_pushed_forward_3rd_derivatives; 5498 } 5499 5500 5501 5502 template <int dim, int spacedim> 5503 inline const std::vector<DerivativeForm<1, spacedim, dim>> & 5504 FEValuesBase<dim, spacedim>::get_inverse_jacobians() const 5505 { 5506 Assert(this->update_flags & update_inverse_jacobians, 5507 ExcAccessToUninitializedField("update_inverse_jacobians")); 5508 Assert(present_cell.get() != nullptr, 5509 ExcMessage("FEValues object is not reinit'ed to any cell")); 5510 return this->mapping_output.inverse_jacobians; 5511 } 5512 5513 5514 5515 template <int dim, int spacedim> 5516 inline std_cxx20::ranges::iota_view<unsigned int, unsigned int> 5517 FEValuesBase<dim, spacedim>::dof_indices() const 5518 { 5519 return {0U, dofs_per_cell}; 5520 } 5521 5522 5523 5524 template <int dim, int spacedim> 5525 inline std_cxx20::ranges::iota_view<unsigned int, unsigned int> 5526 FEValuesBase<dim, spacedim>::dof_indices_starting_at( 5527 const unsigned int start_dof_index) const 5528 { 5529 Assert(start_dof_index <= dofs_per_cell, 5530 ExcIndexRange(start_dof_index, 0, dofs_per_cell + 1)); 5531 return {start_dof_index, dofs_per_cell}; 5532 } 5533 5534 5535 5536 template <int dim, int spacedim> 5537 inline std_cxx20::ranges::iota_view<unsigned int, unsigned int> 5538 FEValuesBase<dim, spacedim>::dof_indices_ending_at( 5539 const unsigned int end_dof_index) const 5540 { 5541 Assert(end_dof_index < dofs_per_cell, 5542 ExcIndexRange(end_dof_index, 0, dofs_per_cell)); 5543 return {0U, end_dof_index + 1}; 5544 } 5545 5546 5547 5548 template <int dim, int spacedim> 5549 inline std_cxx20::ranges::iota_view<unsigned int, unsigned int> 5550 FEValuesBase<dim, spacedim>::quadrature_point_indices() const 5551 { 5552 return {0U, n_quadrature_points}; 5553 } 5554 5555 5556 5557 template <int dim, int spacedim> 5558 inline const Point<spacedim> & 5559 FEValuesBase<dim, spacedim>::quadrature_point(const unsigned int i) const 5560 { 5561 Assert(this->update_flags & update_quadrature_points, 5562 ExcAccessToUninitializedField("update_quadrature_points")); 5563 AssertIndexRange(i, this->mapping_output.quadrature_points.size()); 5564 Assert(present_cell.get() != nullptr, 5565 ExcMessage("FEValues object is not reinit'ed to any cell")); 5566 5567 return this->mapping_output.quadrature_points[i]; 5568 } 5569 5570 5571 5572 template <int dim, int spacedim> 5573 inline double 5574 FEValuesBase<dim, spacedim>::JxW(const unsigned int i) const 5575 { 5576 Assert(this->update_flags & update_JxW_values, 5577 ExcAccessToUninitializedField("update_JxW_values")); 5578 AssertIndexRange(i, this->mapping_output.JxW_values.size()); 5579 Assert(present_cell.get() != nullptr, 5580 ExcMessage("FEValues object is not reinit'ed to any cell")); 5581 5582 return this->mapping_output.JxW_values[i]; 5583 } 5584 5585 5586 5587 template <int dim, int spacedim> 5588 inline const DerivativeForm<1, dim, spacedim> & 5589 FEValuesBase<dim, spacedim>::jacobian(const unsigned int i) const 5590 { 5591 Assert(this->update_flags & update_jacobians, 5592 ExcAccessToUninitializedField("update_jacobians")); 5593 AssertIndexRange(i, this->mapping_output.jacobians.size()); 5594 Assert(present_cell.get() != nullptr, 5595 ExcMessage("FEValues object is not reinit'ed to any cell")); 5596 5597 return this->mapping_output.jacobians[i]; 5598 } 5599 5600 5601 5602 template <int dim, int spacedim> 5603 inline const DerivativeForm<2, dim, spacedim> & 5604 FEValuesBase<dim, spacedim>::jacobian_grad(const unsigned int i) const 5605 { 5606 Assert(this->update_flags & update_jacobian_grads, 5607 ExcAccessToUninitializedField("update_jacobians_grads")); 5608 AssertIndexRange(i, this->mapping_output.jacobian_grads.size()); 5609 Assert(present_cell.get() != nullptr, 5610 ExcMessage("FEValues object is not reinit'ed to any cell")); 5611 5612 return this->mapping_output.jacobian_grads[i]; 5613 } 5614 5615 5616 5617 template <int dim, int spacedim> 5618 inline const DerivativeForm<1, spacedim, dim> & 5619 FEValuesBase<dim, spacedim>::inverse_jacobian(const unsigned int i) const 5620 { 5621 Assert(this->update_flags & update_inverse_jacobians, 5622 ExcAccessToUninitializedField("update_inverse_jacobians")); 5623 AssertIndexRange(i, this->mapping_output.inverse_jacobians.size()); 5624 Assert(present_cell.get() != nullptr, 5625 ExcMessage("FEValues object is not reinit'ed to any cell")); 5626 5627 return this->mapping_output.inverse_jacobians[i]; 5628 } 5629 5630 5631 5632 template <int dim, int spacedim> 5633 inline const Tensor<1, spacedim> & 5634 FEValuesBase<dim, spacedim>::normal_vector(const unsigned int i) const 5635 { 5636 Assert(this->update_flags & update_normal_vectors, 5637 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 5638 "update_normal_vectors"))); 5639 AssertIndexRange(i, this->mapping_output.normal_vectors.size()); 5640 Assert(present_cell.get() != nullptr, 5641 ExcMessage("FEValues object is not reinit'ed to any cell")); 5642 5643 return this->mapping_output.normal_vectors[i]; 5644 } 5645 5646 5647 5648 /*--------------------- Inline functions: FEValues --------------------------*/ 5649 5650 5651 template <int dim, int spacedim> 5652 inline const Quadrature<dim> & 5653 FEValues<dim, spacedim>::get_quadrature() const 5654 { 5655 return quadrature; 5656 } 5657 5658 5659 5660 template <int dim, int spacedim> 5661 inline const FEValues<dim, spacedim> & 5662 FEValues<dim, spacedim>::get_present_fe_values() const 5663 { 5664 return *this; 5665 } 5666 5667 5668 /*---------------------- Inline functions: FEFaceValuesBase -----------------*/ 5669 5670 5671 template <int dim, int spacedim> 5672 inline unsigned int 5673 FEFaceValuesBase<dim, spacedim>::get_face_index() const 5674 { 5675 return present_face_index; 5676 } 5677 5678 5679 /*----------------------- Inline functions: FE*FaceValues -------------------*/ 5680 5681 template <int dim, int spacedim> 5682 inline const Quadrature<dim - 1> & 5683 FEFaceValuesBase<dim, spacedim>::get_quadrature() const 5684 { 5685 return quadrature; 5686 } 5687 5688 5689 5690 template <int dim, int spacedim> 5691 inline const FEFaceValues<dim, spacedim> & 5692 FEFaceValues<dim, spacedim>::get_present_fe_values() const 5693 { 5694 return *this; 5695 } 5696 5697 5698 5699 template <int dim, int spacedim> 5700 inline const FESubfaceValues<dim, spacedim> & 5701 FESubfaceValues<dim, spacedim>::get_present_fe_values() const 5702 { 5703 return *this; 5704 } 5705 5706 5707 5708 template <int dim, int spacedim> 5709 inline const Tensor<1, spacedim> & 5710 FEFaceValuesBase<dim, spacedim>::boundary_form(const unsigned int i) const 5711 { 5712 AssertIndexRange(i, this->mapping_output.boundary_forms.size()); 5713 Assert(this->update_flags & update_boundary_forms, 5714 (typename FEValuesBase<dim, spacedim>::ExcAccessToUninitializedField( 5715 "update_boundary_forms"))); 5716 5717 return this->mapping_output.boundary_forms[i]; 5718 } 5719 5720 #endif // DOXYGEN 5721 5722 DEAL_II_NAMESPACE_CLOSE 5723 5724 #endif 5725