fe_clough_shape_2d.c

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

// $Id: fe_clough_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++ inlcludes

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


// Anonymous namespace for persistant variables.
// This allows us to cache the global-to-local mapping transformation
// FIXME: This should also screw up multithreading royally
namespace
{
  static unsigned int old_elem_id = libMesh::invalid_uint;
  // Coefficient naming: d(1)d(2n) is the coefficient of the
  // global shape function corresponding to value 1 in terms of the
  // local shape function corresponding to normal derivative 2
  static Real d1d2n, d1d3n, d2d3n, d2d1n, d3d1n, d3d2n;
  static Real d1xd1x, d1xd1y, d1xd2n, d1xd3n;
  static Real d1yd1x, d1yd1y, d1yd2n, d1yd3n;
  static Real d2xd2x, d2xd2y, d2xd3n, d2xd1n;
  static Real d2yd2x, d2yd2y, d2yd3n, d2yd1n;
  static Real d3xd3x, d3xd3y, d3xd1n, d3xd2n;
  static Real d3yd3x, d3yd3y, d3yd1n, d3yd2n;
  static Real d1nd1n, d2nd2n, d3nd3n;
  // Normal vector naming: N01x is the x component of the
  // unit vector at point 0 normal to (possibly curved) side 01
  static Real N01x, N01y, N10x, N10y;
  static Real N02x, N02y, N20x, N20y;
  static Real N21x, N21y, N12x, N12y;

Real clough_raw_shape_second_deriv(const unsigned int basis_num,
                                   const unsigned int deriv_type,
                                   const Point& p);
Real clough_raw_shape_deriv(const unsigned int basis_num,
                            const unsigned int deriv_type,
                            const Point& p);
Real clough_raw_shape(const unsigned int basis_num,
                      const Point& p);
unsigned char subtriangle_lookup(const Point& p);


// Compute the static coefficients for an element
void clough_compute_coefs(const Elem* elem)
{
  // Coefficients are cached from old elements
  if (elem->id() == old_elem_id)
    return;

  old_elem_id = elem->id();

  const Order mapping_order        (elem->default_order());
  const ElemType mapping_elem_type (elem->type());
  const int n_mapping_shape_functions =
    FE<2,LAGRANGE>::n_shape_functions(mapping_elem_type,
				      mapping_order);

  // Degrees of freedom are at vertices and edge midpoints
  std::vector<Point> dofpt;
  dofpt.push_back(Point(0,0));
  dofpt.push_back(Point(1,0));
  dofpt.push_back(Point(0,1));
  dofpt.push_back(Point(1/2.,1/2.));
  dofpt.push_back(Point(0,1/2.));
  dofpt.push_back(Point(1/2.,0));

  // Mapping functions - first derivatives at each dofpt
  std::vector<Real> dxdxi(6), dxdeta(6), dydxi(6), dydeta(6);
  std::vector<Real> dxidx(6), detadx(6), dxidy(6), detady(6);

  for (int p = 0; p != 6; ++p)
    {
//      std::cerr << p << ' ' << dofpt[p];
      for (int i = 0; i != n_mapping_shape_functions; ++i)
        {
          const Real ddxi = FE<2,LAGRANGE>::shape_deriv 
            (mapping_elem_type, mapping_order, i, 0, dofpt[p]);
          const Real ddeta = FE<2,LAGRANGE>::shape_deriv 
            (mapping_elem_type, mapping_order, i, 1, dofpt[p]);

//      std::cerr << ddxi << ' ';
//      std::cerr << ddeta << std::endl;

          dxdxi[p] += elem->point(i)(0) * ddxi;
          dydxi[p] += elem->point(i)(1) * ddxi;
          dxdeta[p] += elem->point(i)(0) * ddeta;
          dydeta[p] += elem->point(i)(1) * ddeta;
        }

//      for (int i = 0; i != 12; ++i)
//          std::cerr << i << ' ' << clough_raw_shape(i, dofpt[p]) << std::endl;

//      std::cerr << elem->point(p)(0) << ' ';
//      std::cerr << elem->point(p)(1) << ' ';
//      std::cerr << dxdxi[p] << ' ';
//      std::cerr << dydxi[p] << ' ';
//      std::cerr << dxdeta[p] << ' ';
//      std::cerr << dydeta[p] << std::endl << std::endl;

      const Real inv_jac = 1. / (dxdxi[p]*dydeta[p] - 
        dxdeta[p]*dydxi[p]);
      dxidx[p] = dydeta[p] * inv_jac;
      dxidy[p] = - dxdeta[p] * inv_jac;
      detadx[p] = - dydxi[p] * inv_jac;
      detady[p] = dxdxi[p] * inv_jac;
    }

  // Calculate midpoint normal vectors
  Real N1x = dydeta[3] - dydxi[3];
  Real N1y = dxdxi[3] - dxdeta[3];
  Real Nlength = std::sqrt(static_cast<Real>(N1x*N1x + N1y*N1y));
  N1x /= Nlength; N1y /= Nlength;

  Real N2x = - dydeta[4];
  Real N2y = dxdeta[4];
  Nlength = std::sqrt(static_cast<Real>(N2x*N2x + N2y*N2y));
  N2x /= Nlength; N2y /= Nlength;

  Real N3x = dydxi[5];
  Real N3y = - dxdxi[5];
  Nlength = std::sqrt(static_cast<Real>(N3x*N3x + N3y*N3y));
  N3x /= Nlength; N3y /= Nlength;

  // Calculate corner normal vectors (used for reduced element)
  N01x = dydxi[0];
  N01y = - dxdxi[0];
  Nlength = std::sqrt(static_cast<Real>(N01x*N01x + N01y*N01y));
  N01x /= Nlength; N01y /= Nlength;

  N10x = dydxi[1];
  N10y = - dxdxi[1];
  Nlength = std::sqrt(static_cast<Real>(N10x*N10x + N10y*N10y));
  N10x /= Nlength; N10y /= Nlength;

  N02x = - dydeta[0];
  N02y = dxdeta[0];
  Nlength = std::sqrt(static_cast<Real>(N02x*N02x + N02y*N02y));
  N02x /= Nlength; N02y /= Nlength;

  N20x = - dydeta[2];
  N20y = dxdeta[2];
  Nlength = std::sqrt(static_cast<Real>(N20x*N20x + N20y*N20y));
  N20x /= Nlength; N20y /= Nlength;

  N12x = dydeta[1] - dydxi[1];
  N12y = dxdxi[1] - dxdeta[1];
  Nlength = std::sqrt(static_cast<Real>(N12x*N12x + N12y*N12y));
  N12x /= Nlength; N12y /= Nlength;

  N21x = dydeta[1] - dydxi[1];
  N21y = dxdxi[1] - dxdeta[1];
  Nlength = std::sqrt(static_cast<Real>(N21x*N21x + N21y*N21y));
  N21x /= Nlength; N21y /= Nlength;

//  for (int i=0; i != 6; ++i) {
//    std::cerr << elem->node(i) << ' ';
//  }
//  std::cerr << std::endl;

//  for (int i=0; i != 6; ++i) {
//    std::cerr << elem->point(i)(0) << ' ';
//    std::cerr << elem->point(i)(1) << ' ';
//  }
//  std::cerr << std::endl;


  // give normal vectors a globally unique orientation

  if (elem->point(2) < elem->point(1))
    {
//      std::cerr << "Flipping nodes " << elem->node(2);
//      std::cerr << " and " << elem->node(1);
//      std::cerr << " around node " << elem->node(4);
//      std::cerr << std::endl;
      N1x = -N1x; N1y = -N1y;
      N12x = -N12x; N12y = -N12y;
      N21x = -N21x; N21y = -N21y;
    }
  else
    {
//      std::cerr << "Not flipping nodes " << elem->node(2);
//      std::cerr << " and " << elem->node(1);
//      std::cerr << " around node " << elem->node(4);
//      std::cerr << std::endl;
    }
  if (elem->point(0) < elem->point(2))
    {
//      std::cerr << "Flipping nodes " << elem->node(0);
//      std::cerr << " and " << elem->node(2);
//      std::cerr << " around node " << elem->node(5);
//      std::cerr << std::endl;
//      std::cerr << N2x << ' ' << N2y << std::endl;
      N2x = -N2x; N2y = -N2y;
      N02x = -N02x; N02y = -N02y;
      N20x = -N20x; N20y = -N20y;
//      std::cerr << N2x << ' ' << N2y << std::endl;
    }
  else
    {
//      std::cerr << "Not flipping nodes " << elem->node(0);
//      std::cerr << " and " << elem->node(2);
//      std::cerr << " around node " << elem->node(5);
//      std::cerr << std::endl;
    }
  if (elem->point(1) < elem->point(0))
    {
//      std::cerr << "Flipping nodes " << elem->node(1);
//      std::cerr << " and " << elem->node(0);
//      std::cerr << " around node " << elem->node(3);
//      std::cerr << std::endl;
      N3x = -N3x;
      N3y = -N3y;
      N01x = -N01x; N01y = -N01y;
      N10x = -N10x; N10y = -N10y;
    }
  else
    {
//      std::cerr << "Not flipping nodes " << elem->node(1);
//      std::cerr << " and " << elem->node(0);
//      std::cerr << " around node " << elem->node(3);
//      std::cerr << std::endl;
    }

//  std::cerr << N2x << ' ' << N2y << std::endl;

  // Cache basis function gradients
  // FIXME: the raw_shape calls shouldn't be done on every element!
  // FIXME: I should probably be looping, too...
  // Gradient naming: d(1)d(2n)d(xi) is the xi component of the
  // gradient of the 
  // local basis function corresponding to value 1 at the node
  // corresponding to normal vector 2

  Real d1d2ndxi   = clough_raw_shape_deriv(0, 0, dofpt[4]);
  Real d1d2ndeta  = clough_raw_shape_deriv(0, 1, dofpt[4]);
  Real d1d2ndx = d1d2ndxi * dxidx[4] + d1d2ndeta * detadx[4];
  Real d1d2ndy = d1d2ndxi * dxidy[4] + d1d2ndeta * detady[4];
  Real d1d3ndxi   = clough_raw_shape_deriv(0, 0, dofpt[5]);
  Real d1d3ndeta  = clough_raw_shape_deriv(0, 1, dofpt[5]);
  Real d1d3ndx = d1d3ndxi * dxidx[5] + d1d3ndeta * detadx[5];
  Real d1d3ndy = d1d3ndxi * dxidy[5] + d1d3ndeta * detady[5];
  Real d2d3ndxi   = clough_raw_shape_deriv(1, 0, dofpt[5]);
  Real d2d3ndeta  = clough_raw_shape_deriv(1, 1, dofpt[5]);
  Real d2d3ndx = d2d3ndxi * dxidx[5] + d2d3ndeta * detadx[5];
  Real d2d3ndy = d2d3ndxi * dxidy[5] + d2d3ndeta * detady[5];
  Real d2d1ndxi   = clough_raw_shape_deriv(1, 0, dofpt[3]);
  Real d2d1ndeta  = clough_raw_shape_deriv(1, 1, dofpt[3]);
  Real d2d1ndx = d2d1ndxi * dxidx[3] + d2d1ndeta * detadx[3];
  Real d2d1ndy = d2d1ndxi * dxidy[3] + d2d1ndeta * detady[3];
  Real d3d1ndxi   = clough_raw_shape_deriv(2, 0, dofpt[3]);
  Real d3d1ndeta  = clough_raw_shape_deriv(2, 1, dofpt[3]);
  Real d3d1ndx = d3d1ndxi * dxidx[3] + d3d1ndeta * detadx[3];
  Real d3d1ndy = d3d1ndxi * dxidy[3] + d3d1ndeta * detady[3];
  Real d3d2ndxi   = clough_raw_shape_deriv(2, 0, dofpt[4]);
  Real d3d2ndeta  = clough_raw_shape_deriv(2, 1, dofpt[4]);
  Real d3d2ndx = d3d2ndxi * dxidx[4] + d3d2ndeta * detadx[4];
  Real d3d2ndy = d3d2ndxi * dxidy[4] + d3d2ndeta * detady[4];
  Real d1xd2ndxi  = clough_raw_shape_deriv(3, 0, dofpt[4]);
  Real d1xd2ndeta = clough_raw_shape_deriv(3, 1, dofpt[4]);
  Real d1xd2ndx = d1xd2ndxi * dxidx[4] + d1xd2ndeta * detadx[4];
  Real d1xd2ndy = d1xd2ndxi * dxidy[4] + d1xd2ndeta * detady[4];
  Real d1xd3ndxi  = clough_raw_shape_deriv(3, 0, dofpt[5]);
  Real d1xd3ndeta = clough_raw_shape_deriv(3, 1, dofpt[5]);
  Real d1xd3ndx = d1xd3ndxi * dxidx[5] + d1xd3ndeta * detadx[5];
  Real d1xd3ndy = d1xd3ndxi * dxidy[5] + d1xd3ndeta * detady[5];
  Real d1yd2ndxi  = clough_raw_shape_deriv(4, 0, dofpt[4]);
  Real d1yd2ndeta = clough_raw_shape_deriv(4, 1, dofpt[4]);
  Real d1yd2ndx = d1yd2ndxi * dxidx[4] + d1yd2ndeta * detadx[4];
  Real d1yd2ndy = d1yd2ndxi * dxidy[4] + d1yd2ndeta * detady[4];
  Real d1yd3ndxi  = clough_raw_shape_deriv(4, 0, dofpt[5]);
  Real d1yd3ndeta = clough_raw_shape_deriv(4, 1, dofpt[5]);
  Real d1yd3ndx = d1yd3ndxi * dxidx[5] + d1yd3ndeta * detadx[5];
  Real d1yd3ndy = d1yd3ndxi * dxidy[5] + d1yd3ndeta * detady[5];
  Real d2xd3ndxi  = clough_raw_shape_deriv(5, 0, dofpt[5]);
  Real d2xd3ndeta = clough_raw_shape_deriv(5, 1, dofpt[5]);
  Real d2xd3ndx = d2xd3ndxi * dxidx[5] + d2xd3ndeta * detadx[5];
  Real d2xd3ndy = d2xd3ndxi * dxidy[5] + d2xd3ndeta * detady[5];
  Real d2xd1ndxi  = clough_raw_shape_deriv(5, 0, dofpt[3]);
  Real d2xd1ndeta = clough_raw_shape_deriv(5, 1, dofpt[3]);
  Real d2xd1ndx = d2xd1ndxi * dxidx[3] + d2xd1ndeta * detadx[3];
  Real d2xd1ndy = d2xd1ndxi * dxidy[3] + d2xd1ndeta * detady[3];
  Real d2yd3ndxi  = clough_raw_shape_deriv(6, 0, dofpt[5]);
  Real d2yd3ndeta = clough_raw_shape_deriv(6, 1, dofpt[5]);
  Real d2yd3ndx = d2yd3ndxi * dxidx[5] + d2yd3ndeta * detadx[5];
  Real d2yd3ndy = d2yd3ndxi * dxidy[5] + d2yd3ndeta * detady[5];
  Real d2yd1ndxi  = clough_raw_shape_deriv(6, 0, dofpt[3]);
  Real d2yd1ndeta = clough_raw_shape_deriv(6, 1, dofpt[3]);
  Real d2yd1ndx = d2yd1ndxi * dxidx[3] + d2yd1ndeta * detadx[3];
  Real d2yd1ndy = d2yd1ndxi * dxidy[3] + d2yd1ndeta * detady[3];
  Real d3xd1ndxi  = clough_raw_shape_deriv(7, 0, dofpt[3]);
  Real d3xd1ndeta = clough_raw_shape_deriv(7, 1, dofpt[3]);
  Real d3xd1ndx = d3xd1ndxi * dxidx[3] + d3xd1ndeta * detadx[3];
  Real d3xd1ndy = d3xd1ndxi * dxidy[3] + d3xd1ndeta * detady[3];
  Real d3xd2ndxi  = clough_raw_shape_deriv(7, 0, dofpt[4]);
  Real d3xd2ndeta = clough_raw_shape_deriv(7, 1, dofpt[4]);
  Real d3xd2ndx = d3xd2ndxi * dxidx[4] + d3xd2ndeta * detadx[4];
  Real d3xd2ndy = d3xd2ndxi * dxidy[4] + d3xd2ndeta * detady[4];
  Real d3yd1ndxi  = clough_raw_shape_deriv(8, 0, dofpt[3]);
  Real d3yd1ndeta = clough_raw_shape_deriv(8, 1, dofpt[3]);
  Real d3yd1ndx = d3yd1ndxi * dxidx[3] + d3yd1ndeta * detadx[3];
  Real d3yd1ndy = d3yd1ndxi * dxidy[3] + d3yd1ndeta * detady[3];
  Real d3yd2ndxi  = clough_raw_shape_deriv(8, 0, dofpt[4]);
  Real d3yd2ndeta = clough_raw_shape_deriv(8, 1, dofpt[4]);
  Real d3yd2ndx = d3yd2ndxi * dxidx[4] + d3yd2ndeta * detadx[4];
  Real d3yd2ndy = d3yd2ndxi * dxidy[4] + d3yd2ndeta * detady[4];
  Real d1nd1ndxi  = clough_raw_shape_deriv(9, 0, dofpt[3]);
  Real d1nd1ndeta = clough_raw_shape_deriv(9, 1, dofpt[3]);
  Real d1nd1ndx = d1nd1ndxi * dxidx[3] + d1nd1ndeta * detadx[3];
  Real d1nd1ndy = d1nd1ndxi * dxidy[3] + d1nd1ndeta * detady[3];
  Real d2nd2ndxi  = clough_raw_shape_deriv(10, 0, dofpt[4]);
  Real d2nd2ndeta = clough_raw_shape_deriv(10, 1, dofpt[4]);
  Real d2nd2ndx = d2nd2ndxi * dxidx[4] + d2nd2ndeta * detadx[4];
  Real d2nd2ndy = d2nd2ndxi * dxidy[4] + d2nd2ndeta * detady[4];
  Real d3nd3ndxi  = clough_raw_shape_deriv(11, 0, dofpt[5]);
  Real d3nd3ndeta = clough_raw_shape_deriv(11, 1, dofpt[5]);
  Real d3nd3ndx = d3nd3ndxi * dxidx[3] + d3nd3ndeta * detadx[3];
  Real d3nd3ndy = d3nd3ndxi * dxidy[3] + d3nd3ndeta * detady[3];

  Real d1xd1dxi   = clough_raw_shape_deriv(3, 0, dofpt[0]);
  Real d1xd1deta  = clough_raw_shape_deriv(3, 1, dofpt[0]);
  Real d1xd1dx    = d1xd1dxi * dxidx[0] + d1xd1deta * detadx[0];
  Real d1xd1dy    = d1xd1dxi * dxidy[0] + d1xd1deta * detady[0];
  Real d1yd1dxi   = clough_raw_shape_deriv(4, 0, dofpt[0]);
  Real d1yd1deta  = clough_raw_shape_deriv(4, 1, dofpt[0]);
  Real d1yd1dx    = d1yd1dxi * dxidx[0] + d1yd1deta * detadx[0];
  Real d1yd1dy    = d1yd1dxi * dxidy[0] + d1yd1deta * detady[0];
  Real d2xd2dxi   = clough_raw_shape_deriv(5, 0, dofpt[1]);
  Real d2xd2deta  = clough_raw_shape_deriv(5, 1, dofpt[1]);
  Real d2xd2dx    = d2xd2dxi * dxidx[1] + d2xd2deta * detadx[1];
  Real d2xd2dy    = d2xd2dxi * dxidy[1] + d2xd2deta * detady[1];

//  std::cerr << dofpt[4](0) << ' ';
//  std::cerr << dofpt[4](1) << ' ';
//  std::cerr << (int)subtriangle_lookup(dofpt[5]) << ' ';
//  std::cerr << dxdxi[4] << ' ';
//  std::cerr << dxdeta[4] << ' ';
//  std::cerr << dydxi[4] << ' ';
//  std::cerr << dydeta[4] << ' ';
//  std::cerr << dxidx[4] << ' ';
//  std::cerr << dxidy[4] << ' ';
//  std::cerr << detadx[4] << ' ';
//  std::cerr << detady[4] << ' ';
//  std::cerr << N1x << ' ';
//  std::cerr << N1y << ' ';
//  std::cerr << d2yd1ndxi << ' ';
//  std::cerr << d2yd1ndeta << ' ';
//  std::cerr << d2yd1ndx << ' ';
//  std::cerr << d2yd1ndy << std::endl;

  Real d2yd2dxi   = clough_raw_shape_deriv(6, 0, dofpt[1]);
  Real d2yd2deta  = clough_raw_shape_deriv(6, 1, dofpt[1]);
  Real d2yd2dx    = d2yd2dxi * dxidx[1] + d2yd2deta * detadx[1];
  Real d2yd2dy    = d2yd2dxi * dxidy[1] + d2yd2deta * detady[1];
  Real d3xd3dxi   = clough_raw_shape_deriv(7, 0, dofpt[2]);
  Real d3xd3deta  = clough_raw_shape_deriv(7, 1, dofpt[2]);
  Real d3xd3dx    = d3xd3dxi * dxidx[1] + d3xd3deta * detadx[1];
  Real d3xd3dy    = d3xd3dxi * dxidy[1] + d3xd3deta * detady[1];
  Real d3yd3dxi   = clough_raw_shape_deriv(8, 0, dofpt[2]);
  Real d3yd3deta  = clough_raw_shape_deriv(8, 1, dofpt[2]);
  Real d3yd3dx    = d3yd3dxi * dxidx[1] + d3yd3deta * detadx[1];
  Real d3yd3dy    = d3yd3dxi * dxidy[1] + d3yd3deta * detady[1];

  // Calculate normal derivatives

  Real d1nd1ndn = d1nd1ndx * N1x + d1nd1ndy * N1y;