1 // The libMesh Finite Element Library. 2 // Copyright (C) 2002-2020 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 3 4 // This library is free software; you can redistribute it and/or 5 // modify it under the terms of the GNU Lesser General Public 6 // License as published by the Free Software Foundation; either 7 // version 2.1 of the License, or (at your option) any later version. 8 9 // This library is distributed in the hope that it will be useful, 10 // but WITHOUT ANY WARRANTY; without even the implied warranty of 11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12 // Lesser General Public License for more details. 13 14 // You should have received a copy of the GNU Lesser General Public 15 // License along with this library; if not, write to the Free Software 16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 17 18 19 20 #ifndef LIBMESH_FE_ABSTRACT_H 21 #define LIBMESH_FE_ABSTRACT_H 22 23 // Local includes 24 #include "libmesh/reference_counted_object.h" 25 #include "libmesh/point.h" 26 #include "libmesh/vector_value.h" 27 #include "libmesh/fe_type.h" 28 #include "libmesh/fe_map.h" 29 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 30 #include "libmesh/tensor_value.h" 31 #endif 32 33 #ifdef LIBMESH_FORWARD_DECLARE_ENUMS 34 namespace libMesh 35 { 36 enum ElemType : int; 37 } 38 #else 39 #include "libmesh/enum_elem_type.h" 40 #endif 41 42 // C++ includes 43 #include <cstddef> 44 #include <vector> 45 #include <memory> 46 47 namespace libMesh 48 { 49 50 51 // forward declarations 52 template <typename T> class DenseMatrix; 53 template <typename T> class DenseVector; 54 class BoundaryInfo; 55 class DofConstraints; 56 class DofMap; 57 class Elem; 58 class MeshBase; 59 template <typename T> class NumericVector; 60 class QBase; 61 62 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 63 class NodeConstraints; 64 #endif 65 66 #ifdef LIBMESH_ENABLE_PERIODIC 67 class PeriodicBoundaries; 68 class PointLocatorBase; 69 #endif 70 71 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 72 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map> 73 class InfFE; 74 #endif 75 76 77 78 /** 79 * This class forms the foundation from which generic finite elements 80 * may be derived. In the current implementation, the templated 81 * derived class \p FE offers a wide variety of commonly used finite 82 * element concepts. Check there for details. Use the \p 83 * FEAbstract::build() method to create an object of any of the 84 * derived classes. 85 * 86 * \note In the present design, the number of virtual members is kept 87 * to a minimum for performance reasons, although this is not based on 88 * rigorous profiling. 89 * 90 * All calls to static members of the \p FE classes should be 91 * requested through the \p FEInterface. This interface class 92 * approximates runtime polymorphism for the templated finite element 93 * classes. Even internal library classes, like \p DofMap, request 94 * the number of DOFs through this interface class. This approach 95 * also enables the co-existence of various element-based schemes. 96 * 97 * \author Benjamin S. Kirk 98 * \date 2002 99 */ 100 class FEAbstract : public ReferenceCountedObject<FEAbstract> 101 { 102 protected: 103 104 /** 105 * Constructor. Optionally initializes required data 106 * structures. Protected so that this base class 107 * cannot be explicitly instantiated. 108 */ 109 FEAbstract (const unsigned int dim, 110 const FEType & fet); 111 112 public: 113 114 /** 115 * Destructor. 116 */ 117 virtual ~FEAbstract(); 118 119 /** 120 * Builds a specific finite element type. 121 * 122 * \returns A std::unique_ptr<FEAbstract> to the FE object to prevent 123 * memory leaks. 124 */ 125 static std::unique_ptr<FEAbstract> build (const unsigned int dim, 126 const FEType & type); 127 128 /** 129 * This is at the core of this class. Use this for each 130 * new element in the mesh. Reinitializes the requested physical 131 * element-dependent data based on the current element 132 * \p elem. By default the element data are computed at the quadrature 133 * points specified by the quadrature rule \p qrule, but any set 134 * of points on the reference element may be specified in the optional 135 * argument \p pts. 136 * 137 * \note The FE classes decide which data to initialize based on 138 * which accessor functions such as \p get_phi() or \p get_d2phi() 139 * have been called, so all such accessors should be called before 140 * the first \p reinit(). 141 */ 142 virtual void reinit (const Elem * elem, 143 const std::vector<Point> * const pts = nullptr, 144 const std::vector<Real> * const weights = nullptr) = 0; 145 146 /** 147 * Reinitializes all the physical element-dependent data based on 148 * the \p side of the element \p elem. The \p tolerance parameter 149 * is passed to the involved call to \p inverse_map(). By default the 150 * element data are computed at the quadrature points specified by the 151 * quadrature rule \p qrule, but any set of points on the reference 152 * \em side element may be specified in the optional argument \p pts. 153 */ 154 virtual void reinit (const Elem * elem, 155 const unsigned int side, 156 const Real tolerance = TOLERANCE, 157 const std::vector<Point> * const pts = nullptr, 158 const std::vector<Real> * const weights = nullptr) = 0; 159 160 /** 161 * Reinitializes all the physical element-dependent data based on 162 * the \p edge of the element \p elem. The \p tolerance parameter 163 * is passed to the involved call to \p inverse_map(). By default the 164 * element data are computed at the quadrature points specified by the 165 * quadrature rule \p qrule, but any set of points on the reference 166 * \em edge element may be specified in the optional argument \p pts. 167 */ 168 virtual void edge_reinit (const Elem * elem, 169 const unsigned int edge, 170 const Real tolerance = TOLERANCE, 171 const std::vector<Point> * pts = nullptr, 172 const std::vector<Real> * weights = nullptr) = 0; 173 174 /** 175 * Computes the reference space quadrature points on the side of 176 * an element based on the side quadrature points. 177 */ 178 virtual void side_map (const Elem * elem, 179 const Elem * side, 180 const unsigned int s, 181 const std::vector<Point> & reference_side_points, 182 std::vector<Point> & reference_points) = 0; 183 184 /** 185 * \returns \p true if the point p is located on the reference element 186 * for element type t, false otherwise. Since we are doing floating 187 * point comparisons here the parameter \p eps can be specified to 188 * indicate a tolerance. For example, \f$ x \le 1 \f$ becomes 189 * \f$ x \le 1 + \epsilon \f$. 190 */ 191 static bool on_reference_element(const Point & p, 192 const ElemType t, 193 const Real eps = TOLERANCE); 194 /** 195 * \returns The reference space coordinates of \p nodes based on the 196 * element type. 197 */ 198 static void get_refspace_nodes(const ElemType t, 199 std::vector<Point> & nodes); 200 201 202 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 203 /** 204 * Computes the nodal constraint contributions (for 205 * non-conforming adapted meshes), using Lagrange geometry 206 */ 207 static void compute_node_constraints (NodeConstraints & constraints, 208 const Elem * elem); 209 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS 210 211 #ifdef LIBMESH_ENABLE_PERIODIC 212 213 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 214 /** 215 * Computes the node position constraint equation contributions (for 216 * meshes with periodic boundary conditions) 217 */ 218 static void compute_periodic_node_constraints (NodeConstraints & constraints, 219 const PeriodicBoundaries & boundaries, 220 const MeshBase & mesh, 221 const PointLocatorBase * point_locator, 222 const Elem * elem); 223 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS 224 225 #endif // LIBMESH_ENABLE_PERIODIC 226 227 /** 228 * \returns the dimension of this FE 229 */ get_dim()230 unsigned int get_dim() const 231 { return dim; } 232 233 /** 234 * \returns nothing, but lets the FE know you're explicitly 235 * prerequesting calculations. This is useful when you only want 236 * the FE for n_quadrature_points, n_dofs_on_side, or other methods 237 * that don't require shape function calculations, but you don't 238 * want libMesh "backwards compatibility" mode to assume you've made 239 * no prerequests and need to calculate everything. 240 */ get_nothing()241 void get_nothing() const 242 { calculate_nothing = true; } 243 244 /** 245 * \returns The \p xyz spatial locations of the quadrature 246 * points on the element. 247 */ get_xyz()248 const std::vector<Point> & get_xyz() const 249 { calculate_map = true; return this->_fe_map->get_xyz(); } 250 251 /** 252 * \returns The element Jacobian times the quadrature weight for 253 * each quadrature point. 254 */ get_JxW()255 const std::vector<Real> & get_JxW() const 256 { calculate_map = true; return this->_fe_map->get_JxW(); } 257 258 /** 259 * \returns The element tangents in xi-direction at the quadrature 260 * points. 261 */ get_dxyzdxi()262 const std::vector<RealGradient> & get_dxyzdxi() const 263 { calculate_map = true; return this->_fe_map->get_dxyzdxi(); } 264 265 /** 266 * \returns The element tangents in eta-direction at the quadrature 267 * points. 268 */ get_dxyzdeta()269 const std::vector<RealGradient> & get_dxyzdeta() const 270 { calculate_map = true; return this->_fe_map->get_dxyzdeta(); } 271 272 /** 273 * \returns The element tangents in zeta-direction at the quadrature 274 * points. 275 */ get_dxyzdzeta()276 const std::vector<RealGradient> & get_dxyzdzeta() const 277 { return _fe_map->get_dxyzdzeta(); } 278 279 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 280 281 /** 282 * \returns The second partial derivatives in xi. 283 */ get_d2xyzdxi2()284 const std::vector<RealGradient> & get_d2xyzdxi2() const 285 { calculate_map = true; return this->_fe_map->get_d2xyzdxi2(); } 286 287 /** 288 * \returns The second partial derivatives in eta. 289 */ get_d2xyzdeta2()290 const std::vector<RealGradient> & get_d2xyzdeta2() const 291 { calculate_map = true; return this->_fe_map->get_d2xyzdeta2(); } 292 293 /** 294 * \returns The second partial derivatives in zeta. 295 */ get_d2xyzdzeta2()296 const std::vector<RealGradient> & get_d2xyzdzeta2() const 297 { calculate_map = true; return this->_fe_map->get_d2xyzdzeta2(); } 298 299 /** 300 * \returns The second partial derivatives in xi-eta. 301 */ get_d2xyzdxideta()302 const std::vector<RealGradient> & get_d2xyzdxideta() const 303 { calculate_map = true; return this->_fe_map->get_d2xyzdxideta(); } 304 305 /** 306 * \returns The second partial derivatives in xi-zeta. 307 */ get_d2xyzdxidzeta()308 const std::vector<RealGradient> & get_d2xyzdxidzeta() const 309 { calculate_map = true; return this->_fe_map->get_d2xyzdxidzeta(); } 310 311 /** 312 * \returns The second partial derivatives in eta-zeta. 313 */ get_d2xyzdetadzeta()314 const std::vector<RealGradient> & get_d2xyzdetadzeta() const 315 { calculate_map = true; return this->_fe_map->get_d2xyzdetadzeta(); } 316 317 #endif 318 319 /** 320 * \returns The dxi/dx entry in the transformation 321 * matrix from physical to local coordinates. 322 */ get_dxidx()323 const std::vector<Real> & get_dxidx() const 324 { calculate_map = true; return this->_fe_map->get_dxidx(); } 325 326 /** 327 * \returns The dxi/dy entry in the transformation 328 * matrix from physical to local coordinates. 329 */ get_dxidy()330 const std::vector<Real> & get_dxidy() const 331 { calculate_map = true; return this->_fe_map->get_dxidy(); } 332 333 /** 334 * \returns The dxi/dz entry in the transformation 335 * matrix from physical to local coordinates. 336 */ get_dxidz()337 const std::vector<Real> & get_dxidz() const 338 { calculate_map = true; return this->_fe_map->get_dxidz(); } 339 340 /** 341 * \returns The deta/dx entry in the transformation 342 * matrix from physical to local coordinates. 343 */ get_detadx()344 const std::vector<Real> & get_detadx() const 345 { calculate_map = true; return this->_fe_map->get_detadx(); } 346 347 /** 348 * \returns The deta/dy entry in the transformation 349 * matrix from physical to local coordinates. 350 */ get_detady()351 const std::vector<Real> & get_detady() const 352 { calculate_map = true; return this->_fe_map->get_detady(); } 353 354 /** 355 * \returns The deta/dz entry in the transformation 356 * matrix from physical to local coordinates. 357 */ get_detadz()358 const std::vector<Real> & get_detadz() const 359 { calculate_map = true; return this->_fe_map->get_detadz(); } 360 361 /** 362 * \returns The dzeta/dx entry in the transformation 363 * matrix from physical to local coordinates. 364 */ get_dzetadx()365 const std::vector<Real> & get_dzetadx() const 366 { calculate_map = true; return this->_fe_map->get_dzetadx(); } 367 368 /** 369 * \returns The dzeta/dy entry in the transformation 370 * matrix from physical to local coordinates. 371 */ get_dzetady()372 const std::vector<Real> & get_dzetady() const 373 { calculate_map = true; return this->_fe_map->get_dzetady(); } 374 375 /** 376 * \returns The dzeta/dz entry in the transformation 377 * matrix from physical to local coordinates. 378 */ get_dzetadz()379 const std::vector<Real> & get_dzetadz() const 380 { calculate_map = true; return this->_fe_map->get_dzetadz(); } 381 382 /** 383 * \returns The tangent vectors for face integration. 384 */ get_tangents()385 const std::vector<std::vector<Point>> & get_tangents() const 386 { calculate_map = true; return this->_fe_map->get_tangents(); } 387 388 /** 389 * \returns The outward pointing normal vectors for face integration. 390 */ get_normals()391 const std::vector<Point> & get_normals() const 392 { calculate_map = true; return this->_fe_map->get_normals(); } 393 394 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 395 /** 396 * \returns The curvatures for use in face integration. 397 */ get_curvatures()398 const std::vector<Real> & get_curvatures() const 399 { calculate_map = true; return this->_fe_map->get_curvatures();} 400 401 #endif 402 403 /** 404 * Provides the class with the quadrature rule. Implement 405 * this in derived classes. 406 */ 407 virtual void attach_quadrature_rule (QBase * q) = 0; 408 409 /** 410 * \returns The total number of approximation shape functions 411 * for the current element. Useful during matrix assembly. 412 * Implement this in derived classes. 413 */ 414 virtual unsigned int n_shape_functions () const = 0; 415 416 /** 417 * \returns The total number of quadrature points. Useful 418 * during matrix assembly. Implement this in derived classes. 419 */ 420 virtual unsigned int n_quadrature_points () const = 0; 421 422 /** 423 * \returns The element type that the current shape functions 424 * have been calculated for. Useful in determining when shape 425 * functions must be recomputed. 426 */ get_type()427 ElemType get_type() const { return elem_type; } 428 429 /** 430 * \returns The p refinement level that the current shape 431 * functions have been calculated for. 432 */ get_p_level()433 unsigned int get_p_level() const { return _p_level; } 434 435 /** 436 * \returns The FE Type (approximation order and family) of the finite element. 437 */ get_fe_type()438 FEType get_fe_type() const { return fe_type; } 439 440 /** 441 * \returns The approximation order of the finite element. 442 */ get_order()443 Order get_order() const { return static_cast<Order>(fe_type.order + _p_level); } 444 445 /** 446 * Sets the *base* FE order of the finite element. 447 */ set_fe_order(int new_order)448 void set_fe_order(int new_order) { fe_type.order = new_order; } 449 450 /** 451 * \returns The continuity level of the finite element. 452 */ 453 virtual FEContinuity get_continuity() const = 0; 454 455 /** 456 * \returns \p true if the finite element's higher order shape functions are 457 * hierarchic 458 */ 459 virtual bool is_hierarchic() const = 0; 460 461 /** 462 * \returns The finite element family of this element. 463 */ get_family()464 FEFamily get_family() const { return fe_type.family; } 465 466 /** 467 * \returns The mapping object 468 */ get_fe_map()469 const FEMap & get_fe_map() const { return *_fe_map.get(); } get_fe_map()470 FEMap & get_fe_map() { return *_fe_map.get(); } 471 472 /** 473 * Prints the Jacobian times the weight for each quadrature point. 474 */ 475 void print_JxW(std::ostream & os) const; 476 477 /** 478 * Prints the value of each shape function at each quadrature point. 479 * Implement in derived class since this depends on whether the element 480 * is vector-valued or not. 481 */ 482 virtual void print_phi(std::ostream & os) const =0; 483 virtual void print_dual_phi(std::ostream & os) const =0; 484 485 /** 486 * Prints the value of each shape function's derivative 487 * at each quadrature point. Implement in derived class since this 488 * depends on whether the element is vector-valued or not. 489 */ 490 virtual void print_dphi(std::ostream & os) const =0; 491 virtual void print_dual_dphi(std::ostream & os) const =0; 492 493 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 494 495 /** 496 * Prints the value of each shape function's second derivatives 497 * at each quadrature point. Implement in derived class since this 498 * depends on whether the element is vector-valued or not. 499 */ 500 virtual void print_d2phi(std::ostream & os) const =0; 501 virtual void print_dual_d2phi(std::ostream & os) const =0; 502 503 #endif 504 505 /** 506 * Prints the spatial location of each quadrature point 507 * (on the physical element). 508 */ 509 void print_xyz(std::ostream & os) const; 510 511 /** 512 * Prints all the relevant information about the current element. 513 */ 514 void print_info(std::ostream & os) const; 515 516 /** 517 * Same as above, but allows you to print to a stream. 518 */ 519 friend std::ostream & operator << (std::ostream & os, const FEAbstract & fe); 520 521 /** 522 * request phi calculations 523 */ 524 virtual void request_phi() const = 0; 525 virtual void request_dual_phi() const = 0; 526 527 /** 528 * request dphi calculations 529 */ 530 virtual void request_dphi() const = 0; 531 virtual void request_dual_dphi() const = 0; 532 533 /** 534 * set calculate_dual as needed 535 */ set_calculate_dual(const bool val)536 void set_calculate_dual(const bool val){calculate_dual = val; } 537 538 protected: 539 540 /** 541 * After having updated the jacobian and the transformation 542 * from local to global coordinates in \p FEMap::compute_map(), 543 * the first derivatives of the shape functions are 544 * transformed to global coordinates, giving \p dphi, 545 * \p dphidx, \p dphidy, and \p dphidz. This method 546 * should rarely be re-defined in derived classes, but 547 * still should be usable for children. Therefore, keep 548 * it protected. This needs to be implemented in the 549 * derived class since this function depends on whether 550 * the shape functions are vector-valued or not. 551 */ 552 virtual void compute_shape_functions(const Elem *, const std::vector<Point> & ) =0; 553 554 std::unique_ptr<FEMap> _fe_map; 555 556 557 /** 558 * The dimensionality of the object 559 */ 560 const unsigned int dim; 561 562 /** 563 * Have calculations with this object already been started? 564 * Then all get_* functions should already have been called. 565 */ 566 mutable bool calculations_started; 567 568 /** 569 * Are we calculating dual basis? 570 */ 571 mutable bool calculate_dual; 572 573 /** 574 * Are we potentially deliberately calculating nothing? 575 */ 576 mutable bool calculate_nothing; 577 578 /** 579 * Are we calculating mapping functions? 580 */ 581 mutable bool calculate_map; 582 583 /** 584 * Should we calculate shape functions? 585 */ 586 mutable bool calculate_phi; 587 588 /** 589 * Should we calculate shape function gradients? 590 */ 591 mutable bool calculate_dphi; 592 593 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 594 /** 595 * Should we calculate shape function hessians? 596 */ 597 mutable bool calculate_d2phi; 598 #else 599 // Otherwise several interfaces need to be redone. 600 const bool calculate_d2phi=false; 601 602 #endif 603 604 /** 605 * Should we calculate shape function curls? 606 */ 607 mutable bool calculate_curl_phi; 608 609 /** 610 * Should we calculate shape function divergences? 611 */ 612 mutable bool calculate_div_phi; 613 614 /** 615 * Should we calculate reference shape function gradients? 616 */ 617 mutable bool calculate_dphiref; 618 619 620 /** 621 * The finite element type for this object. 622 * 623 * \note This should be constant for the object. 624 */ 625 FEType fe_type; 626 627 /** 628 * The element type the current data structures are 629 * set up for. 630 */ 631 ElemType elem_type; 632 633 /** 634 * The p refinement level the current data structures are 635 * set up for. 636 */ 637 unsigned int _p_level; 638 639 /** 640 * A pointer to the quadrature rule employed 641 */ 642 QBase * qrule; 643 644 /** 645 * A flag indicating if current data structures 646 * correspond to quadrature rule points 647 */ 648 bool shapes_on_quadrature; 649 650 /** 651 * \returns \p true when the shape functions (for 652 * this \p FEFamily) depend on the particular 653 * element, and therefore needs to be re-initialized 654 * for each new element. \p false otherwise. 655 */ 656 virtual bool shapes_need_reinit() const = 0; 657 658 }; 659 660 } // namespace libMesh 661 662 #endif // LIBMESH_FE_ABSTRACT_H 663