quadrature_monomial_2d.c

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

		  // {0.887188506449625, -0.236496718536120,  0.252528506429544},
		  // {0.494698820670197, -0.698953476086564,  0.361262323882172},
		  // {0.897495818279768, -0.900390774211580,  0.085464254086247}
		  //
		  // Second of 2 degree-8, 16 point rules given by Wissmann.
		  // Either of these will work, I just chose the one with points
		  // further into the element interior.
		  {0.0000000000000000e+00,  6.5956013196034176e-01,  4.5027677630559029e-01},
		  {0.0000000000000000e+00, -9.4914292304312538e-01,  1.6657042677781274e-01},
		  {9.5250946607156228e-01,  7.6505181955768362e-01,  9.8869459933431422e-02},
		  {5.3232745407420624e-01,  9.3697598108841598e-01,  1.5369674714081197e-01},
		  {6.8473629795173504e-01,  3.3365671773574759e-01,  3.9668697607290278e-01},
		  {2.3314324080140552e-01, -7.9583272377396852e-02,  3.5201436794569501e-01},
		  {9.2768331930611748e-01, -2.7224008061253425e-01,  1.8958905457779799e-01},
		  {4.5312068740374942e-01, -6.1373535339802760e-01,  3.7510100114758727e-01},
		  {8.3750364042281223e-01, -8.8847765053597136e-01,  1.2561879164007201e-01}
		};

	      _points.resize(16);
	      _weights.resize(16);

	      wissmann_rule(data, /*10*/ 9);

	      return; 
	    } // end case EIGHTH




	  case NINTH:
	    {
	      // A degree 9, 17-point rule due to Moller.  
	      //
	      // H.M. Moller,  Kubaturformeln mit minimaler Knotenzahl,
	      // Numer. Math.  25 (1976), 185--200.
	      //
	      // This rule is provably minimal in the number of points.
	      // A tensor-product rule accurate for "bi-ninth" degree polynomials would have 25 points.
	      const Real data[5][3] =
		{
		  {0.0000000000000000e+00, 0.0000000000000000e+00, 5.2674897119341563e-01}, // 1
		  {6.3068011973166885e-01, 9.6884996636197772e-01, 8.8879378170198706e-02}, // 4
		  {9.2796164595956966e-01, 7.5027709997890053e-01, 1.1209960212959648e-01}, // 4
		  {4.5333982113564719e-01, 5.2373582021442933e-01, 3.9828243926207009e-01}, // 4
		  {8.5261572933366230e-01, 7.6208328192617173e-02, 2.6905133763978080e-01}  // 4
		};

	      const unsigned int symmetry[5] = {
		0, // Single point
		4, // Rotational Invariant
		4, // Rotational Invariant
		4, // Rotational Invariant
		4  // Rotational Invariant
	      };
	      
	      _points.resize (17);
	      _weights.resize(17);

	      stroud_rule(data, symmetry, 5);

	      return;
	    } // end case NINTH




	  case TENTH:
	  case ELEVENTH:
	    {
	      // A degree 11, 24-point rule due to Cools and Haegemans.  
	      //
	      // R. Cools and A. Haegemans, Another step forward in searching for
	      // cubature formulae with a minimal number of knots for the square,
	      // Computing 40 (1988), 139--146.
	      // 
	      // P. Verlinden and R. Cools, The algebraic construction of a minimal
	      // cubature formula of degree 11 for the square, Cubature Formulas
	      // and their Applications (Russian) (Krasnoyarsk) (M.V. Noskov, ed.),
	      // 1994, pp. 13--23.
	      //
	      // This rule is provably minimal in the number of points.
	      // A tensor-product rule accurate for "bi-tenth" or "bi-eleventh" degree polynomials would have 36 points.
	      const Real data[6][3] =
		{
		  {6.9807610454956756e-01, 9.8263922354085547e-01, 4.8020763350723814e-02}, // 4
		  {9.3948638281673690e-01, 8.2577583590296393e-01, 6.6071329164550595e-02}, // 4
		  {9.5353952820153201e-01, 1.8858613871864195e-01, 9.7386777358668164e-02}, // 4
		  {3.1562343291525419e-01, 8.1252054830481310e-01, 2.1173634999894860e-01}, // 4
		  {7.1200191307533630e-01, 5.2532025036454776e-01, 2.2562606172886338e-01}, // 4
		  {4.2484724884866925e-01, 4.1658071912022368e-02, 3.5115871839824543e-01}  // 4
		};

	      const unsigned int symmetry[6] = {
		4, // Rotational Invariant
		4, // Rotational Invariant
		4, // Rotational Invariant
		4, // Rotational Invariant
		4, // Rotational Invariant
		4  // Rotational Invariant
	      };
	      
	      _points.resize (24);
	      _weights.resize(24);

	      stroud_rule(data, symmetry, 6);

	      return;
	    } // end case TENTH,ELEVENTH




	  case TWELFTH:
	  case THIRTEENTH:
	    {
	      // A degree 13, 33-point rule due to Cools and Haegemans.  
	      //
	      // R. Cools and A. Haegemans, Another step forward in searching for
	      // cubature formulae with a minimal number of knots for the square,
	      // Computing 40 (1988), 139--146.
	      // 
	      // A tensor-product rule accurate for "bi-12" or "bi-13" degree polynomials would have 49 points.
	      const Real data[9][3] =
		{
		  {0.0000000000000000e+00, 0.0000000000000000e+00, 3.0038211543122536e-01}, // 1
		  {9.8348668243987226e-01, 7.7880971155441942e-01, 2.9991838864499131e-02}, // 4
		  {8.5955600564163892e-01, 9.5729769978630736e-01, 3.8174421317083669e-02}, // 4
		  {9.5892517028753485e-01, 1.3818345986246535e-01, 6.0424923817749980e-02}, // 4
		  {3.9073621612946100e-01, 9.4132722587292523e-01, 7.7492738533105339e-02}, // 4
		  {8.5007667369974857e-01, 4.7580862521827590e-01, 1.1884466730059560e-01}, // 4
		  {6.4782163718701073e-01, 7.5580535657208143e-01, 1.2976355037000271e-01}, // 4
		  {7.0741508996444936e-02, 6.9625007849174941e-01, 2.1334158145718938e-01}, // 4
		  {4.0930456169403884e-01, 3.4271655604040678e-01, 2.5687074948196783e-01}  // 4
		};

	      const unsigned int symmetry[9] = {
		0, // Single point
		4, // Rotational Invariant
		4, // Rotational Invariant
		4, // Rotational Invariant
		4, // Rotational Invariant
		4, // Rotational Invariant
		4, // Rotational Invariant
		4, // Rotational Invariant
		4  // Rotational Invariant
	      };
	      
	      _points.resize (33);
	      _weights.resize(33);

	      stroud_rule(data, symmetry, 9);

	      return;
	    } // end case TWELFTH,THIRTEENTH




	  case FOURTEENTH:
	  case FIFTEENTH:
	  case SIXTEENTH:
	  case SEVENTEENTH:
	    {
	      // A degree 17, 60-point rule due to Cools and Haegemans.  
	      //
	      // R. Cools and A. Haegemans, Another step forward in searching for
	      // cubature formulae with a minimal number of knots for the square,
	      // Computing 40 (1988), 139--146.
	      // 
	      // A tensor-product rule accurate for "bi-14" or "bi-15" degree polynomials would have 64 points.
	      // A tensor-product rule accurate for "bi-16" or "bi-17" degree polynomials would have 81 points.
	      const Real data[10][3] =
		{
		  {9.8935307451260049e-01, 0.0000000000000000e+00, 2.0614915919990959e-02}, // 4
		  {3.7628520715797329e-01, 0.0000000000000000e+00, 1.2802571617990983e-01}, // 4
		  {9.7884827926223311e-01, 0.0000000000000000e+00, 5.5117395340318905e-03}, // 4
		  {8.8579472916411612e-01, 0.0000000000000000e+00, 3.9207712457141880e-02}, // 4
		  {1.7175612383834817e-01, 0.0000000000000000e+00, 7.6396945079863302e-02}, // 4
		  {5.9049927380600241e-01, 3.1950503663457394e-01, 1.4151372994997245e-01}, // 8
		  {7.9907913191686325e-01, 5.9797245192945738e-01, 8.3903279363797602e-02}, // 8
		  {8.0374396295874471e-01, 5.8344481776550529e-02, 6.0394163649684546e-02}, // 8
		  {9.3650627612749478e-01, 3.4738631616620267e-01, 5.7387752969212695e-02}, // 8
		  {9.8132117980545229e-01, 7.0600028779864611e-01, 2.1922559481863763e-02}, // 8
		};

	      const unsigned int symmetry[10] = {
		3, // Fully symmetric (x,0)
		3, // Fully symmetric (x,0)
		2, // Fully symmetric (x,x)
		2, // Fully symmetric (x,x)
		2, // Fully symmetric (x,x)
		1, // Fully symmetric (x,y) 
		1, // Fully symmetric (x,y)
		1, // Fully symmetric (x,y)
		1, // Fully symmetric (x,y)
		1  // Fully symmetric (x,y)
	      };
	      
	      _points.resize (60);
	      _weights.resize(60);

	      stroud_rule(data, symmetry, 10);

	      return;
	    } // end case FOURTEENTH through SEVENTEENTH
	    
	    
	    
	    // By default: construct and use a Gauss quadrature rule
	  default:
	    {
	      // Break out and fall down into the default: case for the
	      // outer switch statement.
	      break;
	    }
	    
	  } // end switch(_order + 2*p)
      } // end case QUAD4/8/9

      
      // By default: construct and use a Gauss quadrature rule
    default:
      {
	QGauss gauss_rule(2, _order);
	gauss_rule.init(_type, p);

	// Swap points and weights with the about-to-be destroyed rule.
	_points.swap (gauss_rule.get_points() );
	_weights.swap(gauss_rule.get_weights());

	return;
      }
    } // end switch (_type)
}