quadrature_gauss_3d.c

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

	      // 		const Real b = 0.0909090909090909090909090909090909090909090909090909091;
		
	      // 		_points[5](0) = a;
	      // 		_points[5](1) = b;
	      // 		_points[5](2) = b;
		
	      // 		_points[6](0) = b;
	      // 		_points[6](1) = a;
	      // 		_points[6](2) = b;
		
	      // 		_points[7](0) = b;
	      // 		_points[7](1) = b;
	      // 		_points[7](2) = a;
		
	      // 		_points[8](0) = b;
	      // 		_points[8](1) = b;
	      // 		_points[8](2) = b;
	      // 	      }
	      // 	      {
	      // 		const Real a = 0.066550153573664;
	      // 		const Real b = 0.433449846426336;
		
	      // 		_points[9](0) = b;
	      // 		_points[9](1) = a;
	      // 		_points[9](2) = a;
		
	      // 		_points[10](0) = a;
	      // 		_points[10](1) = a;
	      // 		_points[10](2) = b;
		
	      // 		_points[11](0) = a;
	      // 		_points[11](1) = b;
	      // 		_points[11](2) = b;
		
	      // 		_points[12](0) = b;
	      // 		_points[12](1) = a;
	      // 		_points[12](2) = b;
		
	      // 		_points[13](0) = b;
	      // 		_points[13](1) = b;
	      // 		_points[13](2) = a;
		
	      // 		_points[14](0) = a;
	      // 		_points[14](1) = b;
	      // 		_points[14](2) = a;		
	      // 	      }
	      
	      // 	      _weights[0]  = 0.030283678097089;
	      // 	      _weights[1]  = 0.006026785714286;
	      // 	      _weights[2]  = _weights[1];
	      // 	      _weights[3]  = _weights[1];
	      // 	      _weights[4]  = _weights[1];
	      // 	      _weights[5]  = 0.011645249086029;
	      // 	      _weights[6]  = _weights[5];
	      // 	      _weights[7]  = _weights[5];
	      // 	      _weights[8]  = _weights[5];
	      // 	      _weights[9]  = 0.010949141561386;
	      // 	      _weights[10] = _weights[9];
	      // 	      _weights[11] = _weights[9];
	      // 	      _weights[12] = _weights[9];
	      // 	      _weights[13] = _weights[9];
	      // 	      _weights[14] = _weights[9];
	      
	      return;
	    }



	    
	    // This rule is originally from Keast:
	    //    Patrick Keast,
	    //    Moderate Degree Tetrahedral Quadrature Formulas,
	    //    Computer Methods in Applied Mechanics and Engineering,
	    //    Volume 55, Number 3, May 1986, pages 339-348.
	    //
	    // It is accurate on 6th-degree polynomials and has 24 points
	    // vs. 64 for the comparable conical product rule.
	    //
	    // Values copied 24th June 2008 from:
	    // http://people.scs.fsu.edu/~burkardt/f_src/keast/keast.f90
	  case SIXTH:
	    {
	      _points.resize (24);
	      _weights.resize(24);

	      // The raw data for the quadrature rule.
	      const Real p[4][4] = {
		{0.356191386222544953e+00 , 0.214602871259151684e+00 ,                       0., 0.00665379170969464506e+00}, // 4
		{0.877978124396165982e+00 , 0.0406739585346113397e+00,                       0., 0.00167953517588677620e+00}, // 4
		{0.0329863295731730594e+00, 0.322337890142275646e+00 ,                       0., 0.00922619692394239843e+00}, // 4
		{0.0636610018750175299e+00, 0.269672331458315867e+00 , 0.603005664791649076e+00, 0.00803571428571428248e+00}  // 12    
	      };

	      
	      // Now call the keast routine to generate _points and _weights
	      keast_rule(p, 4);

	      return;
	    }


	    // Keast's 31 point, 7th-order rule contains points on the reference
	    // element boundary, so we've decided not to include it here.
	    //
	    // Keast's 8th-order rule has 45 points.  and a negative
	    // weight, so if you've explicitly disallowed such rules
	    // you will fall through to the conical product rules
	    // below.
	  case SEVENTH:
	  case EIGHTH:
	    {
	      if (allow_rules_with_negative_weights)
		{
		  _points.resize (45);
		  _weights.resize(45);

		  // The raw data for the quadrature rule.
		  const Real p[7][4] = {
		    {0.250000000000000000e+00,                        0.,                        0.,   -0.393270066412926145e-01},  // 1
		    {0.617587190300082967e+00,  0.127470936566639015e+00,                        0.,    0.408131605934270525e-02},  // 4
		    {0.903763508822103123e+00,  0.320788303926322960e-01,                        0.,    0.658086773304341943e-03},  // 4
		    {0.497770956432810185e-01,                        0.,  0.450222904356718978e+00,    0.438425882512284693e-02},  // 6
		    {0.183730447398549945e+00,                        0.,  0.316269552601450060e+00,    0.138300638425098166e-01},  // 6
		    {0.231901089397150906e+00,  0.229177878448171174e-01,  0.513280033360881072e+00,    0.424043742468372453e-02},  // 12
		    {0.379700484718286102e-01,  0.730313427807538396e+00,  0.193746475248804382e+00,    0.223873973961420164e-02}   // 12
		  };

	      
		  // Now call the keast routine to generate _points and _weights
		  keast_rule(p, 7);

		  return;
		} // end if (allow_rules_with_negative_weights)
	      // Note: if !allow_rules_with_negative_weights, fall through to next case.
	    }

	    
	  case NINTH:
	  case TENTH:        
	  case ELEVENTH:
	  case TWELFTH:      
	  case THIRTEENTH:   
	  case FOURTEENTH:   
	  case FIFTEENTH:    
	  case SIXTEENTH:    
	  case SEVENTEENTH:  
	  case EIGHTTEENTH:  
	  case NINTEENTH:    
	  case TWENTIETH:    
	  case TWENTYFIRST:  
	  case TWENTYSECOND: 
	  case TWENTYTHIRD:
	    {
	      if (allow_rules_with_negative_weights)
		{
		  // The Grundmann-Moller rules are defined to arbitrary order and
		  // have significantly fewer evaluation points than concial product
		  // rules.  If you allow rules with negative weights, the GM rules
		  // will be more efficient in general, but may be more susceptible
		  // to round-off error.  Safest is to disallow rules with negative
		  // weights, but this decision should be made on a case-by-case basis.
		  QGrundmann_Moller gm_rule(3, _order);
		  gm_rule.init(_type, p);

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

		  return;
		}

	      else
		{
		  // The following quadrature rules are
		  // generated as conical products.  These
		  // tend to be non-optimal (use too many
		  // points, cluster points in certain
		  // regions of the domain) but they are
		  // quite easy to automatically generate
		  // using a 1D Gauss rule on [0,1] and two 
		  // 1D Jacobi-Gauss rules on [0,1].
	      
		  // Define the quadrature rules...
		  QGauss  gauss1D(1,static_cast<Order>(_order+2*p));
		  QJacobi jacA1D(1,static_cast<Order>(_order+2*p),1,0);
		  QJacobi jacB1D(1,static_cast<Order>(_order+2*p),2,0);

		  // The Gauss rule needs to be scaled to [0,1]
		  std::pair<Real, Real> old_range(-1,1);
		  std::pair<Real, Real> new_range(0,1);
		  gauss1D.scale(old_range,
				new_range);

		  // Compute the tensor product
		  tensor_product_tet(gauss1D, jacA1D, jacB1D);
		  return;
		}
	    }
	  default:
	    {
	      std::cout << "Quadrature rule not supported!" << std::endl;
		    
	      libmesh_error();
	    }
	  }
      }

      
	    
      //---------------------------------------------
      // Prism quadrature rules
    case PRISM6:
    case PRISM15:
    case PRISM18:
      {
	// We compute the 3D quadrature rule as a tensor
	// product of the 1D quadrature rule and a 2D
	// triangle quadrature rule
	    
	QGauss q1D(1,_order);
	QGauss q2D(2,_order);

	// Initialize 
	q1D.init(EDGE2,p);
	q2D.init(TRI3,p);

	tensor_product_prism(q1D, q2D);
	
	return;
      }
      

      
      //---------------------------------------------
      // Pyramid
    case PYRAMID5:
      {
	// We compute the Pyramid rule as a conical
	// product of the interval [0,1] and the
	// reference square [-1,1] x [-1,1] as per
	// Stroud, A.H. "Approximate Calculation of
	// Multiple Integrals."
	// This should be exact for quadratics, (Stroud, 32)

	// Get a rule for the reference quad
	QGauss q2D(2,_order); 
	q2D.init(QUAD4,p);
	
	// Get a rule for the interval [-1,1] 
	QGauss q1D(1,_order);
	q1D.init(EDGE2,p);

	// Storage for our temporary rule
	std::vector<Real> pts(q1D.n_points());
	std::vector<Real> wts(q1D.n_points());

	// Modify the 1D rule to fit in [0,1] instead
	for (unsigned int i=0; i<q1D.n_points(); i++)
	  {
	    pts[i] = 0.5*q1D.qp(i)(0) + 0.5;
	    wts[i] = 0.5*q1D.w(i); 
	  }

	// Allocate space for the new rule
	_points.resize(q1D.n_points() * q2D.n_points());
	_weights.resize(q1D.n_points() * q2D.n_points());
	      
	// Compute the conical product of the 1D and 2D rules
	unsigned int qp = 0;
	for (unsigned int i=0; i<q1D.n_points(); i++)
	  for (unsigned int j=0; j<q2D.n_points(); j++)
	    {
	      _points[qp](0) = q2D.qp(j)(0)*(1.0-pts[i]);
	      _points[qp](1) = q2D.qp(j)(1)*(1.0-pts[i]);
	      _points[qp](2) = pts[i];

	      _weights[qp] = 0.33333333333333 * // Scale factor!
		wts[i] * q2D.w(j);
		    
	      qp++;
	    }

	return;

      }


      
      //---------------------------------------------
      // Unsupported type
    default:
      {
	std::cerr << "ERROR: Unsupported type: " << _type << std::endl;
	libmesh_error();
      }
    }

  libmesh_error();

  return;
  
#endif
}