quadrature_gauss_2d.c

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

// $Id: quadrature_gauss_2D.C 2949 2008-07-31 19:19:11Z jwpeterson $

// 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 "quadrature_gauss.h"
#include "quadrature_jacobi.h"


void QGauss::init_2D(const ElemType _type,
                     unsigned int p)
{
#if DIM > 1
  
  //-----------------------------------------------------------------------
  // 2D quadrature rules
  switch (_type)
    {


      //---------------------------------------------
      // Quadrilateral quadrature rules
    case QUAD4:
    case QUAD8:
    case QUAD9:
      {
	// We compute the 2D quadrature rule as a tensor
	// product of the 1D quadrature rule.
	//
	// For QUADs, a quadrature rule of order 'p' must be able to integrate
	// bilinear (p=1), biquadratic (p=2), bicubic (p=3), etc. polynomials of the form
	//
	// (x^p + x^{p-1} + ... + 1) * (y^p + y^{p-1} + ... + 1)
	//
	// These polynomials have terms *up to* degree 2p but they are *not* complete
	// polynomials of degree 2p. For example, when p=2 we have
	//        1
	//     x      y
	// x^2    xy     y^2
	//    yx^2   xy^2
	//       x^2y^2
	QGauss q1D(1,_order);
	q1D.init(EDGE2,p);
	tensor_product_quad( q1D );
	return;
      }

	    
      //---------------------------------------------
      // Triangle quadrature rules
    case TRI3:
    case TRI6:
      {
	switch(_order + 2*p)
	  {
	  case CONSTANT:
	  case FIRST:
	    {
	      // Exact for linears
	      _points.resize(1);
	      _weights.resize(1);
		  
	      _points[0](0) = 1.0L/3.0L;
	      _points[0](1) = 1.0L/3.0L;

	      _weights[0] = 0.5;

	      return;
	    }
	  case SECOND:
	    {
	      // Exact for quadratics
	      _points.resize(3);
	      _weights.resize(3);

	      // Alternate rule with points on ref. elt. boundaries.
	      // Not ideal for problems with material coefficient discontinuities
	      // aligned along element boundaries.
	      // _points[0](0) = .5;
	      // _points[0](1) = .5;
	      // _points[1](0) = 0.;
	      // _points[1](1) = .5;
	      // _points[2](0) = .5;
	      // _points[2](1) = .0;
	      
	      _points[0](0) = 2.0L/3.0L;
	      _points[0](1) = 1.0L/6.0L;

	      _points[1](0) = 1.0L/6.0L;
	      _points[1](1) = 2.0L/3.0L;

	      _points[2](0) = 1.0L/6.0L;
	      _points[2](1) = 1.0L/6.0L;


	      _weights[0] = 1.0L/6.0L;
	      _weights[1] = 1.0L/6.0L;
	      _weights[2] = 1.0L/6.0L;

	      return;
	    }
	  case THIRD:
	    {
	      // Exact for cubics
	      _points.resize(4);
	      _weights.resize(4);

	      // This rule is formed from a tensor product of appropriately-scaled
	      // Gauss and Jacobi rules.  (See also: the tensor_product_tri() and tensor_product_tet()
	      // routines in quadrature.C  For high orders these rules generally
	      // have too many points, but at extremely low order they are competitive
	      // and have the additional benefit of having all positive weights.
	      _points[0](0)=1.5505102572168219018027159252941e-01L;
	      _points[0](1)=1.7855872826361642311703513337422e-01L;
	      _points[1](0)=6.4494897427831780981972840747059e-01L;
	      _points[1](1)=7.5031110222608118177475598324603e-02L;
	      _points[2](0)=1.5505102572168219018027159252941e-01L;
	      _points[2](1)=6.6639024601470138670269327409637e-01L;
	      _points[3](0)=6.4494897427831780981972840747059e-01L;
	      _points[3](1)=2.8001991549907407200279599420481e-01L;

	      _weights[0]=1.5902069087198858469718450103758e-01L;
	      _weights[1]=9.0979309128011415302815498962418e-02L;
	      _weights[2]=1.5902069087198858469718450103758e-01L;
	      _weights[3]=9.0979309128011415302815498962418e-02L;

	      return;

	      
	      // The following third-order rule is quite commonly cited
	      // in the literature and most likely works fine.  However,
	      // we generally prefer a rule with all positive weights
	      // and an equal number of points, when available.
	      //
	      //  (allow_rules_with_negative_weights)
	      // {
	      //   // Exact for cubics
	      //   _points.resize(4);
	      //   _weights.resize(4);
	      //   
	      //   _points[0](0) = .33333333333333333333333333333333;
	      //   _points[0](1) = .33333333333333333333333333333333;
	      // 
	      //   _points[1](0) = .2;
	      //   _points[1](1) = .6;
	      // 
	      //   _points[2](0) = .2;
	      //   _points[2](1) = .2;
	      // 
	      //   _points[3](0) = .6;
	      //   _points[3](1) = .2;
	      // 
	      // 
	      //   _weights[0] = -27./96.;
	      //   _weights[1] =  25./96.;
	      //   _weights[2] =  25./96.;
	      //   _weights[3] =  25./96.;
	      // 
	      //   return;
	      // } // end if (allow_rules_with_negative_weights)
	      // Note: if !allow_rules_with_negative_weights, fall through to next case.
	    }


	    
	  case FOURTH:
	    {
	      // A degree 4 rule with six points.  This rule can be found in many places
	      // including:
	      //
	      // J.N. Lyness and D. Jespersen, Moderate degree symmetric
	      // quadrature rules for the triangle, J. Inst. Math. Appl.  15 (1975),
	      // 19--32.

	      _points.resize(6);
	      _weights.resize(6);
	      
	      // The points are arranged symmetrically in two sets ('a' and 'b') of three.
	      const Real a = 9.1576213509770743e-02L;  const Real wa = 5.4975871827660933e-02L;
	      const Real b = 4.4594849091596488e-01L;  const Real wb = 1.1169079483900573e-01L;

	      _points[0] = Point(a      ,           a); _weights[0] = wa;
	      _points[1] = Point(a      , 1.0L-2.0L*a); _weights[1] = wa;
	      _points[2] = Point(1.0L-2.0L*a,       a); _weights[2] = wa;

	      _points[3] = Point(b      ,          b); _weights[3] = wb;
	      _points[4] = Point(b      ,1.0L-2.0L*b); _weights[4] = wb;
	      _points[5] = Point(1.0L-2.0L*b,      b); _weights[5] = wb;
	      
	      return;
	    }


	    
	  case FIFTH:
	    {
	      // Exact for quintics
	      // Taken from "Quadrature on Simplices of Arbitrary
	      // Dimension" by Walkington
	      _points.resize(7);
	      _weights.resize(7);
		  
	      const Real b1 = 2.0L/7.0L + std::sqrt(15.0L)/21.0L;
	      const Real a1 = 1.0L - 2.0L*b1;
	      const Real b2 = 2.0L/7.0L - std::sqrt(15.0L)/21.0L;
	      const Real a2 = 1.0L - 2.0L*b2;
		  
	      _points[0](0) = 1.0L/3.0L;
	      _points[0](1) = 1.0L/3.0L;

	      _points[1](0) = a1;
	      _points[1](1) = b1;

	      _points[2](0) = b1;
	      _points[2](1) = a1;

	      _points[3](0) = b1;
	      _points[3](1) = b1;

	      _points[4](0) = a2;
	      _points[4](1) = b2;

	      _points[5](0) = b2;
	      _points[5](1) = a2;

	      _points[6](0) = b2;
	      _points[6](1) = b2;


	      _weights[0] = 9.0L/80.0L;
	      _weights[1] = 31.0L/480.0L + std::sqrt(15.0L)/2400.0L;
	      _weights[2] = _weights[1];
	      _weights[3] = _weights[1];
	      _weights[4] = 31.0L/480.0L - std::sqrt(15.0L)/2400.0L;
	      _weights[5] = _weights[4];
	      _weights[6] = _weights[4];

	      return;
	    }




	    
	    // A degree 7 rule with 12 points.  This rule can be found in:
	    //
	    // K. Gatermann, The construction of symmetric cubature
	    // formulas for the square and the triangle, Computing 40
	    // (1988), 229--240.
	    //
	    // This rule, which is provably minimal in the number of
	    // integration points, is said to be 'Ro3 invariant' which
	    // means that a given set of barycentric coordinates
	    // (z1,z2,z3) implies the quadrature points (z1,z2),
	    // (z3,z1), (z2,z3) which are formed by taking the first
	    // two entries in cyclic permutations of the barycentric
	    // point.  Barycentric coordinates are related in the
	    // sense that: z3 = 1 - z1 - z2.
	    //
	    // The 12-point sixth-order rule for triangles given in
	    // Flaherty's (http://www.cs.rpi.edu/~flaherje/FEM/fem6.ps)
	    // lecture notes has been removed in favor of this rule
	    // which is higher-order (for the same number of
	    // quadrature points) and has a few more digits of
	    // precision in the points and weights.  Some 10-point
	    // degree 6 rules exist for the triangle but they have
	    // quadrature points outside the region of integration.
	  case SIXTH:
	  case SEVENTH:
	    {
	      _points.resize (12);
	      _weights.resize(12);
	      
	      const unsigned int nrows=4;
	      
	      // In each of the rows below, the first two entries are (z1, z2) which imply
	      // z3.  The third entry is the weight for each of the points in the cyclic permutation.
	      const Real p[nrows][3] = {
		{6.2382265094402118e-02, 6.7517867073916085e-02, 2.6517028157436251e-02}, // group A
		{5.5225456656926611e-02, 3.2150249385198182e-01, 4.3881408714446055e-02}, // group B