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_MAP_H 21 #define LIBMESH_FE_MAP_H 22 23 // libMesh 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 29 // C++ includes 30 #include <memory> 31 32 namespace libMesh 33 { 34 35 // forward declarations 36 class Elem; 37 class Node; 38 39 /** 40 * Class contained in FE that encapsulates mapping (i.e. from physical 41 * space to reference space and vice-versa) quantities and operations. 42 * 43 * \author Paul Bauman 44 * \date 2012 45 * \brief Computes finite element mapping function values, gradients, etc. 46 */ 47 class FEMap 48 { 49 public: 50 51 FEMap(Real jtol = 0); ~FEMap()52 virtual ~FEMap(){} 53 54 static std::unique_ptr<FEMap> build(FEType fe_type); 55 56 static FEFamily map_fe_type(const Elem & elem); 57 58 template<unsigned int Dim> 59 void init_reference_to_physical_map(const std::vector<Point> & qp, 60 const Elem * elem); 61 62 /** 63 * Compute the jacobian and some other additional data fields at the 64 * single point with index p. Takes the integration weights as 65 * input, along with a pointer to the element and a list of points 66 * that contribute to the element. Also takes a boolean flag 67 * telling whether second derivatives should actually be computed. 68 */ 69 void compute_single_point_map(const unsigned int dim, 70 const std::vector<Real> & qw, 71 const Elem * elem, 72 unsigned int p, 73 const std::vector<const Node *> & elem_nodes, 74 bool compute_second_derivatives); 75 76 /** 77 * Compute the jacobian and some other additional 78 * data fields. Takes the integration weights 79 * as input, along with a pointer to the element. 80 * The element is assumed to have a constant Jacobian 81 */ 82 virtual void compute_affine_map(const unsigned int dim, 83 const std::vector<Real> & qw, 84 const Elem * elem); 85 86 /** 87 * Assign a fake jacobian and some other additional data fields. 88 * Takes the integration weights as input. For use on non-element 89 * evaluations. 90 */ 91 virtual void compute_null_map(const unsigned int dim, 92 const std::vector<Real> & qw); 93 94 95 /** 96 * Compute the jacobian and some other additional 97 * data fields. Takes the integration weights 98 * as input, along with a pointer to the element. 99 * Also takes a boolean parameter indicating whether second 100 * derivatives need to be calculated, allowing us to potentially 101 * skip unnecessary, expensive computations. 102 */ 103 virtual void compute_map(const unsigned int dim, 104 const std::vector<Real> & qw, 105 const Elem * elem, 106 bool calculate_d2phi); 107 108 /** 109 * Same as compute_map, but for a side. Useful for boundary integration. 110 */ 111 virtual void compute_face_map(int dim, 112 const std::vector<Real> & qw, 113 const Elem * side); 114 115 /** 116 * Same as before, but for an edge. Useful for some projections. 117 */ 118 void compute_edge_map(int dim, 119 const std::vector<Real> & qw, 120 const Elem * side); 121 122 /** 123 * Initializes the reference to physical element map for a side. 124 * This is used for boundary integration. 125 */ 126 template<unsigned int Dim> 127 void init_face_shape_functions(const std::vector<Point> & qp, 128 const Elem * side); 129 130 /** 131 * Same as before, but for an edge. This is used for some projection 132 * operators. 133 */ 134 template<unsigned int Dim> 135 void init_edge_shape_functions(const std::vector<Point> & qp, 136 const Elem * edge); 137 138 /** 139 * \returns The location (in physical space) of the point 140 * \p p located on the reference element. 141 */ 142 static Point map (const unsigned int dim, 143 const Elem * elem, 144 const Point & reference_point); 145 146 /** 147 * \returns component \p j of d(xyz)/d(xi eta zeta) (in physical 148 * space) of the point \p p located on the reference element. 149 */ 150 static Point map_deriv (const unsigned int dim, 151 const Elem * elem, 152 const unsigned int j, 153 const Point & reference_point); 154 155 /** 156 * \returns The location (on the reference element) of the 157 * point \p p located in physical space. This function requires 158 * inverting the (possibly nonlinear) transformation map, so 159 * it is not trivial. The optional parameter \p tolerance defines 160 * how close is "good enough." The map inversion iteration 161 * computes the sequence \f$ \{ p_n \} \f$, and the iteration is 162 * terminated when \f$ \|p - p_n\| < \mbox{\texttt{tolerance}} \f$ 163 * 164 * When secure == true, the following checks are enabled: 165 * 166 * In DEBUG mode only: 167 * .) dim==1,2: throw an error if det(J) <= 0 for any Newton iteration. 168 * .) Print warning for every iteration beyond max_cnt in which the Newton scheme has not converged. 169 * 170 * In !DEBUG mode only: 171 * .) Print a _single_ warning (1 warning for the entire simulation) 172 * if the Newton scheme ever requires more than max_cnt iterations. 173 * 174 * In both DEBUG and !DEBUG modes: 175 * .) dim==3: Throw an exception for singular Jacobian. 176 * .) Throw an error if the Newton iteration has not converged in 2*max_cnt iterations. 177 * 178 * In addition to the checks above, the "extra_checks" parameter can 179 * be used to turn on some additional tests. In particular, when 180 * extra_checks == true *and* compiled in DEBUG mode: 181 * .) Print a warning if p != map(inverse_map(p)) to within tolerance. 182 * .) Print a warning if the inverse-mapped point is not on the reference element to within tolerance. 183 */ 184 static Point inverse_map (const unsigned int dim, 185 const Elem * elem, 186 const Point & p, 187 const Real tolerance = TOLERANCE, 188 const bool secure = true, 189 const bool extra_checks = true); 190 191 /** 192 * Takes a number points in physical space (in the \p 193 * physical_points vector) and finds their location on the reference 194 * element for the input element \p elem. The values on the 195 * reference element are returned in the vector \p 196 * reference_points. The other parameters have the same meaning 197 * as the single Point version of inverse_map() above. 198 */ 199 static void inverse_map (unsigned int dim, 200 const Elem * elem, 201 const std::vector<Point> & physical_points, 202 std::vector<Point> & reference_points, 203 const Real tolerance = TOLERANCE, 204 const bool secure = true, 205 const bool extra_checks = true); 206 207 /** 208 * \returns The \p xyz spatial locations of the quadrature 209 * points on the element. 210 */ get_xyz()211 const std::vector<Point> & get_xyz() const 212 { libmesh_assert(!calculations_started || calculate_xyz); 213 calculate_xyz = true; return xyz; } 214 215 /** 216 * \returns The element Jacobian for each quadrature point. 217 */ get_jacobian()218 const std::vector<Real> & get_jacobian() const 219 { libmesh_assert(!calculations_started || calculate_dxyz); 220 calculate_dxyz = true; return jac; } 221 222 /** 223 * \returns The element Jacobian times the quadrature weight for 224 * each quadrature point. 225 */ get_JxW()226 const std::vector<Real> & get_JxW() const 227 { libmesh_assert(!calculations_started || calculate_dxyz); 228 calculate_dxyz = true; return JxW; } 229 230 /** 231 * \returns The element tangents in xi-direction at the quadrature 232 * points. 233 */ get_dxyzdxi()234 const std::vector<RealGradient> & get_dxyzdxi() const 235 { libmesh_assert(!calculations_started || calculate_dxyz); 236 calculate_dxyz = true; return dxyzdxi_map; } 237 238 /** 239 * \returns The element tangents in eta-direction at the quadrature 240 * points. 241 */ get_dxyzdeta()242 const std::vector<RealGradient> & get_dxyzdeta() const 243 { libmesh_assert(!calculations_started || calculate_dxyz); 244 calculate_dxyz = true; return dxyzdeta_map; } 245 246 /** 247 * \returns The element tangents in zeta-direction at the quadrature 248 * points. 249 */ get_dxyzdzeta()250 const std::vector<RealGradient> & get_dxyzdzeta() const 251 { libmesh_assert(!calculations_started || calculate_dxyz); 252 calculate_dxyz = true; return dxyzdzeta_map; } 253 254 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 255 256 /** 257 * \returns The second partial derivatives in xi. 258 */ get_d2xyzdxi2()259 const std::vector<RealGradient> & get_d2xyzdxi2() const 260 { libmesh_assert(!calculations_started || calculate_d2xyz); 261 calculate_d2xyz = true; return d2xyzdxi2_map; } 262 263 /** 264 * \returns The second partial derivatives in eta. 265 */ get_d2xyzdeta2()266 const std::vector<RealGradient> & get_d2xyzdeta2() const 267 { libmesh_assert(!calculations_started || calculate_d2xyz); 268 calculate_d2xyz = true; return d2xyzdeta2_map; } 269 270 /** 271 * \returns The second partial derivatives in zeta. 272 */ get_d2xyzdzeta2()273 const std::vector<RealGradient> & get_d2xyzdzeta2() const 274 { libmesh_assert(!calculations_started || calculate_d2xyz); 275 calculate_d2xyz = true; return d2xyzdzeta2_map; } 276 277 /** 278 * \returns The second partial derivatives in xi-eta. 279 */ get_d2xyzdxideta()280 const std::vector<RealGradient> & get_d2xyzdxideta() const 281 { libmesh_assert(!calculations_started || calculate_d2xyz); 282 calculate_d2xyz = true; return d2xyzdxideta_map; } 283 284 /** 285 * \returns The second partial derivatives in xi-zeta. 286 */ get_d2xyzdxidzeta()287 const std::vector<RealGradient> & get_d2xyzdxidzeta() const 288 { libmesh_assert(!calculations_started || calculate_d2xyz); 289 calculate_d2xyz = true; return d2xyzdxidzeta_map; } 290 291 /** 292 * \returns The second partial derivatives in eta-zeta. 293 */ get_d2xyzdetadzeta()294 const std::vector<RealGradient> & get_d2xyzdetadzeta() const 295 { libmesh_assert(!calculations_started || calculate_d2xyz); 296 calculate_d2xyz = true; return d2xyzdetadzeta_map; } 297 298 #endif 299 300 /** 301 * \returns The dxi/dx entry in the transformation 302 * matrix from physical to local coordinates. 303 */ get_dxidx()304 const std::vector<Real> & get_dxidx() const 305 { libmesh_assert(!calculations_started || calculate_dxyz); 306 calculate_dxyz = true; return dxidx_map; } 307 308 /** 309 * \returns The dxi/dy entry in the transformation 310 * matrix from physical to local coordinates. 311 */ get_dxidy()312 const std::vector<Real> & get_dxidy() const 313 { libmesh_assert(!calculations_started || calculate_dxyz); 314 calculate_dxyz = true; return dxidy_map; } 315 316 /** 317 * \returns The dxi/dz entry in the transformation 318 * matrix from physical to local coordinates. 319 */ get_dxidz()320 const std::vector<Real> & get_dxidz() const 321 { libmesh_assert(!calculations_started || calculate_dxyz); 322 calculate_dxyz = true; return dxidz_map; } 323 324 /** 325 * \returns The deta/dx entry in the transformation 326 * matrix from physical to local coordinates. 327 */ get_detadx()328 const std::vector<Real> & get_detadx() const 329 { libmesh_assert(!calculations_started || calculate_dxyz); 330 calculate_dxyz = true; return detadx_map; } 331 332 /** 333 * \returns The deta/dy entry in the transformation 334 * matrix from physical to local coordinates. 335 */ get_detady()336 const std::vector<Real> & get_detady() const 337 { libmesh_assert(!calculations_started || calculate_dxyz); 338 calculate_dxyz = true; return detady_map; } 339 340 /** 341 * \returns The deta/dz entry in the transformation 342 * matrix from physical to local coordinates. 343 */ get_detadz()344 const std::vector<Real> & get_detadz() const 345 { libmesh_assert(!calculations_started || calculate_dxyz); 346 calculate_dxyz = true; return detadz_map; } 347 348 /** 349 * \returns The dzeta/dx entry in the transformation 350 * matrix from physical to local coordinates. 351 */ get_dzetadx()352 const std::vector<Real> & get_dzetadx() const 353 { libmesh_assert(!calculations_started || calculate_dxyz); 354 calculate_dxyz = true; return dzetadx_map; } 355 356 /** 357 * \returns The dzeta/dy entry in the transformation 358 * matrix from physical to local coordinates. 359 */ get_dzetady()360 const std::vector<Real> & get_dzetady() const 361 { libmesh_assert(!calculations_started || calculate_dxyz); 362 calculate_dxyz = true; return dzetady_map; } 363 364 /** 365 * \returns The dzeta/dz entry in the transformation 366 * matrix from physical to local coordinates. 367 */ get_dzetadz()368 const std::vector<Real> & get_dzetadz() const 369 { libmesh_assert(!calculations_started || calculate_dxyz); 370 calculate_dxyz = true; return dzetadz_map; } 371 372 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 373 /** 374 * Second derivatives of "xi" reference coordinate wrt physical coordinates. 375 */ get_d2xidxyz2()376 const std::vector<std::vector<Real>> & get_d2xidxyz2() const 377 { libmesh_assert(!calculations_started || calculate_d2xyz); 378 calculate_d2xyz = true; return d2xidxyz2_map; } 379 380 /** 381 * Second derivatives of "eta" reference coordinate wrt physical coordinates. 382 */ get_d2etadxyz2()383 const std::vector<std::vector<Real>> & get_d2etadxyz2() const 384 { libmesh_assert(!calculations_started || calculate_d2xyz); 385 calculate_d2xyz = true; return d2etadxyz2_map; } 386 387 /** 388 * Second derivatives of "zeta" reference coordinate wrt physical coordinates. 389 */ get_d2zetadxyz2()390 const std::vector<std::vector<Real>> & get_d2zetadxyz2() const 391 { libmesh_assert(!calculations_started || calculate_d2xyz); 392 calculate_d2xyz = true; return d2zetadxyz2_map; } 393 #endif 394 395 /** 396 * \returns The reference to physical map for the side/edge 397 */ get_psi()398 const std::vector<std::vector<Real>> & get_psi() const 399 { return psi_map; } 400 401 /** 402 * \returns The reference to physical map for the element 403 */ get_phi_map()404 const std::vector<std::vector<Real>> & get_phi_map() const 405 { libmesh_assert(!calculations_started || calculate_xyz); 406 calculate_xyz = true; return phi_map; } 407 408 /** 409 * \returns The reference to physical map derivative 410 */ get_dphidxi_map()411 const std::vector<std::vector<Real>> & get_dphidxi_map() const 412 { libmesh_assert(!calculations_started || calculate_dxyz); 413 calculate_dxyz = true; return dphidxi_map; } 414 415 /** 416 * \returns The reference to physical map derivative 417 */ get_dphideta_map()418 const std::vector<std::vector<Real>> & get_dphideta_map() const 419 { libmesh_assert(!calculations_started || calculate_dxyz); 420 calculate_dxyz = true; return dphideta_map; } 421 422 /** 423 * \returns The reference to physical map derivative 424 */ get_dphidzeta_map()425 const std::vector<std::vector<Real>> & get_dphidzeta_map() const 426 { libmesh_assert(!calculations_started || calculate_dxyz); 427 calculate_dxyz = true; return dphidzeta_map; } 428 429 /** 430 * \returns The tangent vectors for face integration. 431 */ get_tangents()432 const std::vector<std::vector<Point>> & get_tangents() const 433 { libmesh_assert(!calculations_started || calculate_dxyz); 434 calculate_dxyz = true; return tangents; } 435 436 /** 437 * \returns The outward pointing normal vectors for face integration. 438 */ get_normals()439 const std::vector<Point> & get_normals() const 440 { libmesh_assert(!calculations_started || calculate_dxyz); 441 calculate_dxyz = true; return normals; } 442 443 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 444 445 /** 446 * \returns The curvatures for use in face integration. 447 */ get_curvatures()448 const std::vector<Real> & get_curvatures() const 449 { libmesh_assert(!calculations_started || calculate_d2xyz); 450 calculate_d2xyz = true; return curvatures;} 451 452 #endif 453 454 /** 455 * Prints the Jacobian times the weight for each quadrature point. 456 */ 457 void print_JxW(std::ostream & os) const; 458 459 /** 460 * Prints the spatial location of each quadrature point 461 * (on the physical element). 462 */ 463 void print_xyz(std::ostream & os) const; 464 465 /* FIXME: PB: The public functions below break encapsulation! I'm not happy 466 about it, but I don't have time to redo the infinite element routines. 467 These are needed in inf_fe_boundary.C and inf_fe.C. A proper implementation 468 would subclass this class and hide all the radial function business. */ 469 /** 470 * \returns The reference to physical map for the side/edge 471 */ get_psi()472 std::vector<std::vector<Real>> & get_psi() 473 { return psi_map; } 474 475 /** 476 * \returns The reference to physical map derivative for the side/edge 477 */ get_dpsidxi()478 std::vector<std::vector<Real>> & get_dpsidxi() 479 { libmesh_assert(!calculations_started || calculate_dxyz); 480 calculate_dxyz = true; return dpsidxi_map; } 481 get_dpsidxi()482 const std::vector<std::vector<Real>> & get_dpsidxi() const 483 { libmesh_assert(!calculations_started || calculate_dxyz); 484 calculate_dxyz = true; return dpsidxi_map; } 485 486 /** 487 * \returns The reference to physical map derivative for the side/edge 488 */ get_dpsideta()489 std::vector<std::vector<Real>> & get_dpsideta() 490 { libmesh_assert(!calculations_started || calculate_dxyz); 491 calculate_dxyz = true; return dpsideta_map; } 492 get_dpsideta()493 const std::vector<std::vector<Real>> & get_dpsideta() const 494 { libmesh_assert(!calculations_started || calculate_dxyz); 495 calculate_dxyz = true; return dpsideta_map; } 496 497 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 498 499 /** 500 * \returns The reference to physical map 2nd derivative for the side/edge 501 */ get_d2psidxi2()502 std::vector<std::vector<Real>> & get_d2psidxi2() 503 { libmesh_assert(!calculations_started || calculate_d2xyz); 504 calculate_d2xyz = true; return d2psidxi2_map; } 505 506 /** 507 * \returns const reference to physical map 2nd derivative for the side/edge 508 */ get_d2psidxi2()509 const std::vector<std::vector<Real>> & get_d2psidxi2() const 510 { libmesh_assert(!calculations_started || calculate_d2xyz); 511 calculate_d2xyz = true; return d2psidxi2_map; } 512 513 /** 514 * \returns The reference to physical map 2nd derivative for the side/edge 515 */ get_d2psidxideta()516 std::vector<std::vector<Real>> & get_d2psidxideta() 517 { libmesh_assert(!calculations_started || calculate_d2xyz); 518 calculate_d2xyz = true; return d2psidxideta_map; } 519 520 /** 521 * \returns const reference to physical map 2nd derivative for the side/edge 522 */ get_d2psidxideta()523 const std::vector<std::vector<Real>> & get_d2psidxideta() const 524 { libmesh_assert(!calculations_started || calculate_d2xyz); 525 calculate_d2xyz = true; return d2psidxideta_map; } 526 527 /** 528 * \returns The reference to physical map 2nd derivative for the side/edge 529 */ get_d2psideta2()530 std::vector<std::vector<Real>> & get_d2psideta2() 531 { libmesh_assert(!calculations_started || calculate_d2xyz); 532 calculate_d2xyz = true; return d2psideta2_map; } 533 534 535 /** 536 * \returns const reference to physical map 2nd derivative for the side/edge 537 */ get_d2psideta2()538 const std::vector<std::vector<Real>> & get_d2psideta2() const 539 { libmesh_assert(!calculations_started || calculate_d2xyz); 540 calculate_d2xyz = true; return d2psideta2_map; } 541 542 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES 543 544 /** 545 * \returns The reference to physical map for the element 546 */ get_phi_map()547 std::vector<std::vector<Real>> & get_phi_map() 548 { libmesh_assert(!calculations_started || calculate_xyz); 549 calculate_xyz = true; return phi_map; } 550 551 /** 552 * \returns The reference to physical map derivative 553 */ get_dphidxi_map()554 std::vector<std::vector<Real>> & get_dphidxi_map() 555 { libmesh_assert(!calculations_started || calculate_dxyz); 556 calculate_dxyz = true; return dphidxi_map; } 557 558 /** 559 * \returns The reference to physical map derivative 560 */ get_dphideta_map()561 std::vector<std::vector<Real>> & get_dphideta_map() 562 { libmesh_assert(!calculations_started || calculate_dxyz); 563 calculate_dxyz = true; return dphideta_map; } 564 565 /** 566 * \returns The reference to physical map derivative 567 */ get_dphidzeta_map()568 std::vector<std::vector<Real>> & get_dphidzeta_map() 569 { libmesh_assert(!calculations_started || calculate_dxyz); 570 calculate_dxyz = true; return dphidzeta_map; } 571 572 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 573 /** 574 * \returns The reference to physical map 2nd derivative 575 */ get_d2phidxi2_map()576 std::vector<std::vector<Real>> & get_d2phidxi2_map() 577 { libmesh_assert(!calculations_started || calculate_d2xyz); 578 calculate_d2xyz = true; return d2phidxi2_map; } 579 580 /** 581 * \returns The reference to physical map 2nd derivative 582 */ get_d2phidxideta_map()583 std::vector<std::vector<Real>> & get_d2phidxideta_map() 584 { libmesh_assert(!calculations_started || calculate_d2xyz); 585 calculate_d2xyz = true; return d2phidxideta_map; } 586 587 /** 588 * \returns The reference to physical map 2nd derivative 589 */ get_d2phidxidzeta_map()590 std::vector<std::vector<Real>> & get_d2phidxidzeta_map() 591 { libmesh_assert(!calculations_started || calculate_d2xyz); 592 calculate_d2xyz = true; return d2phidxidzeta_map; } 593 594 /** 595 * \returns The reference to physical map 2nd derivative 596 */ get_d2phideta2_map()597 std::vector<std::vector<Real>> & get_d2phideta2_map() 598 { libmesh_assert(!calculations_started || calculate_d2xyz); 599 calculate_d2xyz = true; return d2phideta2_map; } 600 601 /** 602 * \returns The reference to physical map 2nd derivative 603 */ get_d2phidetadzeta_map()604 std::vector<std::vector<Real>> & get_d2phidetadzeta_map() 605 { libmesh_assert(!calculations_started || calculate_d2xyz); 606 calculate_d2xyz = true; return d2phidetadzeta_map; } 607 608 /** 609 * \returns The reference to physical map 2nd derivative 610 */ get_d2phidzeta2_map()611 std::vector<std::vector<Real>> & get_d2phidzeta2_map() 612 { libmesh_assert(!calculations_started || calculate_d2xyz); 613 calculate_d2xyz = true; return d2phidzeta2_map; } 614 #endif 615 616 /* FIXME: PB: This function breaks encapsulation! Needed in FE<>::reinit and 617 InfFE<>::reinit. Not sure yet if the algorithm can be redone to avoid this. */ 618 /** 619 * \returns Writable reference to the element Jacobian times 620 * the quadrature weight for each quadrature point. 621 */ get_JxW()622 std::vector<Real> & get_JxW() 623 { libmesh_assert(!calculations_started || calculate_dxyz); 624 calculate_dxyz = true; return JxW; } 625 626 /** 627 * Allows the user to prerequest additional calculations in between 628 * two calls to reinit(); 629 */ 630 void add_calculations(); 631 632 /** 633 * Set the Jacobian tolerance used for determining when the mapping fails. The mapping is 634 * determined to fail if jac <= jacobian_tolerance. 635 */ set_jacobian_tolerance(Real tol)636 void set_jacobian_tolerance(Real tol) { jacobian_tolerance = tol; } 637 638 protected: 639 640 /** 641 * Determine which values are to be calculated 642 */ determine_calculations()643 void determine_calculations() { 644 calculations_started = true; 645 646 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 647 // Second derivative calculations currently have first derivative 648 // calculations as a prerequisite 649 if (calculate_d2xyz) 650 calculate_dxyz = true; 651 #endif 652 } 653 654 /** 655 * A utility function for use by compute_*_map 656 */ 657 void resize_quadrature_map_vectors(const unsigned int dim, unsigned int n_qp); 658 659 /** 660 * Used in \p FEMap::compute_map(), which should be 661 * be usable in derived classes, and therefore protected. 662 * 663 * \returns The x value of the pth entry of the dxzydxi_map. 664 */ dxdxi_map(const unsigned int p)665 Real dxdxi_map(const unsigned int p) const { return dxyzdxi_map[p](0); } 666 667 /** 668 * Used in \p FEMap::compute_map(), which should be 669 * be usable in derived classes, and therefore protected. 670 * 671 * \returns The y value of the pth entry of the dxzydxi_map. 672 */ dydxi_map(const unsigned int p)673 Real dydxi_map(const unsigned int p) const { return dxyzdxi_map[p](1); } 674 675 /** 676 * Used in \p FEMap::compute_map(), which should be 677 * be usable in derived classes, and therefore protected. 678 * 679 * \returns The z value of the pth entry of the dxzydxi_map. 680 */ dzdxi_map(const unsigned int p)681 Real dzdxi_map(const unsigned int p) const { return dxyzdxi_map[p](2); } 682 683 /** 684 * Used in \p FEMap::compute_map(), which should be 685 * be usable in derived classes, and therefore protected. 686 * 687 * \returns The x value of the pth entry of the dxzydeta_map. 688 */ dxdeta_map(const unsigned int p)689 Real dxdeta_map(const unsigned int p) const { return dxyzdeta_map[p](0); } 690 691 /** 692 * Used in \p FEMap::compute_map(), which should be 693 * be usable in derived classes, and therefore protected. 694 * 695 * \returns The y value of the pth entry of the dxzydeta_map. 696 */ dydeta_map(const unsigned int p)697 Real dydeta_map(const unsigned int p) const { return dxyzdeta_map[p](1); } 698 699 /** 700 * Used in \p FEMap::compute_map(), which should be 701 * be usable in derived classes, and therefore protected. 702 * 703 * \returns The z value of the pth entry of the dxzydeta_map. 704 */ dzdeta_map(const unsigned int p)705 Real dzdeta_map(const unsigned int p) const { return dxyzdeta_map[p](2); } 706 707 /** 708 * Used in \p FEMap::compute_map(), which should be 709 * be usable in derived classes, and therefore protected. 710 * 711 * \returns The x value of the pth entry of the dxzydzeta_map. 712 */ dxdzeta_map(const unsigned int p)713 Real dxdzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](0); } 714 715 /** 716 * Used in \p FEMap::compute_map(), which should be 717 * be usable in derived classes, and therefore protected. 718 * 719 * \returns The y value of the pth entry of the dxzydzeta_map. 720 */ dydzeta_map(const unsigned int p)721 Real dydzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](1); } 722 723 /** 724 * Used in \p FEMap::compute_map(), which should be 725 * be usable in derived classes, and therefore protected. 726 * 727 * \returns The z value of the pth entry of the dxzydzeta_map. 728 */ dzdzeta_map(const unsigned int p)729 Real dzdzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](2); } 730 731 /** 732 * The spatial locations of the quadrature points 733 */ 734 std::vector<Point> xyz; 735 736 /** 737 * Vector of partial derivatives: 738 * d(x)/d(xi), d(y)/d(xi), d(z)/d(xi) 739 */ 740 std::vector<RealGradient> dxyzdxi_map; 741 742 /** 743 * Vector of partial derivatives: 744 * d(x)/d(eta), d(y)/d(eta), d(z)/d(eta) 745 */ 746 std::vector<RealGradient> dxyzdeta_map; 747 748 /** 749 * Vector of partial derivatives: 750 * d(x)/d(zeta), d(y)/d(zeta), d(z)/d(zeta) 751 */ 752 std::vector<RealGradient> dxyzdzeta_map; 753 754 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 755 756 /** 757 * Vector of second partial derivatives in xi: 758 * d^2(x)/d(xi)^2, d^2(y)/d(xi)^2, d^2(z)/d(xi)^2 759 */ 760 std::vector<RealGradient> d2xyzdxi2_map; 761 762 /** 763 * Vector of mixed second partial derivatives in xi-eta: 764 * d^2(x)/d(xi)d(eta) d^2(y)/d(xi)d(eta) d^2(z)/d(xi)d(eta) 765 */ 766 std::vector<RealGradient> d2xyzdxideta_map; 767 768 /** 769 * Vector of second partial derivatives in eta: 770 * d^2(x)/d(eta)^2 771 */ 772 std::vector<RealGradient> d2xyzdeta2_map; 773 774 /** 775 * Vector of second partial derivatives in xi-zeta: 776 * d^2(x)/d(xi)d(zeta), d^2(y)/d(xi)d(zeta), d^2(z)/d(xi)d(zeta) 777 */ 778 std::vector<RealGradient> d2xyzdxidzeta_map; 779 780 /** 781 * Vector of mixed second partial derivatives in eta-zeta: 782 * d^2(x)/d(eta)d(zeta) d^2(y)/d(eta)d(zeta) d^2(z)/d(eta)d(zeta) 783 */ 784 std::vector<RealGradient> d2xyzdetadzeta_map; 785 786 /** 787 * Vector of second partial derivatives in zeta: 788 * d^2(x)/d(zeta)^2 789 */ 790 std::vector<RealGradient> d2xyzdzeta2_map; 791 792 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES 793 794 /** 795 * Map for partial derivatives: 796 * d(xi)/d(x). Needed for the Jacobian. 797 */ 798 std::vector<Real> dxidx_map; 799 800 /** 801 * Map for partial derivatives: 802 * d(xi)/d(y). Needed for the Jacobian. 803 */ 804 std::vector<Real> dxidy_map; 805 806 /** 807 * Map for partial derivatives: 808 * d(xi)/d(z). Needed for the Jacobian. 809 */ 810 std::vector<Real> dxidz_map; 811 812 813 /** 814 * Map for partial derivatives: 815 * d(eta)/d(x). Needed for the Jacobian. 816 */ 817 std::vector<Real> detadx_map; 818 819 /** 820 * Map for partial derivatives: 821 * d(eta)/d(y). Needed for the Jacobian. 822 */ 823 std::vector<Real> detady_map; 824 825 /** 826 * Map for partial derivatives: 827 * d(eta)/d(z). Needed for the Jacobian. 828 */ 829 std::vector<Real> detadz_map; 830 831 832 /** 833 * Map for partial derivatives: 834 * d(zeta)/d(x). Needed for the Jacobian. 835 */ 836 std::vector<Real> dzetadx_map; 837 838 /** 839 * Map for partial derivatives: 840 * d(zeta)/d(y). Needed for the Jacobian. 841 */ 842 std::vector<Real> dzetady_map; 843 844 /** 845 * Map for partial derivatives: 846 * d(zeta)/d(z). Needed for the Jacobian. 847 */ 848 std::vector<Real> dzetadz_map; 849 850 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 851 /** 852 * Second derivatives of "xi" reference coordinate wrt physical coordinates. 853 * At each qp: (xi_{xx}, xi_{xy}, xi_{xz}, xi_{yy}, xi_{yz}, xi_{zz}) 854 */ 855 std::vector<std::vector<Real>> d2xidxyz2_map; 856 857 /** 858 * Second derivatives of "eta" reference coordinate wrt physical coordinates. 859 * At each qp: (eta_{xx}, eta_{xy}, eta_{xz}, eta_{yy}, eta_{yz}, eta_{zz}) 860 */ 861 std::vector<std::vector<Real>> d2etadxyz2_map; 862 863 /** 864 * Second derivatives of "zeta" reference coordinate wrt physical coordinates. 865 * At each qp: (zeta_{xx}, zeta_{xy}, zeta_{xz}, zeta_{yy}, zeta_{yz}, zeta_{zz}) 866 */ 867 std::vector<std::vector<Real>> d2zetadxyz2_map; 868 #endif 869 870 /** 871 * Map for the shape function phi. 872 */ 873 std::vector<std::vector<Real>> phi_map; 874 875 /** 876 * Map for the derivative, d(phi)/d(xi). 877 */ 878 std::vector<std::vector<Real>> dphidxi_map; 879 880 /** 881 * Map for the derivative, d(phi)/d(eta). 882 */ 883 std::vector<std::vector<Real>> dphideta_map; 884 885 /** 886 * Map for the derivative, d(phi)/d(zeta). 887 */ 888 std::vector<std::vector<Real>> dphidzeta_map; 889 890 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 891 892 /** 893 * Map for the second derivative, d^2(phi)/d(xi)^2. 894 */ 895 std::vector<std::vector<Real>> d2phidxi2_map; 896 897 /** 898 * Map for the second derivative, d^2(phi)/d(xi)d(eta). 899 */ 900 std::vector<std::vector<Real>> d2phidxideta_map; 901 902 /** 903 * Map for the second derivative, d^2(phi)/d(xi)d(zeta). 904 */ 905 std::vector<std::vector<Real>> d2phidxidzeta_map; 906 907 /** 908 * Map for the second derivative, d^2(phi)/d(eta)^2. 909 */ 910 std::vector<std::vector<Real>> d2phideta2_map; 911 912 /** 913 * Map for the second derivative, d^2(phi)/d(eta)d(zeta). 914 */ 915 std::vector<std::vector<Real>> d2phidetadzeta_map; 916 917 /** 918 * Map for the second derivative, d^2(phi)/d(zeta)^2. 919 */ 920 std::vector<std::vector<Real>> d2phidzeta2_map; 921 922 #endif 923 924 /** 925 * Map for the side shape functions, psi. 926 */ 927 std::vector<std::vector<Real>> psi_map; 928 929 /** 930 * Map for the derivative of the side functions, 931 * d(psi)/d(xi). 932 */ 933 std::vector<std::vector<Real>> dpsidxi_map; 934 935 /** 936 * Map for the derivative of the side function, 937 * d(psi)/d(eta). 938 */ 939 std::vector<std::vector<Real>> dpsideta_map; 940 941 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 942 943 /** 944 * Map for the second derivatives (in xi) of the 945 * side shape functions. Useful for computing 946 * the curvature at the quadrature points. 947 */ 948 std::vector<std::vector<Real>> d2psidxi2_map; 949 950 /** 951 * Map for the second (cross) derivatives in xi, eta 952 * of the side shape functions. Useful for 953 * computing the curvature at the quadrature points. 954 */ 955 std::vector<std::vector<Real>> d2psidxideta_map; 956 957 /** 958 * Map for the second derivatives (in eta) of the 959 * side shape functions. Useful for computing the 960 * curvature at the quadrature points. 961 */ 962 std::vector<std::vector<Real>> d2psideta2_map; 963 964 #endif 965 966 /** 967 * Tangent vectors on boundary at quadrature points. 968 */ 969 std::vector<std::vector<Point>> tangents; 970 971 /** 972 * Normal vectors on boundary at quadrature points 973 */ 974 std::vector<Point> normals; 975 976 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 977 /** 978 * The mean curvature (= one half the sum of the principal 979 * curvatures) on the boundary at the quadrature points. 980 * The mean curvature is a scalar value. 981 */ 982 std::vector<Real> curvatures; 983 984 #endif 985 986 /** 987 * Jacobian values at quadrature points 988 */ 989 std::vector<Real> jac; 990 991 /** 992 * Jacobian*Weight values at quadrature points 993 */ 994 std::vector<Real> JxW; 995 996 /** 997 * Have calculations with this object already been started? 998 * Then all get_* functions should already have been called. 999 */ 1000 mutable bool calculations_started; 1001 1002 /** 1003 * Should we calculate physical point locations? 1004 */ 1005 mutable bool calculate_xyz; 1006 1007 /** 1008 * Should we calculate mapping gradients? 1009 */ 1010 mutable bool calculate_dxyz; 1011 1012 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 1013 1014 /** 1015 * Should we calculate mapping hessians? 1016 */ 1017 mutable bool calculate_d2xyz; 1018 1019 #endif 1020 1021 /** 1022 * The Jacobian tolerance used for determining when the mapping fails. The mapping is 1023 * determined to fail if jac <= jacobian_tolerance. If not set by the user, this number 1024 * defaults to 0 1025 */ 1026 Real jacobian_tolerance; 1027 1028 private: 1029 /** 1030 * A helper function used by FEMap::compute_single_point_map() to 1031 * compute second derivatives of the inverse map. 1032 */ 1033 void compute_inverse_map_second_derivs(unsigned p); 1034 1035 /** 1036 * Work vector for compute_affine_map() 1037 */ 1038 std::vector<const Node *> _elem_nodes; 1039 }; 1040 1041 } 1042 1043 #endif // LIBMESH_FE_MAP_H 1044