fe_monomial_shape_3D.C
Go to the documentation of this file.
00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 00019 // C++ inlcludes 00020 00021 // Local includes 00022 #include "libmesh/fe.h" 00023 #include "libmesh/elem.h" 00024 00025 namespace libMesh 00026 { 00027 00028 00029 00030 00031 template <> 00032 Real FE<3,MONOMIAL>::shape(const ElemType, 00033 const Order libmesh_dbg_var(order), 00034 const unsigned int i, 00035 const Point& p) 00036 { 00037 #if LIBMESH_DIM == 3 00038 00039 const Real xi = p(0); 00040 const Real eta = p(1); 00041 const Real zeta = p(2); 00042 00043 libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)* 00044 (static_cast<unsigned int>(order)+2)* 00045 (static_cast<unsigned int>(order)+3)/6); 00046 00047 // monomials. since they are hierarchic we only need one case block. 00048 switch (i) 00049 { 00050 // constant 00051 case 0: 00052 return 1.; 00053 00054 // linears 00055 case 1: 00056 return xi; 00057 00058 case 2: 00059 return eta; 00060 00061 case 3: 00062 return zeta; 00063 00064 // quadratics 00065 case 4: 00066 return xi*xi; 00067 00068 case 5: 00069 return xi*eta; 00070 00071 case 6: 00072 return eta*eta; 00073 00074 case 7: 00075 return xi*zeta; 00076 00077 case 8: 00078 return zeta*eta; 00079 00080 case 9: 00081 return zeta*zeta; 00082 00083 // cubics 00084 case 10: 00085 return xi*xi*xi; 00086 00087 case 11: 00088 return xi*xi*eta; 00089 00090 case 12: 00091 return xi*eta*eta; 00092 00093 case 13: 00094 return eta*eta*eta; 00095 00096 case 14: 00097 return xi*xi*zeta; 00098 00099 case 15: 00100 return xi*eta*zeta; 00101 00102 case 16: 00103 return eta*eta*zeta; 00104 00105 case 17: 00106 return xi*zeta*zeta; 00107 00108 case 18: 00109 return eta*zeta*zeta; 00110 00111 case 19: 00112 return zeta*zeta*zeta; 00113 00114 // quartics 00115 case 20: 00116 return xi*xi*xi*xi; 00117 00118 case 21: 00119 return xi*xi*xi*eta; 00120 00121 case 22: 00122 return xi*xi*eta*eta; 00123 00124 case 23: 00125 return xi*eta*eta*eta; 00126 00127 case 24: 00128 return eta*eta*eta*eta; 00129 00130 case 25: 00131 return xi*xi*xi*zeta; 00132 00133 case 26: 00134 return xi*xi*eta*zeta; 00135 00136 case 27: 00137 return xi*eta*eta*zeta; 00138 00139 case 28: 00140 return eta*eta*eta*zeta; 00141 00142 case 29: 00143 return xi*xi*zeta*zeta; 00144 00145 case 30: 00146 return xi*eta*zeta*zeta; 00147 00148 case 31: 00149 return eta*eta*zeta*zeta; 00150 00151 case 32: 00152 return xi*zeta*zeta*zeta; 00153 00154 case 33: 00155 return eta*zeta*zeta*zeta; 00156 00157 case 34: 00158 return zeta*zeta*zeta*zeta; 00159 00160 default: 00161 unsigned int o = 0; 00162 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 00163 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 00164 unsigned int block=o, nz = 0; 00165 for (; block < i2; block += (o-nz+1)) { nz++; } 00166 const unsigned int nx = block - i2; 00167 const unsigned int ny = o - nx - nz; 00168 Real val = 1.; 00169 for (unsigned int index=0; index != nx; index++) 00170 val *= xi; 00171 for (unsigned int index=0; index != ny; index++) 00172 val *= eta; 00173 for (unsigned int index=0; index != nz; index++) 00174 val *= zeta; 00175 return val; 00176 } 00177 00178 #endif 00179 00180 libmesh_error_msg("We'll never get here!"); 00181 return 0.; 00182 } 00183 00184 00185 00186 template <> 00187 Real FE<3,MONOMIAL>::shape(const Elem* elem, 00188 const Order order, 00189 const unsigned int i, 00190 const Point& p) 00191 { 00192 libmesh_assert(elem); 00193 00194 // call the orientation-independent shape functions 00195 return FE<3,MONOMIAL>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p); 00196 } 00197 00198 00199 00200 template <> 00201 Real FE<3,MONOMIAL>::shape_deriv(const ElemType, 00202 const Order libmesh_dbg_var(order), 00203 const unsigned int i, 00204 const unsigned int j, 00205 const Point& p) 00206 { 00207 #if LIBMESH_DIM == 3 00208 00209 libmesh_assert_less (j, 3); 00210 00211 libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)* 00212 (static_cast<unsigned int>(order)+2)* 00213 (static_cast<unsigned int>(order)+3)/6); 00214 00215 00216 const Real xi = p(0); 00217 const Real eta = p(1); 00218 const Real zeta = p(2); 00219 00220 // monomials. since they are hierarchic we only need one case block. 00221 switch (j) 00222 { 00223 // d()/dxi 00224 case 0: 00225 { 00226 switch (i) 00227 { 00228 // constant 00229 case 0: 00230 return 0.; 00231 00232 // linear 00233 case 1: 00234 return 1.; 00235 00236 case 2: 00237 return 0.; 00238 00239 case 3: 00240 return 0.; 00241 00242 // quadratic 00243 case 4: 00244 return 2.*xi; 00245 00246 case 5: 00247 return eta; 00248 00249 case 6: 00250 return 0.; 00251 00252 case 7: 00253 return zeta; 00254 00255 case 8: 00256 return 0.; 00257 00258 case 9: 00259 return 0.; 00260 00261 // cubic 00262 case 10: 00263 return 3.*xi*xi; 00264 00265 case 11: 00266 return 2.*xi*eta; 00267 00268 case 12: 00269 return eta*eta; 00270 00271 case 13: 00272 return 0.; 00273 00274 case 14: 00275 return 2.*xi*zeta; 00276 00277 case 15: 00278 return eta*zeta; 00279 00280 case 16: 00281 return 0.; 00282 00283 case 17: 00284 return zeta*zeta; 00285 00286 case 18: 00287 return 0.; 00288 00289 case 19: 00290 return 0.; 00291 00292 // quartics 00293 case 20: 00294 return 4.*xi*xi*xi; 00295 00296 case 21: 00297 return 3.*xi*xi*eta; 00298 00299 case 22: 00300 return 2.*xi*eta*eta; 00301 00302 case 23: 00303 return eta*eta*eta; 00304 00305 case 24: 00306 return 0.; 00307 00308 case 25: 00309 return 3.*xi*xi*zeta; 00310 00311 case 26: 00312 return 2.*xi*eta*zeta; 00313 00314 case 27: 00315 return eta*eta*zeta; 00316 00317 case 28: 00318 return 0.; 00319 00320 case 29: 00321 return 2.*xi*zeta*zeta; 00322 00323 case 30: 00324 return eta*zeta*zeta; 00325 00326 case 31: 00327 return 0.; 00328 00329 case 32: 00330 return zeta*zeta*zeta; 00331 00332 case 33: 00333 return 0.; 00334 00335 case 34: 00336 return 0.; 00337 00338 default: 00339 unsigned int o = 0; 00340 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 00341 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 00342 unsigned int block=o, nz = 0; 00343 for (; block < i2; block += (o-nz+1)) { nz++; } 00344 const unsigned int nx = block - i2; 00345 const unsigned int ny = o - nx - nz; 00346 Real val = nx; 00347 for (unsigned int index=1; index < nx; index++) 00348 val *= xi; 00349 for (unsigned int index=0; index != ny; index++) 00350 val *= eta; 00351 for (unsigned int index=0; index != nz; index++) 00352 val *= zeta; 00353 return val; 00354 } 00355 } 00356 00357 00358 // d()/deta 00359 case 1: 00360 { 00361 switch (i) 00362 { 00363 // constant 00364 case 0: 00365 return 0.; 00366 00367 // linear 00368 case 1: 00369 return 0.; 00370 00371 case 2: 00372 return 1.; 00373 00374 case 3: 00375 return 0.; 00376 00377 // quadratic 00378 case 4: 00379 return 0.; 00380 00381 case 5: 00382 return xi; 00383 00384 case 6: 00385 return 2.*eta; 00386 00387 case 7: 00388 return 0.; 00389 00390 case 8: 00391 return zeta; 00392 00393 case 9: 00394 return 0.; 00395 00396 // cubic 00397 case 10: 00398 return 0.; 00399 00400 case 11: 00401 return xi*xi; 00402 00403 case 12: 00404 return 2.*xi*eta; 00405 00406 case 13: 00407 return 3.*eta*eta; 00408 00409 case 14: 00410 return 0.; 00411 00412 case 15: 00413 return xi*zeta; 00414 00415 case 16: 00416 return 2.*eta*zeta; 00417 00418 case 17: 00419 return 0.; 00420 00421 case 18: 00422 return zeta*zeta; 00423 00424 case 19: 00425 return 0.; 00426 00427 // quartics 00428 case 20: 00429 return 0.; 00430 00431 case 21: 00432 return xi*xi*xi; 00433 00434 case 22: 00435 return 2.*xi*xi*eta; 00436 00437 case 23: 00438 return 3.*xi*eta*eta; 00439 00440 case 24: 00441 return 4.*eta*eta*eta; 00442 00443 case 25: 00444 return 0.; 00445 00446 case 26: 00447 return xi*xi*zeta; 00448 00449 case 27: 00450 return 2.*xi*eta*zeta; 00451 00452 case 28: 00453 return 3.*eta*eta*zeta; 00454 00455 case 29: 00456 return 0.; 00457 00458 case 30: 00459 return xi*zeta*zeta; 00460 00461 case 31: 00462 return 2.*eta*zeta*zeta; 00463 00464 case 32: 00465 return 0.; 00466 00467 case 33: 00468 return zeta*zeta*zeta; 00469 00470 case 34: 00471 return 0.; 00472 00473 default: 00474 unsigned int o = 0; 00475 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 00476 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 00477 unsigned int block=o, nz = 0; 00478 for (; block < i2; block += (o-nz+1)) { nz++; } 00479 const unsigned int nx = block - i2; 00480 const unsigned int ny = o - nx - nz; 00481 Real val = ny; 00482 for (unsigned int index=0; index != nx; index++) 00483 val *= xi; 00484 for (unsigned int index=1; index < ny; index++) 00485 val *= eta; 00486 for (unsigned int index=0; index != nz; index++) 00487 val *= zeta; 00488 return val; 00489 } 00490 } 00491 00492 00493 // d()/dzeta 00494 case 2: 00495 { 00496 switch (i) 00497 { 00498 // constant 00499 case 0: 00500 return 0.; 00501 00502 // linear 00503 case 1: 00504 return 0.; 00505 00506 case 2: 00507 return 0.; 00508 00509 case 3: 00510 return 1.; 00511 00512 // quadratic 00513 case 4: 00514 return 0.; 00515 00516 case 5: 00517 return 0.; 00518 00519 case 6: 00520 return 0.; 00521 00522 case 7: 00523 return xi; 00524 00525 case 8: 00526 return eta; 00527 00528 case 9: 00529 return 2.*zeta; 00530 00531 // cubic 00532 case 10: 00533 return 0.; 00534 00535 case 11: 00536 return 0.; 00537 00538 case 12: 00539 return 0.; 00540 00541 case 13: 00542 return 0.; 00543 00544 case 14: 00545 return xi*xi; 00546 00547 case 15: 00548 return xi*eta; 00549 00550 case 16: 00551 return eta*eta; 00552 00553 case 17: 00554 return 2.*xi*zeta; 00555 00556 case 18: 00557 return 2.*eta*zeta; 00558 00559 case 19: 00560 return 3.*zeta*zeta; 00561 00562 // quartics 00563 case 20: 00564 return 0.; 00565 00566 case 21: 00567 return 0.; 00568 00569 case 22: 00570 return 0.; 00571 00572 case 23: 00573 return 0.; 00574 00575 case 24: 00576 return 0.; 00577 00578 case 25: 00579 return xi*xi*xi; 00580 00581 case 26: 00582 return xi*xi*eta; 00583 00584 case 27: 00585 return xi*eta*eta; 00586 00587 case 28: 00588 return eta*eta*eta; 00589 00590 case 29: 00591 return 2.*xi*xi*zeta; 00592 00593 case 30: 00594 return 2.*xi*eta*zeta; 00595 00596 case 31: 00597 return 2.*eta*eta*zeta; 00598 00599 case 32: 00600 return 3.*xi*zeta*zeta; 00601 00602 case 33: 00603 return 3.*eta*zeta*zeta; 00604 00605 case 34: 00606 return 4.*zeta*zeta*zeta; 00607 00608 default: 00609 unsigned int o = 0; 00610 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 00611 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 00612 unsigned int block=o, nz = 0; 00613 for (; block < i2; block += (o-nz+1)) { nz++; } 00614 const unsigned int nx = block - i2; 00615 const unsigned int ny = o - nx - nz; 00616 Real val = nz; 00617 for (unsigned int index=0; index != nx; index++) 00618 val *= xi; 00619 for (unsigned int index=0; index != ny; index++) 00620 val *= eta; 00621 for (unsigned int index=1; index < nz; index++) 00622 val *= zeta; 00623 return val; 00624 } 00625 } 00626 00627 default: 00628 libmesh_error_msg("Invalid shape function derivative j = " << j); 00629 } 00630 00631 #endif 00632 00633 libmesh_error_msg("We'll never get here!"); 00634 return 0.; 00635 } 00636 00637 00638 00639 template <> 00640 Real FE<3,MONOMIAL>::shape_deriv(const Elem* elem, 00641 const Order order, 00642 const unsigned int i, 00643 const unsigned int j, 00644 const Point& p) 00645 { 00646 libmesh_assert(elem); 00647 00648 // call the orientation-independent shape function derivatives 00649 return FE<3,MONOMIAL>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 00650 } 00651 00652 00653 00654 template <> 00655 Real FE<3,MONOMIAL>::shape_second_deriv(const ElemType, 00656 const Order libmesh_dbg_var(order), 00657 const unsigned int i, 00658 const unsigned int j, 00659 const Point& p) 00660 { 00661 #if LIBMESH_DIM == 3 00662 00663 libmesh_assert_less (j, 6); 00664 00665 libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)* 00666 (static_cast<unsigned int>(order)+2)* 00667 (static_cast<unsigned int>(order)+3)/6); 00668 00669 const Real xi = p(0); 00670 const Real eta = p(1); 00671 const Real zeta = p(2); 00672 00673 // monomials. since they are hierarchic we only need one case block. 00674 switch (j) 00675 { 00676 // d^2()/dxi^2 00677 case 0: 00678 { 00679 switch (i) 00680 { 00681 // constant 00682 case 0: 00683 00684 // linear 00685 case 1: 00686 case 2: 00687 case 3: 00688 return 0.; 00689 00690 // quadratic 00691 case 4: 00692 return 2.; 00693 00694 case 5: 00695 case 6: 00696 case 7: 00697 case 8: 00698 case 9: 00699 return 0.; 00700 00701 // cubic 00702 case 10: 00703 return 6.*xi; 00704 00705 case 11: 00706 return 2.*eta; 00707 00708 case 12: 00709 case 13: 00710 return 0.; 00711 00712 case 14: 00713 return 2.*zeta; 00714 00715 case 15: 00716 case 16: 00717 case 17: 00718 case 18: 00719 case 19: 00720 return 0.; 00721 00722 // quartics 00723 case 20: 00724 return 12.*xi*xi; 00725 00726 case 21: 00727 return 6.*xi*eta; 00728 00729 case 22: 00730 return 2.*eta*eta; 00731 00732 case 23: 00733 case 24: 00734 return 0.; 00735 00736 case 25: 00737 return 6.*xi*zeta; 00738 00739 case 26: 00740 return 2.*eta*zeta; 00741 00742 case 27: 00743 case 28: 00744 return 0.; 00745 00746 case 29: 00747 return 2.*zeta*zeta; 00748 00749 case 30: 00750 case 31: 00751 case 32: 00752 case 33: 00753 case 34: 00754 return 0.; 00755 00756 default: 00757 unsigned int o = 0; 00758 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 00759 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 00760 unsigned int block=o, nz = 0; 00761 for (; block < i2; block += (o-nz+1)) { nz++; } 00762 const unsigned int nx = block - i2; 00763 const unsigned int ny = o - nx - nz; 00764 Real val = nx * (nx - 1); 00765 for (unsigned int index=2; index < nx; index++) 00766 val *= xi; 00767 for (unsigned int index=0; index != ny; index++) 00768 val *= eta; 00769 for (unsigned int index=0; index != nz; index++) 00770 val *= zeta; 00771 return val; 00772 } 00773 } 00774 00775 00776 // d^2()/dxideta 00777 case 1: 00778 { 00779 switch (i) 00780 { 00781 // constant 00782 case 0: 00783 00784 // linear 00785 case 1: 00786 case 2: 00787 case 3: 00788 return 0.; 00789 00790 // quadratic 00791 case 4: 00792 return 0.; 00793 00794 case 5: 00795 return 1.; 00796 00797 case 6: 00798 case 7: 00799 case 8: 00800 case 9: 00801 return 0.; 00802 00803 // cubic 00804 case 10: 00805 return 0.; 00806 00807 case 11: 00808 return 2.*xi; 00809 00810 case 12: 00811 return 2.*eta; 00812 00813 case 13: 00814 case 14: 00815 return 0.; 00816 00817 case 15: 00818 return zeta; 00819 00820 case 16: 00821 case 17: 00822 case 18: 00823 case 19: 00824 return 0.; 00825 00826 // quartics 00827 case 20: 00828 return 0.; 00829 00830 case 21: 00831 return 3.*xi*xi; 00832 00833 case 22: 00834 return 4.*xi*eta; 00835 00836 case 23: 00837 return 3.*eta*eta; 00838 00839 case 24: 00840 case 25: 00841 return 0.; 00842 00843 case 26: 00844 return 2.*xi*zeta; 00845 00846 case 27: 00847 return 2.*eta*zeta; 00848 00849 case 28: 00850 case 29: 00851 return 0.; 00852 00853 case 30: 00854 return zeta*zeta; 00855 00856 case 31: 00857 case 32: 00858 case 33: 00859 case 34: 00860 return 0.; 00861 00862 default: 00863 unsigned int o = 0; 00864 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 00865 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 00866 unsigned int block=o, nz = 0; 00867 for (; block < i2; block += (o-nz+1)) { nz++; } 00868 const unsigned int nx = block - i2; 00869 const unsigned int ny = o - nx - nz; 00870 Real val = nx * ny; 00871 for (unsigned int index=1; index < nx; index++) 00872 val *= xi; 00873 for (unsigned int index=1; index < ny; index++) 00874 val *= eta; 00875 for (unsigned int index=0; index != nz; index++) 00876 val *= zeta; 00877 return val; 00878 } 00879 } 00880 00881 00882 // d^2()/deta^2 00883 case 2: 00884 { 00885 switch (i) 00886 { 00887 // constant 00888 case 0: 00889 00890 // linear 00891 case 1: 00892 case 2: 00893 case 3: 00894 return 0.; 00895 00896 // quadratic 00897 case 4: 00898 case 5: 00899 return 0.; 00900 00901 case 6: 00902 return 2.; 00903 00904 case 7: 00905 case 8: 00906 case 9: 00907 return 0.; 00908 00909 // cubic 00910 case 10: 00911 case 11: 00912 return 0.; 00913 00914 case 12: 00915 return 2.*xi; 00916 case 13: 00917 return 6.*eta; 00918 00919 case 14: 00920 case 15: 00921 return 0.; 00922 00923 case 16: 00924 return 2.*zeta; 00925 00926 case 17: 00927 case 18: 00928 case 19: 00929 return 0.; 00930 00931 // quartics 00932 case 20: 00933 case 21: 00934 return 0.; 00935 00936 case 22: 00937 return 2.*xi*xi; 00938 00939 case 23: 00940 return 6.*xi*eta; 00941 00942 case 24: 00943 return 12.*eta*eta; 00944 00945 case 25: 00946 case 26: 00947 return 0.; 00948 00949 case 27: 00950 return 2.*xi*zeta; 00951 00952 case 28: 00953 return 6.*eta*zeta; 00954 00955 case 29: 00956 case 30: 00957 return 0.; 00958 00959 case 31: 00960 return 2.*zeta*zeta; 00961 00962 case 32: 00963 case 33: 00964 case 34: 00965 return 0.; 00966 00967 default: 00968 unsigned int o = 0; 00969 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 00970 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 00971 unsigned int block=o, nz = 0; 00972 for (; block < i2; block += (o-nz+1)) { nz++; } 00973 const unsigned int nx = block - i2; 00974 const unsigned int ny = o - nx - nz; 00975 Real val = ny * (ny - 1); 00976 for (unsigned int index=0; index != nx; index++) 00977 val *= xi; 00978 for (unsigned int index=2; index < ny; index++) 00979 val *= eta; 00980 for (unsigned int index=0; index != nz; index++) 00981 val *= zeta; 00982 return val; 00983 } 00984 } 00985 00986 00987 // d^2()/dxidzeta 00988 case 3: 00989 { 00990 switch (i) 00991 { 00992 // constant 00993 case 0: 00994 00995 // linear 00996 case 1: 00997 case 2: 00998 case 3: 00999 return 0.; 01000 01001 // quadratic 01002 case 4: 01003 case 5: 01004 case 6: 01005 return 0.; 01006 01007 case 7: 01008 return 1.; 01009 01010 case 8: 01011 case 9: 01012 return 0.; 01013 01014 // cubic 01015 case 10: 01016 case 11: 01017 case 12: 01018 case 13: 01019 return 0.; 01020 01021 case 14: 01022 return 2.*xi; 01023 01024 case 15: 01025 return eta; 01026 01027 case 16: 01028 return 0.; 01029 01030 case 17: 01031 return 2.*zeta; 01032 01033 case 18: 01034 case 19: 01035 return 0.; 01036 01037 // quartics 01038 case 20: 01039 case 21: 01040 case 22: 01041 case 23: 01042 case 24: 01043 return 0.; 01044 01045 case 25: 01046 return 3.*xi*xi; 01047 01048 case 26: 01049 return 2.*xi*eta; 01050 01051 case 27: 01052 return eta*eta; 01053 01054 case 28: 01055 return 0.; 01056 01057 case 29: 01058 return 4.*xi*zeta; 01059 01060 case 30: 01061 return 2.*eta*zeta; 01062 01063 case 31: 01064 return 0.; 01065 01066 case 32: 01067 return 3.*zeta*zeta; 01068 01069 case 33: 01070 case 34: 01071 return 0.; 01072 01073 default: 01074 unsigned int o = 0; 01075 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 01076 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 01077 unsigned int block=o, nz = 0; 01078 for (; block < i2; block += (o-nz+1)) { nz++; } 01079 const unsigned int nx = block - i2; 01080 const unsigned int ny = o - nx - nz; 01081 Real val = nx * nz; 01082 for (unsigned int index=1; index < nx; index++) 01083 val *= xi; 01084 for (unsigned int index=0; index != ny; index++) 01085 val *= eta; 01086 for (unsigned int index=1; index < nz; index++) 01087 val *= zeta; 01088 return val; 01089 } 01090 } 01091 01092 // d^2()/detadzeta 01093 case 4: 01094 { 01095 switch (i) 01096 { 01097 // constant 01098 case 0: 01099 01100 // linear 01101 case 1: 01102 case 2: 01103 case 3: 01104 return 0.; 01105 01106 // quadratic 01107 case 4: 01108 case 5: 01109 case 6: 01110 case 7: 01111 return 0.; 01112 01113 case 8: 01114 return 1.; 01115 01116 case 9: 01117 return 0.; 01118 01119 // cubic 01120 case 10: 01121 case 11: 01122 case 12: 01123 case 13: 01124 case 14: 01125 return 0.; 01126 01127 case 15: 01128 return xi; 01129 01130 case 16: 01131 return 2.*eta; 01132 01133 case 17: 01134 return 0.; 01135 01136 case 18: 01137 return 2.*zeta; 01138 01139 case 19: 01140 return 0.; 01141 01142 // quartics 01143 case 20: 01144 case 21: 01145 case 22: 01146 case 23: 01147 case 24: 01148 case 25: 01149 return 0.; 01150 01151 case 26: 01152 return xi*xi; 01153 01154 case 27: 01155 return 2.*xi*eta; 01156 01157 case 28: 01158 return 3.*eta*eta; 01159 01160 case 29: 01161 return 0.; 01162 01163 case 30: 01164 return 2.*xi*zeta; 01165 01166 case 31: 01167 return 4.*eta*zeta; 01168 01169 case 32: 01170 return 0.; 01171 01172 case 33: 01173 return 3.*zeta*zeta; 01174 01175 case 34: 01176 return 0.; 01177 01178 default: 01179 unsigned int o = 0; 01180 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 01181 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 01182 unsigned int block=o, nz = 0; 01183 for (; block < i2; block += (o-nz+1)) { nz++; } 01184 const unsigned int nx = block - i2; 01185 const unsigned int ny = o - nx - nz; 01186 Real val = ny * nz; 01187 for (unsigned int index=0; index != nx; index++) 01188 val *= xi; 01189 for (unsigned int index=1; index < ny; index++) 01190 val *= eta; 01191 for (unsigned int index=1; index < nz; index++) 01192 val *= zeta; 01193 return val; 01194 } 01195 } 01196 01197 01198 // d^2()/dzeta^2 01199 case 5: 01200 { 01201 switch (i) 01202 { 01203 // constant 01204 case 0: 01205 01206 // linear 01207 case 1: 01208 case 2: 01209 case 3: 01210 return 0.; 01211 01212 // quadratic 01213 case 4: 01214 case 5: 01215 case 6: 01216 case 7: 01217 case 8: 01218 return 0.; 01219 01220 case 9: 01221 return 2.; 01222 01223 // cubic 01224 case 10: 01225 case 11: 01226 case 12: 01227 case 13: 01228 case 14: 01229 case 15: 01230 case 16: 01231 return 0.; 01232 01233 case 17: 01234 return 2.*xi; 01235 01236 case 18: 01237 return 2.*eta; 01238 01239 case 19: 01240 return 6.*zeta; 01241 01242 // quartics 01243 case 20: 01244 case 21: 01245 case 22: 01246 case 23: 01247 case 24: 01248 case 25: 01249 case 26: 01250 case 27: 01251 case 28: 01252 return 0.; 01253 01254 case 29: 01255 return 2.*xi*xi; 01256 01257 case 30: 01258 return 2.*xi*eta; 01259 01260 case 31: 01261 return 2.*eta*eta; 01262 01263 case 32: 01264 return 6.*xi*zeta; 01265 01266 case 33: 01267 return 6.*eta*zeta; 01268 01269 case 34: 01270 return 12.*zeta*zeta; 01271 01272 default: 01273 unsigned int o = 0; 01274 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 01275 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 01276 unsigned int block=o, nz = 0; 01277 for (; block < i2; block += (o-nz+1)) { nz++; } 01278 const unsigned int nx = block - i2; 01279 const unsigned int ny = o - nx - nz; 01280 Real val = nz * (nz - 1); 01281 for (unsigned int index=0; index != nx; index++) 01282 val *= xi; 01283 for (unsigned int index=0; index != ny; index++) 01284 val *= eta; 01285 for (unsigned int index=2; index < nz; index++) 01286 val *= zeta; 01287 return val; 01288 } 01289 } 01290 01291 default: 01292 libmesh_error_msg("Invalid j = " << j); 01293 } 01294 01295 #endif 01296 01297 libmesh_error_msg("We'll never get here!"); 01298 return 0.; 01299 } 01300 01301 01302 01303 template <> 01304 Real FE<3,MONOMIAL>::shape_second_deriv(const Elem* elem, 01305 const Order order, 01306 const unsigned int i, 01307 const unsigned int j, 01308 const Point& p) 01309 { 01310 libmesh_assert(elem); 01311 01312 // call the orientation-independent shape function derivatives 01313 return FE<3,MONOMIAL>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 01314 } 01315 01316 } // namespace libMesh
generated by
