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