fe_hierarchic_shape_2d.c

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

// $Id: fe_hierarchic_shape_2D.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


// C++ includes

// Local includes
#include "fe.h"
#include "elem.h"
#include "number_lookups.h"
#include "utility.h"


template <>
Real FE<2,HIERARCHIC>::shape(const ElemType,
			     const Order,
			     const unsigned int,
			     const Point&)
{
  std::cerr << "Hierarchic polynomials require the element type\n"
	    << "because edge orientation is needed."
	    << std::endl;
  
  libmesh_error();
  return 0.;
}



template <>
Real FE<2,HIERARCHIC>::shape(const Elem* elem,
			     const Order order,
			     const unsigned int i,
			     const Point& p)
{
  libmesh_assert (elem != NULL);

  // Declare that we are using our own special power function
  // from the Utility namespace.  This saves typing later.
  using Utility::pow;

  const Order totalorder = static_cast<Order>(order+elem->p_level());
  libmesh_assert(totalorder > 0);
  
  switch (elem->type())
    {      
    case TRI3:
    case TRI6:
      {
	const Real zeta1 = p(0);
	const Real zeta2 = p(1);
	const Real zeta0 = 1. - zeta1 - zeta2;
	      
	libmesh_assert (i<(totalorder+1u)*(totalorder+2u)/2);
	libmesh_assert (elem->type() == TRI6 || totalorder < 2);

        // Vertex DoFs
        if (i == 0)
          return zeta0;
        else if (i == 1)
          return zeta1;
        else if (i == 2)
          return zeta2;
        // Edge DoFs
        else if (i < totalorder + 2u)
          {
            // Avoid returning NaN on vertices!
            if (zeta0 + zeta1 == 0.)
              return 0.;

            const unsigned int basisorder = i - 1;
	    // Get factors to account for edge-flipping
	    Real f0 = 1;
	    if (basisorder%2 && (elem->point(0) > elem->point(1)))
	      f0 = -1.;
	      
            Real edgeval = (zeta1 - zeta0) / (zeta1 + zeta0);
            Real crossfunc = zeta0 + zeta1;
            for (unsigned int n=1; n != basisorder; ++n)
              crossfunc *= (zeta0 + zeta1);

            return f0 * crossfunc *
                   FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
                                           basisorder, edgeval);
          }
        else if (i < 2u*totalorder + 1)
          {
            // Avoid returning NaN on vertices!
            if (zeta1 + zeta2 == 0.)
              return 0.;

            const unsigned int basisorder = i - totalorder;
	    // Get factors to account for edge-flipping
	    Real f1 = 1;
	    if (basisorder%2 && (elem->point(1) > elem->point(2)))
	      f1 = -1.;
	      
            Real edgeval = (zeta2 - zeta1) / (zeta2 + zeta1);
            Real crossfunc = zeta2 + zeta1;
            for (unsigned int n=1; n != basisorder; ++n)
              crossfunc *= (zeta2 + zeta1);

            return f1 * crossfunc *
                   FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
                                           basisorder, edgeval);
          }
        else if (i < 3u*totalorder)
          {
            // Avoid returning NaN on vertices!
            if (zeta0 + zeta2 == 0.)
              return 0.;

            const unsigned int basisorder = i - (2u*totalorder) + 1;
	    // Get factors to account for edge-flipping
	    Real f2 = 1;
	    if (basisorder%2 && (elem->point(2) > elem->point(0)))
	      f2 = -1.;

            Real edgeval = (zeta0 - zeta2) / (zeta0 + zeta2);
            Real crossfunc = zeta0 + zeta2;
            for (unsigned int n=1; n != basisorder; ++n)
              crossfunc *= (zeta0 + zeta2);

            return f2 * crossfunc *
                   FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
                                           basisorder, edgeval);
          }
        // Interior DoFs
        else
          {
            const unsigned int basisnum = i - (3u*totalorder);
            unsigned int exp0 = triangular_number_column[basisnum] + 1;
            unsigned int exp1 = triangular_number_row[basisnum] + 1 -
                                triangular_number_column[basisnum];

            Real returnval = 1;
            for (unsigned int n = 0; n != exp0; ++n)
              returnval *= zeta0;
            for (unsigned int n = 0; n != exp1; ++n)
              returnval *= zeta1;
            returnval *= zeta2;
            return returnval;
          }
      }

    // Hierarchic shape functions on the quadrilateral.
    case QUAD4:
      libmesh_assert (totalorder < 2);
    case QUAD8:
    case QUAD9:
      {
        // Compute quad shape functions as a tensor-product
        const Real xi  = p(0);
        const Real eta = p(1);
      
        libmesh_assert (i < (totalorder+1u)*(totalorder+1u));

// Example i, i0, i1 values for totalorder = 5:
//                                    0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
//  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 3, 2, 4, 4, 4, 3, 2, 5, 5, 5, 5, 4, 3, 2};
//  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 3, 3, 2, 3, 4, 4, 4, 2, 3, 4, 5, 5, 5, 5};
	      
        unsigned int i0, i1;

        // Vertex DoFs
        if (i == 0)
          { i0 = 0; i1 = 0; }
        else if (i == 1)
          { i0 = 1; i1 = 0; }
        else if (i == 2)
          { i0 = 1; i1 = 1; }
        else if (i == 3)
          { i0 = 0; i1 = 1; }
        // Edge DoFs
        else if (i < totalorder + 3u)
          { i0 = i - 2; i1 = 0; }
        else if (i < 2u*totalorder + 2)
          { i0 = 1; i1 = i - totalorder - 1; }
        else if (i < 3u*totalorder + 1)
          { i0 = i - 2u*totalorder; i1 = 1; }
        else if (i < 4u*totalorder)
          { i0 = 0; i1 = i - 3u*totalorder + 1; }
        // Interior DoFs
        else
          {
            unsigned int basisnum = i - 4*totalorder;
            i0 = square_number_column[basisnum] + 2;
            i1 = square_number_row[basisnum] + 2;
          }

        // Flip odd degree of freedom values if necessary
        // to keep continuity on sides
        Real f = 1.;

        if ((i0%2) && (i0 > 2) && (i1 == 0))
	  f = (elem->point(0) > elem->point(1))?-1.:1.;
        else if ((i0%2) && (i0>2) && (i1 == 1))
	  f = (elem->point(3) > elem->point(2))?-1.:1.;
        else if ((i0 == 0) && (i1%2) && (i1>2))
	  f = (elem->point(0) > elem->point(3))?-1.:1.;
        else if ((i0 == 1) && (i1%2) && (i1>2))
	  f = (elem->point(1) > elem->point(2))?-1.:1.;
	      
        return f*(FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i0, xi)*
		  FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i1, eta));	      
      }

    default:
      std::cerr << "ERROR: Unsupported element type!" << std::endl;
      libmesh_error();
    }

  return 0.;
}



