quadrature_monomial_2d.c

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

// $Id: quadrature_monomial_2D.C 2932 2008-07-14 22:29:43Z 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_monomial.h"
#include "quadrature_gauss.h"


void QMonomial::init_2D(const ElemType _type,
			unsigned int p)
{

  switch (_type)
    {
      //---------------------------------------------
      // Quadrilateral quadrature rules
    case QUAD4:
    case QUAD8:
    case QUAD9:
      {
	switch(_order + 2*p)
	  {
	  case SECOND:
	    {
	      // A degree=2 rule for the QUAD with 3 points.
	      // A tensor product degree-2 Gauss would have 4 points.
	      // This rule (or a variation on it) is probably available in
	      //
	      // A.H. Stroud, Approximate calculation of multiple integrals,
	      // Prentice-Hall, Englewood Cliffs, N.J., 1971.
	      //
	      // though I have never actually seen a reference for it.
	      // Luckily it's fairly easy to derive, which is what I've done
	      // here [JWP].
	      const Real
		s=std::sqrt(1./3.),
		t=std::sqrt(2./3.);
	      
	      const Real data[2][3] =
		{
		  {0.0,  s,  2.0},
		  {  t, -s,  1.0}
		};

	      _points.resize(3);
	      _weights.resize(3);

	      wissmann_rule(data, 2);
	      
	      return;
	    } // end case SECOND


	    
	  case FOURTH:
	    {
	      // A pair of degree=4 rules for the QUAD "C2" due to
	      // Wissmann and Becker. These rules both have six points.
	      // A tensor product degree-4 Gauss would have 9 points.
	      //
	      // J. W. Wissmann and T. Becker, Partially symmetric cubature 
	      // formulas for even degrees of exactness, SIAM J. Numer. Anal.  23
	      // (1986), 676--685.                                           
	      const Real data[4][3] =
		{
		  // First of 2 degree-4 rules given by Wissmann
		  {0.0000000000000000e+00,  0.0000000000000000e+00,  1.1428571428571428e+00},
		  {0.0000000000000000e+00,  9.6609178307929590e-01,  4.3956043956043956e-01},
		  {8.5191465330460049e-01,  4.5560372783619284e-01,  5.6607220700753210e-01},
		  {6.3091278897675402e-01, -7.3162995157313452e-01,  6.4271900178367668e-01}
		  //
		  // Second of 2 degree-4 rules given by Wissmann.  These both
		  // yield 4th-order accurate rules, I just chose the one that
		  // happened to contain the origin.
		  // {0.000000000000000, -0.356822089773090,  1.286412084888852}, 
		  // {0.000000000000000,  0.934172358962716,  0.491365692888926},
		  // {0.774596669241483,  0.390885162530071,  0.761883709085613},
		  // {0.774596669241483, -0.852765377881771,  0.349227402025498} 
		};

	      _points.resize(6);
	      _weights.resize(6);

	      wissmann_rule(data, 4);
	      
	      return;
	    } // end case FOURTH



	    
	  case FIFTH:
	    {
	      // A degree 5, 7-point rule due to Stroud.  
	      //
	      // A.H. Stroud, Approximate calculation of multiple integrals,
	      // Prentice-Hall, Englewood Cliffs, N.J., 1971.
	      //
	      // This rule is provably minimal in the number of points.
	      // A tensor-product rule accurate for "bi-quintic" polynomials would have 9 points.
	      const Real data[3][3] =
		{
		  {0.0000000000000000e+00, 0.0000000000000000e+00, 1.1428571428571428e+00}, // 1
		  {0.0000000000000000e+00, 9.6609178307929590e-01, 3.1746031746031746e-01}, // 2
		  {7.7459666924148337e-01, 5.7735026918962576e-01, 5.5555555555555555e-01}  // 4
		};

	      const unsigned int symmetry[3] = {
		0, // Origin
		7, // Central Symmetry
		6  // Rectangular
	      };
	      
	      _points.resize (7);
	      _weights.resize(7);

	      stroud_rule(data, symmetry, 3);
	      
	      return;
	    } // end case FIFTH



	    
	  case SIXTH:
	    {
	      // A pair of degree=6 rules for the QUAD "C2" due to
	      // Wissmann and Becker. These rules both have 10 points.
	      // A tensor product degree-6 Gauss would have 16 points.
	      //
	      // J. W. Wissmann and T. Becker, Partially symmetric cubature 
	      // formulas for even degrees of exactness, SIAM J. Numer. Anal.  23
	      // (1986), 676--685.                                           
	      const Real data[6][3] =
		{
		  // First of 2 degree-6, 10 point rules given by Wissmann
		  // {0.000000000000000,  0.836405633697626,  0.455343245714174},
		  // {0.000000000000000, -0.357460165391307,  0.827395973202966},
		  // {0.888764014654765,  0.872101531193131,  0.144000884599645},
		  // {0.604857639464685,  0.305985162155427,  0.668259104262665},
		  // {0.955447506641064, -0.410270899466658,  0.225474004890679},
		  // {0.565459993438754, -0.872869311156879,  0.320896396788441} 
		  //
		  // Second of 2 degree-6, 10 point rules given by Wissmann.
		  // Either of these will work, I just chose the one with points
		  // slightly further into the element interior.
		  {0.0000000000000000e+00,  8.6983337525005900e-01,  3.9275059096434794e-01},
		  {0.0000000000000000e+00, -4.7940635161211124e-01,  7.5476288124261053e-01},
		  {8.6374282634615388e-01,  8.0283751620765670e-01,  2.0616605058827902e-01},
		  {5.1869052139258234e-01,  2.6214366550805818e-01,  6.8999213848986375e-01},
		  {9.3397254497284950e-01, -3.6309658314806653e-01,  2.6051748873231697e-01},
		  {6.0897753601635630e-01, -8.9660863276245265e-01,  2.6956758608606100e-01}
		};

	      _points.resize(10);
	      _weights.resize(10);

	      wissmann_rule(data, 6);

	      return; 
	    } // end case SIXTH



	    
	  case SEVENTH:
	    {
	      // A degree 7, 12-point rule due to Stroud.  
	      //
	      // A.H. Stroud, Approximate calculation of multiple integrals,
	      // Prentice-Hall, Englewood Cliffs, N.J., 1971.
	      //
	      // This rule is fully-symmetric and provably minimal in the number of points.
	      // A tensor-product rule accurate for "bi-septic" polynomials would have 16 points.
	      const Real data[3][3] =
		{
		  {9.2582009977255146e-01, 0.0000000000000000e+00, 2.4197530864197530e-01}, // 4
		  {3.8055443320831565e-01, 0.0000000000000000e+00, 5.2059291666739445e-01}, // 4
		  {8.0597978291859874e-01, 0.0000000000000000e+00, 2.3743177469063023e-01}  // 4
		};

	      const unsigned int symmetry[3] = {
		3, // Full Symmetry, (x,0)
		2, // Full Symmetry, (x,x)
		2  // Full Symmetry, (x,x)
	      };
	      
	      _points.resize (12);
	      _weights.resize(12);

	      stroud_rule(data, symmetry, 3);

	      return;
	    } // end case SEVENTH



	    
	  case EIGHTH:
	    {
	      // A pair of degree=8 rules for the QUAD "C2" due to
	      // Wissmann and Becker. These rules both have 16 points.
	      // A tensor product degree-6 Gauss would have 25 points.
	      //
	      // J. W. Wissmann and T. Becker, Partially symmetric cubature 
	      // formulas for even degrees of exactness, SIAM J. Numer. Anal.  23
	      // (1986), 676--685.                                           
	      const Real data[10][3] =
		{
		  // First of 2 degree-8, 16 point rules given by Wissmann
		  // {0.000000000000000,  0.000000000000000,  0.055364705621440},
		  // {0.000000000000000,  0.757629177660505,  0.404389368726076},
		  // {0.000000000000000, -0.236871842255702,  0.533546604952635},
		  // {0.000000000000000, -0.989717929044527,  0.117054188786739},
		  // {0.639091304900370,  0.950520955645667,  0.125614417613747},
		  // {0.937069076924990,  0.663882736885633,  0.136544584733588},
		  // {0.537083530541494,  0.304210681724104,  0.483408479211257},