fe_bernstein_shape_1D.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 
00020 // Local includes
00021 #include "libmesh/libmesh_config.h"
00022 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
00023 
00024 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
00025 #include <cmath>
00026 
00027 #include "libmesh/libmesh_common.h"
00028 #include "libmesh/fe.h"
00029 #include "libmesh/elem.h"
00030 #include "libmesh/utility.h"
00031 
00032 
00033 namespace libMesh
00034 {
00035 
00036 
00037 template <>
00038 Real FE<1,BERNSTEIN>::shape(const ElemType,
00039                             const Order order,
00040                             const unsigned int i,
00041                             const Point& p)
00042 {
00043   const Real xi = p(0);
00044   using Utility::pow;
00045 
00046   switch (order)
00047     {
00048     case FIRST:
00049 
00050       switch(i)
00051         {
00052         case 0:
00053           return (1.-xi)/2.;
00054         case 1:
00055           return (1.+xi)/2.;
00056         default:
00057           libMesh::err << "Invalid shape function index!" << std::endl;
00058           libmesh_error();
00059         }
00060 
00061     case SECOND:
00062 
00063       switch(i)
00064         {
00065         case 0:
00066           return (1./4.)*pow<2>(1.-xi);
00067         case 1:
00068           return (1./4.)*pow<2>(1.+xi);
00069         case 2:
00070           return (1./2.)*(1.-xi)*(1.+xi);
00071         default:
00072           libMesh::err << "Invalid shape function index!" << std::endl;
00073           libmesh_error();
00074         }
00075 
00076     case THIRD:
00077 
00078       switch(i)
00079         {
00080         case 0:
00081           return (1./8.)*pow<3>(1.-xi);
00082         case 1:
00083           return (1./8.)*pow<3>(1.+xi);
00084         case 2:
00085           return (3./8.)*(1.+xi)*pow<2>(1.-xi);
00086         case 3:
00087           return (3./8.)*pow<2>(1.+xi)*(1.-xi);
00088         default:
00089           libMesh::err << "Invalid shape function index!" << std::endl;
00090           libmesh_error();
00091         }
00092 
00093     case FOURTH:
00094 
00095       switch(i)
00096         {
00097         case 0:
00098           return (1./16.)*pow<4>(1.-xi);
00099         case 1:
00100           return (1./16.)*pow<4>(1.+xi);
00101         case 2:
00102           return (1./ 4.)*(1.+xi)*pow<3>(1.-xi);
00103         case 3:
00104           return (3./ 8.)*pow<2>(1.+xi)*pow<2>(1.-xi);
00105         case 4:
00106           return (1./ 4.)*pow<3>(1.+xi)*(1.-xi);
00107         default:
00108           libMesh::err << "Invalid shape function index!" << std::endl;
00109           libmesh_error();
00110         }
00111 
00112 
00113     case FIFTH:
00114 
00115       switch(i)
00116         {
00117         case 0:
00118           return (1./32.)*pow<5>(1.-xi);
00119         case 1:
00120           return (1./32.)*pow<5>(1.+xi);
00121         case 2:
00122           return (5./32.)*(1.+xi)*pow<4>(1.-xi);
00123         case 3:
00124           return (5./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
00125         case 4:
00126           return (5./16.)*pow<3>(1.+xi)*pow<2>(1.-xi);
00127         case 5:
00128           return (5./32.)*pow<4>(1.+xi)*(1.-xi);
00129         default:
00130           libMesh::err << "Invalid shape function index!" << std::endl;
00131           libmesh_error();
00132         }
00133 
00134 
00135     case SIXTH:
00136 
00137       switch (i)
00138         {
00139         case 0:
00140           return ( 1./64.)*pow<6>(1.-xi);
00141         case 1:
00142           return ( 1./64.)*pow<6>(1.+xi);
00143         case 2:
00144           return ( 3./32.)*(1.+xi)*pow<5>(1.-xi);
00145         case 3:
00146           return (15./64.)*pow<2>(1.+xi)*pow<4>(1.-xi);
00147         case 4:
00148           return ( 5./16.)*pow<3>(1.+xi)*pow<3>(1.-xi);
00149         case 5:
00150           return (15./64.)*pow<4>(1.+xi)*pow<2>(1.-xi);
00151         case 6:
00152           return ( 3./32.)*pow<5>(1.+xi)*(1.-xi);
00153         default:
00154           libMesh::err << "Invalid shape function index!" << std::endl;
00155           libmesh_error();
00156         }
00157 
00158     default:
00159       {
00160         libmesh_assert (order>6);
00161 
00162         // Use this for arbitrary orders.
00163         // Note that this implementation is less efficient.
00164         const int p_order = static_cast<unsigned int>(order);
00165         const int m       = p_order-i+1;
00166         const int n       = (i-1);
00167 
00168         unsigned int binomial_p_i = 1;
00169 
00170         // the binomial coefficient (p choose n)
00171         if (i>1)
00172           binomial_p_i = Utility::factorial(p_order)
00173             / (Utility::factorial(n)*Utility::factorial(p_order-n));
00174 
00175 
00176         switch(i)
00177           {
00178           case 0:
00179             return binomial_p_i * std::pow((1-xi)/2,static_cast<int>(p_order));
00180           case 1:
00181             return binomial_p_i * std::pow((1+xi)/2,static_cast<int>(p_order));
00182           default:
00183             {
00184               return binomial_p_i * std::pow((1+xi)/2,n)
00185                                   * std::pow((1-xi)/2,m);
00186             }
00187           }
00188 
00189         // we should never get here
00190         libmesh_error();
00191       }
00192     }
00193 
00194   libmesh_error();
00195   return 0.;
00196 }
00197 
00198 
00199 
00200 template <>
00201 Real FE<1,BERNSTEIN>::shape(const Elem* elem,
00202                             const Order order,
00203                             const unsigned int i,
00204                             const Point& p)
00205 {
00206   libmesh_assert(elem);
00207 
00208   return FE<1,BERNSTEIN>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p);
00209 }
00210 
00211 
00212 
00213 template <>
00214 Real FE<1,BERNSTEIN>::shape_deriv(const ElemType,
00215                                   const Order order,
00216                                   const unsigned int i,
00217                                   const unsigned int libmesh_dbg_var(j),
00218                                   const Point& p)
00219 {
00220   // only d()/dxi in 1D!
00221 
00222   libmesh_assert_equal_to (j, 0);
00223 
00224   const Real xi = p(0);
00225 
00226   using Utility::pow;
00227 
00228   switch (order)
00229     {
00230     case FIRST:
00231 
00232       switch(i)
00233         {
00234         case 0:
00235           return -.5;
00236         case 1:
00237           return .5;
00238         default:
00239           libMesh::err << "Invalid shape function index " << i << std::endl;
00240           libmesh_error();
00241         }
00242 
00243     case SECOND:
00244 
00245       switch(i)
00246         {
00247         case 0:
00248           return (xi-1.)*.5;
00249         case 1:
00250           return (xi+1.)*.5;
00251         case 2:
00252           return -xi;
00253         default:
00254           libMesh::err << "Invalid shape function index!" << std::endl;
00255           libmesh_error();
00256         }
00257 
00258     case THIRD:
00259 
00260       switch(i)
00261         {
00262         case 0:
00263           return -0.375*pow<2>(1.-xi);
00264         case 1:
00265           return  0.375*pow<2>(1.+xi);
00266         case 2:
00267           return -0.375 -.75*xi +1.125*pow<2>(xi);
00268         case 3:
00269           return  0.375 -.75*xi -1.125*pow<2>(xi);
00270         default:
00271           libMesh::err << "Invalid shape function index!" << std::endl;
00272           libmesh_error();
00273         }
00274 
00275     case FOURTH:
00276 
00277       switch(i)
00278         {
00279         case 0:
00280           return -0.25*pow<3>(1.-xi);
00281         case 1:
00282           return  0.25*pow<3>(1.+xi);
00283         case 2:
00284           return -0.5 +1.5*pow<2>(xi)-pow<3>(xi);
00285         case 3:
00286           return  1.5*(pow<3>(xi)-xi);
00287         case 4:
00288           return  0.5 -1.5*pow<2>(xi)-pow<3>(xi);
00289         default:
00290           libMesh::err << "Invalid shape function index!" << std::endl;
00291           libmesh_error();
00292         }
00293 
00294     case FIFTH:
00295 
00296       switch(i)
00297         {
00298         case 0:
00299           return -(5./32.)*pow<4>(xi-1.);
00300         case 1:
00301           return  (5./32.)*pow<4>(xi+1.);
00302         case 2:
00303           return  (5./32.)*pow<4>(1.-xi)         -(5./8.)*(1.+xi)*pow<3>(1.-xi);
00304         case 3:
00305           return  (5./ 8.)*(1.+xi)*pow<3>(1.-xi) -(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
00306         case 4:
00307           return -(5./ 8.)*pow<3>(1.+xi)*(1.-xi) +(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
00308         case 5:
00309           return  (5./ 8.)*pow<3>(1.+xi)*(1.-xi) -(5./32.)*pow<4>(1.+xi);
00310         default:
00311           libMesh::err << "Invalid shape function index!" << std::endl;
00312           libmesh_error();
00313         }
00314 
00315     case SIXTH:
00316 
00317       switch(i)
00318         {
00319         case 0:
00320           return -( 3./32.)*pow<5>(1.-xi);
00321         case 1:
00322           return  ( 3./32.)*pow<5>(1.+xi);
00323         case 2:
00324           return  ( 3./32.)*pow<5>(1.-xi)-(15./32.)*(1.+xi)*pow<4>(1.-xi);
00325         case 3:
00326           return  (15./32.)*(1.+xi)*pow<4>(1.-xi)-(15./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
00327         case 4:
00328           return -(15./ 8.)*xi +(15./4.)*pow<3>(xi)-(15./8.)*pow<5>(xi);
00329         case 5:
00330           return -(15./32.)*(1.-xi)*pow<4>(1.+xi)+(15./16.)*pow<2>(1.-xi)*pow<3>(1.+xi);
00331         case 6:
00332           return  (15./32.)*pow<4>(1.+xi)*(1.-xi)-(3./32.)*pow<5>(1.+xi);
00333         default:
00334           libMesh::err << "Invalid shape function index!" << std::endl;
00335           libmesh_error();
00336         }
00337 
00338 
00339     default:
00340       {
00341         libmesh_assert (order>6);
00342 
00343         // Use this for arbitrary orders
00344         const int p_order = static_cast<unsigned int>(order);
00345         const int m       = p_order-(i-1);
00346         const int n       = (i-1);
00347 
00348         unsigned int binomial_p_i = 1;
00349 
00350         // the binomial coefficient (p choose n)
00351         if (i>1)
00352           binomial_p_i = Utility::factorial(p_order)
00353             / (Utility::factorial(n)*Utility::factorial(p_order-n));
00354 
00355 
00356 
00357         switch(i)
00358           {
00359           case 0:
00360             return binomial_p_i * (-1./2.) * p_order * std::pow((1-xi)/2,static_cast<int>(p_order-1));
00361           case 1:
00362             return binomial_p_i * ( 1./2.) * p_order * std::pow((1+xi)/2,static_cast<int>(p_order-1));
00363 
00364           default:
00365             {
00366               return binomial_p_i * (1./2. * n * std::pow((1+xi)/2,n-1) * std::pow((1-xi)/2,m)
00367                                    - 1./2. * m * std::pow((1+xi)/2,n)   * std::pow((1-xi)/2,m-1));
00368             }
00369           }
00370 
00371         // we should never get here
00372         libmesh_error();
00373       }
00374 
00375     }
00376 
00377   libmesh_error();
00378   return 0.;
00379 }
00380 
00381 
00382 
00383 template <>
00384 Real FE<1,BERNSTEIN>::shape_deriv(const Elem* elem,
00385                                   const Order order,
00386                                   const unsigned int i,
00387                                   const unsigned int j,
00388                                   const Point& p)
00389 {
00390   libmesh_assert(elem);
00391 
00392   return FE<1,BERNSTEIN>::shape_deriv(elem->type(),
00393                                       static_cast<Order>(order + elem->p_level()), i, j, p);
00394 }
00395 
00396 
00397 
00398 template <>
00399 Real FE<1,BERNSTEIN>::shape_second_deriv(const ElemType,
00400                                          const Order,
00401                                          const unsigned int,
00402                                          const unsigned int,
00403                                          const Point&)
00404 {
00405   static bool warning_given = false;
00406 
00407   if (!warning_given)
00408   libMesh::err << "Second derivatives for Bernstein elements "
00409                 << "are not yet implemented!"
00410                 << std::endl;
00411 
00412   warning_given = true;
00413   return 0.;
00414 }
00415 
00416 
00417 
00418 
00419 template <>
00420 Real FE<1,BERNSTEIN>::shape_second_deriv(const Elem*,
00421                                          const Order,
00422                                          const unsigned int,
00423                                          const unsigned int,
00424                                          const Point&)
00425 {
00426   static bool warning_given = false;
00427 
00428   if (!warning_given)
00429   libMesh::err << "Second derivatives for Bernstein elements "
00430                 << "are not yet implemented!"
00431                 << std::endl;
00432 
00433   warning_given = true;
00434   return 0.;
00435 }
00436 
00437 } // namespace libMesh
00438 
00439 
00440 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:46 UTC

Hosted By: