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