fe_bernstein_shape_1D.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



// Local includes
#include "libmesh/libmesh_config.h"
#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES

#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
#include <cmath>

#include "libmesh/libmesh_common.h"
#include "libmesh/fe.h"
#include "libmesh/elem.h"
#include "libmesh/utility.h"


namespace libMesh
{


template <>
Real FE<1,BERNSTEIN>::shape(const ElemType,
                            const Order order,
                            const unsigned int i,
                            const Point& p)
{
  const Real xi = p(0);
  using Utility::pow;

  switch (order)
    {
    case FIRST:

      switch(i)
        {
        case 0:
          return (1.-xi)/2.;
        case 1:
          return (1.+xi)/2.;
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }

    case SECOND:

      switch(i)
        {
        case 0:
          return (1./4.)*pow<2>(1.-xi);
        case 1:
          return (1./4.)*pow<2>(1.+xi);
        case 2:
          return (1./2.)*(1.-xi)*(1.+xi);
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }

    case THIRD:

      switch(i)
        {
        case 0:
          return (1./8.)*pow<3>(1.-xi);
        case 1:
          return (1./8.)*pow<3>(1.+xi);
        case 2:
          return (3./8.)*(1.+xi)*pow<2>(1.-xi);
        case 3:
          return (3./8.)*pow<2>(1.+xi)*(1.-xi);
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }

    case FOURTH:

      switch(i)
        {
        case 0:
          return (1./16.)*pow<4>(1.-xi);
        case 1:
          return (1./16.)*pow<4>(1.+xi);
        case 2:
          return (1./ 4.)*(1.+xi)*pow<3>(1.-xi);
        case 3:
          return (3./ 8.)*pow<2>(1.+xi)*pow<2>(1.-xi);
        case 4:
          return (1./ 4.)*pow<3>(1.+xi)*(1.-xi);
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }


    case FIFTH:

      switch(i)
        {
        case 0:
          return (1./32.)*pow<5>(1.-xi);
        case 1:
          return (1./32.)*pow<5>(1.+xi);
        case 2:
          return (5./32.)*(1.+xi)*pow<4>(1.-xi);
        case 3:
          return (5./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
        case 4:
          return (5./16.)*pow<3>(1.+xi)*pow<2>(1.-xi);
        case 5:
          return (5./32.)*pow<4>(1.+xi)*(1.-xi);
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }


    case SIXTH:

      switch (i)
        {
        case 0:
          return ( 1./64.)*pow<6>(1.-xi);
        case 1:
          return ( 1./64.)*pow<6>(1.+xi);
        case 2:
          return ( 3./32.)*(1.+xi)*pow<5>(1.-xi);
        case 3:
          return (15./64.)*pow<2>(1.+xi)*pow<4>(1.-xi);
        case 4:
          return ( 5./16.)*pow<3>(1.+xi)*pow<3>(1.-xi);
        case 5:
          return (15./64.)*pow<4>(1.+xi)*pow<2>(1.-xi);
        case 6:
          return ( 3./32.)*pow<5>(1.+xi)*(1.-xi);
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }

    default:
      {
        libmesh_assert (order>6);

        // Use this for arbitrary orders.
        // Note that this implementation is less efficient.
        const int p_order = static_cast<unsigned int>(order);
        const int m       = p_order-i+1;
        const int n       = (i-1);

        unsigned int binomial_p_i = 1;

        // the binomial coefficient (p choose n)
        if (i>1)
          binomial_p_i = Utility::factorial(p_order)
            / (Utility::factorial(n)*Utility::factorial(p_order-n));


        switch(i)
          {
          case 0:
            return binomial_p_i * std::pow((1-xi)/2,static_cast<int>(p_order));
          case 1:
            return binomial_p_i * std::pow((1+xi)/2,static_cast<int>(p_order));
          default:
            {
              return binomial_p_i * std::pow((1+xi)/2,n)
                * std::pow((1-xi)/2,m);
            }
          }
      }
    }

  libmesh_error_msg("We'll never get here!");
  return 0.;
}



template <>
Real FE<1,BERNSTEIN>::shape(const Elem* elem,
                            const Order order,
                            const unsigned int i,
                            const Point& p)
{
  libmesh_assert(elem);

  return FE<1,BERNSTEIN>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p);
}



template <>
Real FE<1,BERNSTEIN>::shape_deriv(const ElemType,
                                  const Order order,
                                  const unsigned int i,
                                  const unsigned int libmesh_dbg_var(j),
                                  const Point& p)
{
  // only d()/dxi in 1D!

  libmesh_assert_equal_to (j, 0);

  const Real xi = p(0);

  using Utility::pow;

  switch (order)
    {
    case FIRST:

      switch(i)
        {
        case 0:
          return -.5;
        case 1:
          return .5;
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }

    case SECOND:

      switch(i)
        {
        case 0:
          return (xi-1.)*.5;
        case 1:
          return (xi+1.)*.5;
        case 2:
          return -xi;
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }

    case THIRD:

      switch(i)
        {
        case 0:
          return -0.375*pow<2>(1.-xi);
        case 1:
          return  0.375*pow<2>(1.+xi);
        case 2:
          return -0.375 -.75*xi +1.125*pow<2>(xi);
        case 3:
          return  0.375 -.75*xi -1.125*pow<2>(xi);
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }

    case FOURTH:

      switch(i)
        {
        case 0:
          return -0.25*pow<3>(1.-xi);
        case 1:
          return  0.25*pow<3>(1.+xi);
        case 2:
          return -0.5 +1.5*pow<2>(xi)-pow<3>(xi);
        case 3:
          return  1.5*(pow<3>(xi)-xi);
        case 4:
          return  0.5 -1.5*pow<2>(xi)-pow<3>(xi);
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }

    case FIFTH:

      switch(i)
        {
        case 0:
          return -(5./32.)*pow<4>(xi-1.);
        case 1:
          return  (5./32.)*pow<4>(xi+1.);
        case 2:
          return  (5./32.)*pow<4>(1.-xi)         -(5./8.)*(1.+xi)*pow<3>(1.-xi);
        case 3:
          return  (5./ 8.)*(1.+xi)*pow<3>(1.-xi) -(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
        case 4:
          return -(5./ 8.)*pow<3>(1.+xi)*(1.-xi) +(15./16.)*pow<2>(1.+xi)*pow<2>(1.-xi);
        case 5:
          return  (5./ 8.)*pow<3>(1.+xi)*(1.-xi) -(5./32.)*pow<4>(1.+xi);
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }

    case SIXTH:

      switch(i)
        {
        case 0:
          return -( 3./32.)*pow<5>(1.-xi);
        case 1:
          return  ( 3./32.)*pow<5>(1.+xi);
        case 2:
          return  ( 3./32.)*pow<5>(1.-xi)-(15./32.)*(1.+xi)*pow<4>(1.-xi);
        case 3:
          return  (15./32.)*(1.+xi)*pow<4>(1.-xi)-(15./16.)*pow<2>(1.+xi)*pow<3>(1.-xi);
        case 4:
          return -(15./ 8.)*xi +(15./4.)*pow<3>(xi)-(15./8.)*pow<5>(xi);
        case 5:
          return -(15./32.)*(1.-xi)*pow<4>(1.+xi)+(15./16.)*pow<2>(1.-xi)*pow<3>(1.+xi);
        case 6:
          return  (15./32.)*pow<4>(1.+xi)*(1.-xi)-(3./32.)*pow<5>(1.+xi);
        default:
          libmesh_error_msg("Invalid shape function index i = " << i);
        }


    default:
      {
        libmesh_assert (order>6);

        // Use this for arbitrary orders
        const int p_order = static_cast<unsigned int>(order);
        const int m       = p_order-(i-1);
        const int n       = (i-1);

        unsigned int binomial_p_i = 1;

        // the binomial coefficient (p choose n)
        if (i>1)
          binomial_p_i = Utility::factorial(p_order)
            / (Utility::factorial(n)*Utility::factorial(p_order-n));



        switch(i)
          {
          case 0:
            return binomial_p_i * (-1./2.) * p_order * std::pow((1-xi)/2,static_cast<int>(p_order-1));
          case 1:
            return binomial_p_i * ( 1./2.) * p_order * std::pow((1+xi)/2,static_cast<int>(p_order-1));

          default:
            {
              return binomial_p_i * (1./2. * n * std::pow((1+xi)/2,n-1) * std::pow((1-xi)/2,m)
                                     - 1./2. * m * std::pow((1+xi)/2,n)   * std::pow((1-xi)/2,m-1));
            }
          }
      }

    }

  libmesh_error_msg("We'll never get here!");
  return 0.;
}



template <>
Real FE<1,BERNSTEIN>::shape_deriv(const Elem* elem,
                                  const Order order,
                                  const unsigned int i,
                                  const unsigned int j,
                                  const Point& p)
{
  libmesh_assert(elem);

  return FE<1,BERNSTEIN>::shape_deriv(elem->type(),
                                      static_cast<Order>(order + elem->p_level()), i, j, p);
}



template <>
Real FE<1,BERNSTEIN>::shape_second_deriv(const ElemType,
                                         const Order,
                                         const unsigned int,
                                         const unsigned int,
                                         const Point&)
{
  static bool warning_given = false;

  if (!warning_given)
    libMesh::err << "Second derivatives for Bernstein elements "
                 << "are not yet implemented!"
                 << std::endl;

  warning_given = true;
  return 0.;
}




template <>
Real FE<1,BERNSTEIN>::shape_second_deriv(const Elem*,
                                         const Order,
                                         const unsigned int,
                                         const unsigned int,
                                         const Point&)
{
  static bool warning_given = false;

  if (!warning_given)
    libMesh::err << "Second derivatives for Bernstein elements "
                 << "are not yet implemented!"
                 << std::endl;

  warning_given = true;
  return 0.;
}

} // namespace libMesh


#endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES