Home
Publications
Presentations
Support
Download
Examples
Developers
Class Docs
Namespaces
Classes
Files
Mailing Lists
Gallery
Wiki

fe_szabab_shape_1D.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2012 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



// C++ includes

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

#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES

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

namespace libMesh
{


template <>
Real FE<1,SZABAB>::shape(const ElemType,
                         const Order libmesh_dbg_var(order),
                         const unsigned int i,
                         const Point& p)
{
  const Real xi  = p(0);
  const Real xi2 = xi*xi;


  // Use this libmesh_assert rather than a switch with a single entry...
  // It will go away in optimized mode, essentially has the same effect.
  libmesh_assert_less_equal (order, SEVENTH);

//   switch (order)
//     {
//     case FIRST:
//     case SECOND:
//     case THIRD:
//     case FOURTH:
//     case FIFTH:
//     case SIXTH:
//     case SEVENTH:

      switch(i)
        {
          //nodal shape functions
        case 0: return 1./2.-1./2.*xi;
        case 1: return 1./2.+1./2.*xi;
        case 2: return 1./4.  *2.4494897427831780982*(xi2-1.);
        case 3: return 1./4.  *3.1622776601683793320*(xi2-1.)*xi;
        case 4: return 1./16. *3.7416573867739413856*((5.*xi2-6.)*xi2+1.);
        case 5: return 3./16. *1.4142135623730950488*(3.+(-10.+7.*xi2)*xi2)*xi;
        case 6: return 1./32. *4.6904157598234295546*(-1.+(15.+(-35.+21.*xi2)*xi2)*xi2);
        case 7: return 1./32. *5.0990195135927848300*(-5.+(35.+(-63.+33.*xi2)*xi2)*xi2)*xi;
        case 8: return 1./256.*5.4772255750516611346*(5.+(-140.+(630.+(-924.+429.*xi2)*xi2)*xi2)*xi2);

        default:
          libMesh::err << "Invalid shape function index!" << std::endl;
          libmesh_error();
        }

//     default:
//       {
//      libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
//      libmesh_error();
//       }
//     }

  libmesh_error();
  return 0.;
}



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

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



template <>
Real FE<1,SZABAB>::shape_deriv(const ElemType,
                               const Order libmesh_dbg_var(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);
  const Real xi2 = xi*xi;

  // Use this libmesh_assert rather than a switch with a single entry...
  // It will go away in optimized mode, essentially has the same effect.
  libmesh_assert_less_equal (order, SEVENTH);

//   switch (order)
//     {
//     case FIRST:
//     case SECOND:
//     case THIRD:
//     case FOURTH:
//     case FIFTH:
//     case SIXTH:
//     case SEVENTH:

      switch(i)
        {
        case 0: return -1./2.;
        case 1: return 1./2.;
        case 2: return 1./2.*2.4494897427831780982*xi;
        case 3: return -1./4.*3.1622776601683793320+3./4.*3.1622776601683793320*xi2;
        case 4: return 1./16.*3.7416573867739413856*(-12.+20*xi2)*xi;
        case 5: return 9./16.*1.4142135623730950488+(-45./8.*1.4142135623730950488+105./16.*1.4142135623730950488*xi2)*xi2;
        case 6: return 1./32.*4.6904157598234295546*(30.+(-140.+126.*xi2)*xi2)*xi;
        case 7: return -5./32.*5.0990195135927848300+(105./32.*5.0990195135927848300+(-315./32.*5.0990195135927848300+231./32.*5.0990195135927848300*xi2)*xi2)*xi2;
        case 8: return 1./256.*5.4772255750516611346*(-280.+(2520.+(-5544.+3432.*xi2)*xi2)*xi2)*xi;

        default:
          libMesh::err << "Invalid shape function index!" << std::endl;
          libmesh_error();
        }

//     default:
//       {
//      libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
//      libmesh_error();
//       }
//     }

  libmesh_error();
  return 0.;
}



template <>
Real FE<1,SZABAB>::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,SZABAB>::shape_deriv(elem->type(),
                                       static_cast<Order>(order + elem->p_level()), i, j, p);
}



template <>
Real FE<1,SZABAB>::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 Szabab elements "
            << " are not yet implemented!"
            << std::endl;

  warning_given = true;
  return 0.;
}



template <>
Real FE<1,SZABAB>::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 Szabab elements "
            << " are not yet implemented!"
            << std::endl;

  warning_given = true;
  return 0.;
}

} // namespace libMesh


#endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES

---

Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:46 UTC

Hosted By: