inf_fe_lagrange_eval.c

一个用来实现偏微分方程中网格的计算库

// $Id: inf_fe_lagrange_eval.C 2789 2008-04-13 02:24:40Z roystgnr $

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

// 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_config.h"
#ifdef ENABLE_INFINITE_ELEMENTS
#include "inf_fe.h"



template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
Real InfFE<Dim,T_radial,T_map>::eval(const Real v,
				     const Order o_radial,
				     const unsigned int i)
{
  libmesh_assert (-1.-1.e-5 <= v && v < 1.);

  switch (o_radial)
    {
    case CONSTANT:
      switch (i)
        {
        case 0:
          return 1.;

        default:
          std::cerr << "bad index i = " << i << std::endl;
          libmesh_error();
        }


    case FIRST:
      switch (i)
        {
        case 0:
          return -1.*v;

        case 1:
          return v+1.;

        default:
          std::cerr << "bad index i = " << i << std::endl;
          libmesh_error();
        }


    case SECOND:
      switch (i)
        {
        case 0:
          return 2.*v*(v+.5);

        case 1:
          return 2.*(v+1.)*(v+.5);

        case 2:
          return -4.*(v+1.)*v;

        default:
          std::cerr << "bad index i = " << i << std::endl;
          libmesh_error();
        }


    case THIRD:
      switch (i)
        {
        case 0:
          return -4.5*v*(v+.6666666666666667)*(v+.3333333333333333);

        case 1:
          return 4.5*(v+1.)*(v+.6666666666666667)*(v+.3333333333333333);

        case 2:
          return 13.5*(v+1.)*v*(v+.3333333333333333);

        case 3:
          return -13.5*(v+1.)*v*(v+.6666666666666667);

        default:
          std::cerr << "bad index i = " << i << std::endl;
          libmesh_error();
        }


    case FOURTH:
      switch (i)
        {
        case 0:
          return 10.66666666666667*v*(v+.75)*(v+.5)*(v+.25);

        case 1:
          return 10.66666666666667*(v+1.)*(v+.75)*(v+.5)*(v+.25);

        case 2:
          return -42.66666666666667*(v+1.)*v*(v+.5)*(v+.25);

        case 3:
          return 64.*(v+1.)*v*(v+.75)*(v+.25);

        case 4:
          return -42.66666666666667*(v+1.)*v*(v+.75)*(v+.5);

        default:
          std::cerr << "bad index i = " << i << std::endl;
          libmesh_error();
        }


    case FIFTH:
      switch (i)
        {
        case 0:
          return -26.04166666666667*v*(v+.8)*(v+.6)*(v+.4)*(v+.2);

        case 1:
          return 26.04166666666667*(v+1.)*(v+.8)*(v+.6)*(v+.4)*(v+.2);

        case 2:
          return 130.2083333333333*(v+1.)*v*(v+.6)*(v+.4)*(v+.2);

        case 3:
          return -260.4166666666667*(v+1.)*v*(v+.8)*(v+.4)*(v+.2);

        case 4:
          return 260.4166666666667*(v+1.)*v*(v+.8)*(v+.6)*(v+.2);

        case 5:
          return -130.2083333333333*(v+1.)*v*(v+.8)*(v+.6)*(v+.4);

        default:
          std::cerr << "bad index i = " << i << std::endl;
          libmesh_error();
        }


    case SIXTH:
      switch (i)
        {
        case 0:
          return 64.8*v*(v+.8333333333333333)*(v+.6666666666666667)*(v+.5)*(v+.3333333333333333)*(v+.1666666666666667);

        case 1:
          return 64.8*(v+1.)*(v+.8333333333333333)*(v+.6666666666666667)*(v+.5)*(v+.3333333333333333)*(v+.1666666666666667);

        case 2:
          return -388.8*(v+1.)*v*(v+.6666666666666667)*(v+.5)*(v+.3333333333333333)*(v+.1666666666666667);

        case 3:
          return 972.*(v+1.)*v*(v+.8333333333333333)*(v+.5)*(v+.3333333333333333)*(v+.1666666666666667);

        case 4:
          return -1296.*(v+1.)*v*(v+.8333333333333333)*(v+.6666666666666667)*(v+.3333333333333333)*(v+.1666666666666667);

        case 5:
          return 972.*(v+1.)*v*(v+.8333333333333333)*(v+.6666666666666667)*(v+.5)*(v+.1666666666666667);

        case 6:
          return -388.8*(v+1.)*v*(v+.8333333333333333)*(v+.6666666666666667)*(v+.5)*(v+.3333333333333333);

        default:
          std::cerr << "bad index i = " << i << std::endl;
          libmesh_error();
        }


    case SEVENTH:
      switch (i)
        {
        case 0:
          return -163.401388888889*v*(v+.8571428571428571)*(v+.7142857142857143)*(v+.5714285714285714)*(v+.4285714285714286)*(v+.2857142857142857)*(v+.1428571428571429);

        case 1:
          return 163.4013888888888*(v+1.)*(v+.8571428571428571)*(v+.7142857142857143)*(v+.5714285714285714)*(v+.4285714285714286)*(v+.2857142857142857)*(v+.1428571428571429);

        case 2:
          return 1143.809722222222*(v+1.)*v*(v+.7142857142857143)*(v+.5714285714285714)*(v+.4285714285714286)*(v+.2857142857142857)*(v+.1428571428571429);

        case 3:
          return -3431.429166666665*(v+1.)*v*(v+.8571428571428571)*(v+.5714285714285714)*(v+.4285714285714286)*(v+.2857142857142857)*(v+.1428571428571429);

        case 4:
          return 5719.048611111112*(v+1.)*v*(v+.8571428571428571)*(v+.7142857142857143)*(v+.4285714285714286)*(v+.2857142857142857)*(v+.1428571428571429);

        case 5:
          return -5719.048611111112*(v+1.)*v*(v+.8571428571428571)*(v+.7142857142857143)*(v+.5714285714285714)*(v+.2857142857142857)*(v+.1428571428571429);

        case 6:
          return 3431.429166666666*(v+1.)*v*(v+.8571428571428571)*(v+.7142857142857143)*(v+.5714285714285714)*(v+.4285714285714286)*(v+.1428571428571429);

        case 7:
          return -1143.809722222223*(v+1.)*v*(v+.8571428571428571)*(v+.7142857142857143)*(v+.5714285714285714)*(v+.4285714285714286)*(v+.2857142857142857);

        default:
          std::cerr << "bad index i = " << i << std::endl;
          libmesh_error();
        }


    case EIGHTH:
      switch (i)
        {
        case 0:
          return 416.1015873015873*v*(v+.875)*(v+.75)*(v+.625)*(v+.5)*(v+.375)*(v+.25)*(v+.125);

        case 1:
          return 416.1015873015873*(v+1.)*(v+.875)*(v+.75)*(v+.625)*(v+.5)*(v+.375)*(v+.25)*(v+.125);

        case 2:
          return -3328.812698412698*(v+1.)*v*(v+.75)*(v+.625)*(v+.5)*(v+.375)*(v+.25)*(v+.125);

        case 3:
          return 11650.84444444444*(v+1.)*v*(v+.875)*(v+.625)*(v+.5)*(v+.375)*(v+.25)*(v+.125);

        case 4:
          return -23301.68888888889*(v+1.)*v*(v+.875)*(v+.75)*(v+.5)*(v+.375)*(v+.25)*(v+.125);

        case 5:
          return 29127.11111111111*(v+1.)*v*(v+.875)*(v+.75)*(v+.625)*(v+.375)*(v+.25)*(v+.125);

        case 6:
          return -23301.68888888889*(v+1.)*v*(v+.875)*(v+.75)*(v+.625)*(v+.5)*(v+.25)*(v+.125);

        case 7:
          return 11650.84444444444*(v+1.)*v*(v+.875)*(v+.75)*(v+.625)*(v+.5)*(v+.375)*(v+.125);

        case 8:
          return -3328.812698412698*(v+1.)*v*(v+.875)*(v+.75)*(v+.625)*(v+.5)*(v+.375)*(v+.25);

        default:
          std::cerr << "bad index i = " << i << std::endl;
          libmesh_error();
        }


    case NINTH:
      switch (i)
        {
        case 0:
          return -1067.627008928571*v*(v+.8888888888888889)*(v+.7777777777777778)*(v+.6666666666666667)*(v+.5555555555555556)*(v+.4444444444444444)*(v+.3333333333333333)*(v+.2222222222222222)*(v+.1111111111111111);

        case 1:
          return 1067.627008928571*(v+1.)*(v+.8888888888888889)*(v+.7777777777777778)*(v+.6666666666666667)*(v+.5555555555555556)*(v+.4444444444444444)*(v+.3333333333333333)*(v+.2222222222222222)*(v+.1111111111111111);

        case 2:
          return 9608.643080357158*(v+1.)*v*(v+.7777777777777778)*(v+.6666666666666667)*(v+.5555555555555556)*(v+.4444444444444444)*(v+.3333333333333333)*(v+.2222222222222222)*(v+.1111111111111111);

        case 3:
          return -38434.57232142862*(v+1.)*v*(v+.8888888888888889)*(v+.6666666666666667)*(v+.5555555555555556)*(v+.4444444444444444)*(v+.3333333333333333)*(v+.2222222222222222)*(v+.1111111111111111);

        case 4:
          return 89680.66874999999*(v+1.)*v*(v+.8888888888888889)*(v+.7777777777777778)*(v+.5555555555555556)*(v+.4444444444444444)*(v+.3333333333333333)*(v+.2222222222222222)*(v+.1111111111111111);

        case 5:
          return -134521.0031249999*(v+1.)*v*(v+.8888888888888889)*(v+.7777777777777778)*(v+.6666666666666667)*(v+.4444444444444444)*(v+.3333333333333333)*(v+.2222222222222222)*(v+.1111111111111111);