template <>
Real FE<2,HIERARCHIC>::shape_deriv(const ElemType,
				   const Order,			    
				   const unsigned int,
				   const unsigned int,
				   const Point&)
{
  std::cerr << "Hierarchic polynomials require the element type\n"
	    << "because edge orientation is needed."
	    << std::endl;

  libmesh_error();
  return 0.;
}



template <>
Real FE<2,HIERARCHIC>::shape_deriv(const Elem* elem,
				   const Order order,
				   const unsigned int i,
				   const unsigned int j,
				   const Point& p)
{
  libmesh_assert (elem != NULL);

  const ElemType type = elem->type();

  const Order totalorder = static_cast<Order>(order+elem->p_level());

  libmesh_assert (totalorder > 0);
  
  switch (type)
    {
      // 1st & 2nd-order Hierarchics.
    case TRI3:
    case TRI6:
      {
	const Real eps = 1.e-6;
	      
	libmesh_assert (j < 2);
	      
	switch (j)
	  {
	    //  d()/dxi
	  case 0:
	    {
	      const Point pp(p(0)+eps, p(1));
	      const Point pm(p(0)-eps, p(1));

	      return (FE<2,HIERARCHIC>::shape(elem, order, i, pp) -
		      FE<2,HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
	    }

	     // d()/deta
	  case 1:
	    {
	      const Point pp(p(0), p(1)+eps);
	      const Point pm(p(0), p(1)-eps);

	      return (FE<2,HIERARCHIC>::shape(elem, order, i, pp) -
		      FE<2,HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
	    }
		  

	  default:
	    libmesh_error();
	  }
      }

    case QUAD4:
      libmesh_assert (totalorder < 2);
    case QUAD8:
    case QUAD9:
      {
	// Compute quad shape functions as a tensor-product
	const Real xi  = p(0);
	const Real eta = p(1);

        libmesh_assert (i < (totalorder+1u)*(totalorder+1u));

// Example i, i0, i1 values for totalorder = 5:
//                                    0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
//  static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
//  static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
	      
        unsigned int i0, i1;

        // Vertex DoFs
        if (i == 0)
          { i0 = 0; i1 = 0; }
        else if (i == 1)
          { i0 = 1; i1 = 0; }
        else if (i == 2)
          { i0 = 1; i1 = 1; }
        else if (i == 3)
          { i0 = 0; i1 = 1; }
        // Edge DoFs
        else if (i < totalorder + 3u)
          { i0 = i - 2; i1 = 0; }
        else if (i < 2u*totalorder + 2)
          { i0 = 1; i1 = i - totalorder - 1; }
        else if (i < 3u*totalorder + 1u)
          { i0 = i - 2u*totalorder; i1 = 1; }
        else if (i < 4u*totalorder)
          { i0 = 0; i1 = i - 3u*totalorder + 1; }
        // Interior DoFs
        else
          {
            unsigned int basisnum = i - 4*totalorder;
            i0 = square_number_column[basisnum] + 2;
            i1 = square_number_row[basisnum] + 2;
          }

        // Flip odd degree of freedom values if necessary
        // to keep continuity on sides
        Real f = 1.;

        if ((i0%2) && (i0 > 2) && (i1 == 0))
	  f = (elem->point(0) > elem->point(1))?-1.:1.;
        else if ((i0%2) && (i0>2) && (i1 == 1))
	  f = (elem->point(3) > elem->point(2))?-1.:1.;
        else if ((i0 == 0) && (i1%2) && (i1>2))
	  f = (elem->point(0) > elem->point(3))?-1.:1.;
        else if ((i0 == 1) && (i1%2) && (i1>2))
	  f = (elem->point(1) > elem->point(2))?-1.:1.;
	      
	switch (j)
	  {
	    // d()/dxi
	  case 0:		  		  
	    return f*(FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
		      FE<1,HIERARCHIC>::shape      (EDGE3, totalorder, i1,    eta));
	      
	    // d()/deta
	  case 1:		  		  
	    return f*(FE<1,HIERARCHIC>::shape      (EDGE3, totalorder, i0,    xi)*
		      FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i1, 0, eta));

	  default:
	    libmesh_error();
	  }

      }
      
    default:
      std::cerr << "ERROR: Unsupported element type!" << std::endl;
      libmesh_error();
    }

  return 0.;
}



template <>
Real FE<2,HIERARCHIC>::shape_second_deriv(const ElemType,
				          const Order,			    
				          const unsigned int,
				          const unsigned int,
				          const Point&)
{
  std::cerr << "Hierarchic polynomials require the element type\n"
	    << "because edge orientation is needed."
	    << std::endl;

  libmesh_error();
  return 0.;
}



template <>
Real FE<2,HIERARCHIC>::shape_second_deriv(const Elem* elem,
				          const Order order,
				          const unsigned int i,
				          const unsigned int j,
				          const Point& p)
{
  libmesh_assert (elem != NULL);

  // I have been lazy here and am using finite differences
  // to compute the derivatives!
  const Real eps = 1.e-6;
  Point pp, pm;
  unsigned int prevj = libMesh::invalid_uint;
	      
  switch (j)
  {
    //  d^2()/dxi^2
    case 0:
      {
        pp = Point(p(0)+eps, p(1));
        pm = Point(p(0)-eps, p(1));
        prevj = 0;
	break;
      }
  
    // d^2()/dxideta
    case 1:
      {
        pp = Point(p(0), p(1)+eps);
        pm = Point(p(0), p(1)-eps);
        prevj = 0;
	break;
      }

    // d^2()/deta^2
    case 2:
      {
        pp = Point(p(0), p(1)+eps);
        pm = Point(p(0), p(1)-eps);
        prevj = 1;
	break;
      }
    default:
      libmesh_error();
  }
  return (FE<2,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pp) -
	  FE<2,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pm)
	  )/2./eps;
}