fe_monomial_shape_3D.C
Go to the documentation of this file.00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2012 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(); 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 00628 #endif 00629 00630 libmesh_error(); 00631 return 0.; 00632 } 00633 00634 00635 00636 template <> 00637 Real FE<3,MONOMIAL>::shape_deriv(const Elem* elem, 00638 const Order order, 00639 const unsigned int i, 00640 const unsigned int j, 00641 const Point& p) 00642 { 00643 libmesh_assert(elem); 00644 00645 // call the orientation-independent shape function derivatives 00646 return FE<3,MONOMIAL>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 00647 } 00648 00649 00650 00651 template <> 00652 Real FE<3,MONOMIAL>::shape_second_deriv(const ElemType, 00653 const Order libmesh_dbg_var(order), 00654 const unsigned int i, 00655 const unsigned int j, 00656 const Point& p) 00657 { 00658 #if LIBMESH_DIM == 3 00659 00660 libmesh_assert_less (j, 6); 00661 00662 libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)* 00663 (static_cast<unsigned int>(order)+2)* 00664 (static_cast<unsigned int>(order)+3)/6); 00665 00666 const Real xi = p(0); 00667 const Real eta = p(1); 00668 const Real zeta = p(2); 00669 00670 // monomials. since they are hierarchic we only need one case block. 00671 switch (j) 00672 { 00673 // d^2()/dxi^2 00674 case 0: 00675 { 00676 switch (i) 00677 { 00678 // constant 00679 case 0: 00680 00681 // linear 00682 case 1: 00683 case 2: 00684 case 3: 00685 return 0.; 00686 00687 // quadratic 00688 case 4: 00689 return 2.; 00690 00691 case 5: 00692 case 6: 00693 case 7: 00694 case 8: 00695 case 9: 00696 return 0.; 00697 00698 // cubic 00699 case 10: 00700 return 6.*xi; 00701 00702 case 11: 00703 return 2.*eta; 00704 00705 case 12: 00706 case 13: 00707 return 0.; 00708 00709 case 14: 00710 return 2.*zeta; 00711 00712 case 15: 00713 case 16: 00714 case 17: 00715 case 18: 00716 case 19: 00717 return 0.; 00718 00719 // quartics 00720 case 20: 00721 return 12.*xi*xi; 00722 00723 case 21: 00724 return 6.*xi*eta; 00725 00726 case 22: 00727 return 2.*eta*eta; 00728 00729 case 23: 00730 case 24: 00731 return 0.; 00732 00733 case 25: 00734 return 6.*xi*zeta; 00735 00736 case 26: 00737 return 2.*eta*zeta; 00738 00739 case 27: 00740 case 28: 00741 return 0.; 00742 00743 case 29: 00744 return 2.*zeta*zeta; 00745 00746 case 30: 00747 case 31: 00748 case 32: 00749 case 33: 00750 case 34: 00751 return 0.; 00752 00753 default: 00754 unsigned int o = 0; 00755 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 00756 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 00757 unsigned int block=o, nz = 0; 00758 for (; block < i2; block += (o-nz+1)) { nz++; } 00759 const unsigned int nx = block - i2; 00760 const unsigned int ny = o - nx - nz; 00761 Real val = nx * (nx - 1); 00762 for (unsigned int index=2; index < nx; index++) 00763 val *= xi; 00764 for (unsigned int index=0; index != ny; index++) 00765 val *= eta; 00766 for (unsigned int index=0; index != nz; index++) 00767 val *= zeta; 00768 return val; 00769 } 00770 } 00771 00772 00773 // d^2()/dxideta 00774 case 1: 00775 { 00776 switch (i) 00777 { 00778 // constant 00779 case 0: 00780 00781 // linear 00782 case 1: 00783 case 2: 00784 case 3: 00785 return 0.; 00786 00787 // quadratic 00788 case 4: 00789 return 0.; 00790 00791 case 5: 00792 return 1.; 00793 00794 case 6: 00795 case 7: 00796 case 8: 00797 case 9: 00798 return 0.; 00799 00800 // cubic 00801 case 10: 00802 return 0.; 00803 00804 case 11: 00805 return 2.*xi; 00806 00807 case 12: 00808 return 2.*eta; 00809 00810 case 13: 00811 case 14: 00812 return 0.; 00813 00814 case 15: 00815 return zeta; 00816 00817 case 16: 00818 case 17: 00819 case 18: 00820 case 19: 00821 return 0.; 00822 00823 // quartics 00824 case 20: 00825 return 0.; 00826 00827 case 21: 00828 return 3.*xi*xi; 00829 00830 case 22: 00831 return 4.*xi*eta; 00832 00833 case 23: 00834 return 3.*eta*eta; 00835 00836 case 24: 00837 case 25: 00838 return 0.; 00839 00840 case 26: 00841 return 2.*xi*zeta; 00842 00843 case 27: 00844 return 2.*eta*zeta; 00845 00846 case 28: 00847 case 29: 00848 return 0.; 00849 00850 case 30: 00851 return zeta*zeta; 00852 00853 case 31: 00854 case 32: 00855 case 33: 00856 case 34: 00857 return 0.; 00858 00859 default: 00860 unsigned int o = 0; 00861 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 00862 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 00863 unsigned int block=o, nz = 0; 00864 for (; block < i2; block += (o-nz+1)) { nz++; } 00865 const unsigned int nx = block - i2; 00866 const unsigned int ny = o - nx - nz; 00867 Real val = nx * ny; 00868 for (unsigned int index=1; index < nx; index++) 00869 val *= xi; 00870 for (unsigned int index=1; index < ny; index++) 00871 val *= eta; 00872 for (unsigned int index=0; index != nz; index++) 00873 val *= zeta; 00874 return val; 00875 } 00876 } 00877 00878 00879 // d^2()/deta^2 00880 case 2: 00881 { 00882 switch (i) 00883 { 00884 // constant 00885 case 0: 00886 00887 // linear 00888 case 1: 00889 case 2: 00890 case 3: 00891 return 0.; 00892 00893 // quadratic 00894 case 4: 00895 case 5: 00896 return 0.; 00897 00898 case 6: 00899 return 2.; 00900 00901 case 7: 00902 case 8: 00903 case 9: 00904 return 0.; 00905 00906 // cubic 00907 case 10: 00908 case 11: 00909 return 0.; 00910 00911 case 12: 00912 return 2.*xi; 00913 case 13: 00914 return 6.*eta; 00915 00916 case 14: 00917 case 15: 00918 return 0.; 00919 00920 case 16: 00921 return 2.*zeta; 00922 00923 case 17: 00924 case 18: 00925 case 19: 00926 return 0.; 00927 00928 // quartics 00929 case 20: 00930 case 21: 00931 return 0.; 00932 00933 case 22: 00934 return 2.*xi*xi; 00935 00936 case 23: 00937 return 6.*xi*eta; 00938 00939 case 24: 00940 return 12.*eta*eta; 00941 00942 case 25: 00943 case 26: 00944 return 0.; 00945 00946 case 27: 00947 return 2.*xi*zeta; 00948 00949 case 28: 00950 return 6.*eta*zeta; 00951 00952 case 29: 00953 case 30: 00954 return 0.; 00955 00956 case 31: 00957 return 2.*zeta*zeta; 00958 00959 case 32: 00960 case 33: 00961 case 34: 00962 return 0.; 00963 00964 default: 00965 unsigned int o = 0; 00966 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 00967 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 00968 unsigned int block=o, nz = 0; 00969 for (; block < i2; block += (o-nz+1)) { nz++; } 00970 const unsigned int nx = block - i2; 00971 const unsigned int ny = o - nx - nz; 00972 Real val = ny * (ny - 1); 00973 for (unsigned int index=0; index != nx; index++) 00974 val *= xi; 00975 for (unsigned int index=2; index < ny; index++) 00976 val *= eta; 00977 for (unsigned int index=0; index != nz; index++) 00978 val *= zeta; 00979 return val; 00980 } 00981 } 00982 00983 00984 // d^2()/dxidzeta 00985 case 3: 00986 { 00987 switch (i) 00988 { 00989 // constant 00990 case 0: 00991 00992 // linear 00993 case 1: 00994 case 2: 00995 case 3: 00996 return 0.; 00997 00998 // quadratic 00999 case 4: 01000 case 5: 01001 case 6: 01002 return 0.; 01003 01004 case 7: 01005 return 1.; 01006 01007 case 8: 01008 case 9: 01009 return 0.; 01010 01011 // cubic 01012 case 10: 01013 case 11: 01014 case 12: 01015 case 13: 01016 return 0.; 01017 01018 case 14: 01019 return 2.*xi; 01020 01021 case 15: 01022 return eta; 01023 01024 case 16: 01025 return 0.; 01026 01027 case 17: 01028 return 2.*zeta; 01029 01030 case 18: 01031 case 19: 01032 return 0.; 01033 01034 // quartics 01035 case 20: 01036 case 21: 01037 case 22: 01038 case 23: 01039 case 24: 01040 return 0.; 01041 01042 case 25: 01043 return 3.*xi*xi; 01044 01045 case 26: 01046 return 2.*xi*eta; 01047 01048 case 27: 01049 return eta*eta; 01050 01051 case 28: 01052 return 0.; 01053 01054 case 29: 01055 return 4.*xi*zeta; 01056 01057 case 30: 01058 return 2.*eta*zeta; 01059 01060 case 31: 01061 return 0.; 01062 01063 case 32: 01064 return 3.*zeta*zeta; 01065 01066 case 33: 01067 case 34: 01068 return 0.; 01069 01070 default: 01071 unsigned int o = 0; 01072 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 01073 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 01074 unsigned int block=o, nz = 0; 01075 for (; block < i2; block += (o-nz+1)) { nz++; } 01076 const unsigned int nx = block - i2; 01077 const unsigned int ny = o - nx - nz; 01078 Real val = nx * nz; 01079 for (unsigned int index=1; index < nx; index++) 01080 val *= xi; 01081 for (unsigned int index=0; index != ny; index++) 01082 val *= eta; 01083 for (unsigned int index=1; index < nz; index++) 01084 val *= zeta; 01085 return val; 01086 } 01087 } 01088 01089 // d^2()/detadzeta 01090 case 4: 01091 { 01092 switch (i) 01093 { 01094 // constant 01095 case 0: 01096 01097 // linear 01098 case 1: 01099 case 2: 01100 case 3: 01101 return 0.; 01102 01103 // quadratic 01104 case 4: 01105 case 5: 01106 case 6: 01107 case 7: 01108 return 0.; 01109 01110 case 8: 01111 return 1.; 01112 01113 case 9: 01114 return 0.; 01115 01116 // cubic 01117 case 10: 01118 case 11: 01119 case 12: 01120 case 13: 01121 case 14: 01122 return 0.; 01123 01124 case 15: 01125 return xi; 01126 01127 case 16: 01128 return 2.*eta; 01129 01130 case 17: 01131 return 0.; 01132 01133 case 18: 01134 return 2.*zeta; 01135 01136 case 19: 01137 return 0.; 01138 01139 // quartics 01140 case 20: 01141 case 21: 01142 case 22: 01143 case 23: 01144 case 24: 01145 case 25: 01146 return 0.; 01147 01148 case 26: 01149 return xi*xi; 01150 01151 case 27: 01152 return 2.*xi*eta; 01153 01154 case 28: 01155 return 3.*eta*eta; 01156 01157 case 29: 01158 return 0.; 01159 01160 case 30: 01161 return 2.*xi*zeta; 01162 01163 case 31: 01164 return 4.*eta*zeta; 01165 01166 case 32: 01167 return 0.; 01168 01169 case 33: 01170 return 3.*zeta*zeta; 01171 01172 case 34: 01173 return 0.; 01174 01175 default: 01176 unsigned int o = 0; 01177 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 01178 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 01179 unsigned int block=o, nz = 0; 01180 for (; block < i2; block += (o-nz+1)) { nz++; } 01181 const unsigned int nx = block - i2; 01182 const unsigned int ny = o - nx - nz; 01183 Real val = ny * nz; 01184 for (unsigned int index=0; index != nx; index++) 01185 val *= xi; 01186 for (unsigned int index=1; index < ny; index++) 01187 val *= eta; 01188 for (unsigned int index=1; index < nz; index++) 01189 val *= zeta; 01190 return val; 01191 } 01192 } 01193 01194 01195 // d^2()/dzeta^2 01196 case 5: 01197 { 01198 switch (i) 01199 { 01200 // constant 01201 case 0: 01202 01203 // linear 01204 case 1: 01205 case 2: 01206 case 3: 01207 return 0.; 01208 01209 // quadratic 01210 case 4: 01211 case 5: 01212 case 6: 01213 case 7: 01214 case 8: 01215 return 0.; 01216 01217 case 9: 01218 return 2.; 01219 01220 // cubic 01221 case 10: 01222 case 11: 01223 case 12: 01224 case 13: 01225 case 14: 01226 case 15: 01227 case 16: 01228 return 0.; 01229 01230 case 17: 01231 return 2.*xi; 01232 01233 case 18: 01234 return 2.*eta; 01235 01236 case 19: 01237 return 6.*zeta; 01238 01239 // quartics 01240 case 20: 01241 case 21: 01242 case 22: 01243 case 23: 01244 case 24: 01245 case 25: 01246 case 26: 01247 case 27: 01248 case 28: 01249 return 0.; 01250 01251 case 29: 01252 return 2.*xi*xi; 01253 01254 case 30: 01255 return 2.*xi*eta; 01256 01257 case 31: 01258 return 2.*eta*eta; 01259 01260 case 32: 01261 return 6.*xi*zeta; 01262 01263 case 33: 01264 return 6.*eta*zeta; 01265 01266 case 34: 01267 return 12.*zeta*zeta; 01268 01269 default: 01270 unsigned int o = 0; 01271 for (; i >= (o+1)*(o+2)*(o+3)/6; o++) { } 01272 unsigned int i2 = i - (o*(o+1)*(o+2)/6); 01273 unsigned int block=o, nz = 0; 01274 for (; block < i2; block += (o-nz+1)) { nz++; } 01275 const unsigned int nx = block - i2; 01276 const unsigned int ny = o - nx - nz; 01277 Real val = nz * (nz - 1); 01278 for (unsigned int index=0; index != nx; index++) 01279 val *= xi; 01280 for (unsigned int index=0; index != ny; index++) 01281 val *= eta; 01282 for (unsigned int index=2; index < nz; index++) 01283 val *= zeta; 01284 return val; 01285 } 01286 } 01287 01288 01289 default: 01290 libmesh_error(); 01291 } 01292 01293 #endif 01294 01295 libmesh_error(); 01296 return 0.; 01297 } 01298 01299 01300 01301 template <> 01302 Real FE<3,MONOMIAL>::shape_second_deriv(const Elem* elem, 01303 const Order order, 01304 const unsigned int i, 01305 const unsigned int j, 01306 const Point& p) 01307 { 01308 libmesh_assert(elem); 01309 01310 // call the orientation-independent shape function derivatives 01311 return FE<3,MONOMIAL>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p); 01312 } 01313 01314 } // namespace libMesh
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:46 UTC
Hosted By: