quadratureformular.cpp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

/*
Region: Simplex 
Dimension: 3 
Degree: 3 
Points: 5 
Structure: Fully symmetric 
Rule struct: 0 0 0 0 0 0 1 1 0 0 0 
Generator: [ Fully symmetric ] 
( 0.25, 0.25, 0.25, ) 
Corresponding weight: 
-0.133333333333333333333333333333333,
Generator: [ Fully symmetric ] 
( 0.166666666666666666666666666666666, 
0.166666666666666666666666666666666, 
0.166666666666666666666666666666666, 
) 
Corresponding weight: 
0.075,

Region: Simplex 
Dimension: 3 
Degree: 4 
Points: 11 
Structure: Fully symmetric 
Rule struct: 0 0 0 0 0 0 1 1 1 0 0 
Generator: [ Fully symmetric ] 
( 0.25, 0.25, 0.25, ) 
Corresponding weight: 
-0.0131555555555555555555555555555555,
Generator: [ Fully symmetric ] 
( 0.0714285714285714285714285714285714, 
0.0714285714285714285714285714285714, 
0.0714285714285714285714285714285714, 
) 
Corresponding weight: 
10 ^ -3 x 7.62222222222222222222222222222222,

Generator: [ Fully symmetric ] 
( 0.399403576166799204996102147461640, 
0.399403576166799204996102147461640, 
0.100596423833200795003897852538359, 
) 
Corresponding weight: 
0.0248888888888888888888888888888888,


Region: Simplex 
Dimension: 3 
Degree: 6 
Points: 24 
Structure: Fully symmetric 
Rule struct: 0 0 0 0 0 0 0 3 0 1 0 
Generator: [ Fully symmetric ] 
( 0.214602871259152029288839219386284, 
0.214602871259152029288839219386284, 
0.214602871259152029288839219386284, 
) 
Corresponding weight: 
10 ^ -3 x 6.65379170969458201661510459291332,
Generator: [ Fully symmetric ] 
( 0.0406739585346113531155794489564100, 
0.0406739585346113531155794489564100, 
0.0406739585346113531155794489564100, 
) 
Corresponding weight: 
10 ^ -3 x 1.67953517588677382466887290765614,

Generator: [ Fully symmetric ] 
( 0.322337890142275510343994470762492, 
0.322337890142275510343994470762492, 
0.322337890142275510343994470762492, 
) 
Corresponding weight: 
10 ^ -3 x 9.22619692394245368252554630895433,

Generator: [ Fully symmetric ] 
( 0.0636610018750175252992355276057269, 
0.0636610018750175252992355276057269, 
0.269672331458315808034097805727606, 
) 
Corresponding weight: 
10 ^ -3 x 8.03571428571428571428571428571428,



Region: Simplex 
Dimension: 3 
Degree: 7 
Points: 31 
Structure: Fully symmetric 
Rule struct: 0 1 0 0 0 0 1 3 0 1 0 
Generator: [ Fully symmetric ] 
( 0.5, 0.5, 0., ) 
Corresponding weight: 
10 ^ -4 x 9.70017636684303350970017636684303,
Generator: [ Fully symmetric ] 
( 0.25, 0.25, 0.25, ) 
Corresponding weight: 
0.0182642234661088202912015685649462,

Generator: [ Fully symmetric ] 
( 0.0782131923303180643739942508375545, 
0.0782131923303180643739942508375545, 
0.0782131923303180643739942508375545, 
) 
Corresponding weight: 
0.0105999415244136869164138748545257,

Generator: [ Fully symmetric ] 
( 0.121843216663905174652156372684818, 
0.121843216663905174652156372684818, 
0.121843216663905174652156372684818, 
) 
Corresponding weight: 
-0.0625177401143318516914703474927900,

Generator: [ Fully symmetric ] 
( 0.332539164446420624152923823157707, 
0.332539164446420624152923823157707, 
0.332539164446420624152923823157707, 
) 
Corresponding weight: 
10 ^ -3 x 4.89142526307349938479576303671027,

Generator: [ Fully symmetric ] 
( 0.1, 0.1, 0.2, ) 
Corresponding weight: 
0.0275573192239858906525573192239858,



Region: Simplex 
Dimension: 3 
Degree: 8 
Points: 43 
Structure: Fully symmetric 
Rule struct: 0 0 0 0 0 0 1 3 1 2 0 
Generator: [ Fully symmetric ] 
( 0.25, 0.25, 0.25, ) 
Corresponding weight: 
-0.0205001886586399158405865177642941,
Generator: [ Fully symmetric ] 
( 0.206829931610673204083980900024961, 
0.206829931610673204083980900024961, 
0.206829931610673204083980900024961, 
) 
Corresponding weight: 
0.0142503058228669012484397415358704,

Generator: [ Fully symmetric ] 
( 0.0821035883105467230906058078714215, 
0.0821035883105467230906058078714215, 
0.0821035883105467230906058078714215, 
) 
Corresponding weight: 
10 ^ -3 x 1.96703331313390098756280342445466,

Generator: [ Fully symmetric ] 
( 10 ^ -3 x 5.78195050519799725317663886414270, 
10 ^ -3 x 5.78195050519799725317663886414270, 
10 ^ -3 x 5.78195050519799725317663886414270, 
) 
Corresponding weight: 
10 ^ -4 x 1.69834109092887379837744566704016,

Generator: [ Fully symmetric ] 
( 0.0505327400188942244256245285579071, 
0.0505327400188942244256245285579071, 
0.449467259981105775574375471442092, 
) 
Corresponding weight: 
10 ^ -3 x 4.57968382446728180074351446297276,

Generator: [ Fully symmetric ] 
( 0.229066536116811139600408854554753, 
0.229066536116811139600408854554753, 
0.0356395827885340437169173969506114, 
) 
Corresponding weight: 
10 ^ -3 x 5.70448580868191850680255862783040,

Generator: [ Fully symmetric ] 
( 0.0366077495531974236787738546327104, 
0.0366077495531974236787738546327104, 
0.190486041934633455699433285315099, 
) 
Corresponding weight: 
10 ^ -3 x 2.14051914116209259648335300092023,

Region: Simplex 
Dimension: 3 
Degree: 9 
Points: 53 
Structure: Fully symmetric 
Rule struct: 0 0 0 0 0 0 1 4 0 3 0 
Generator: [ Fully symmetric ] 
( 0.25, 0.25, 0.25, ) 
Corresponding weight: 
-0.137799038326108641245781827743502,
Generator: [ Fully symmetric ] 
( 0.0483510385497367408791638973335652, 
0.0483510385497367408791638973335652, 
0.0483510385497367408791638973335652, 
) 
Corresponding weight: 
10 ^ -3 x 1.86533656908528954751732956352791,

Generator: [ Fully symmetric ] 
( 0.324579280117882365858796772348839, 
0.324579280117882365858796772348839, 
0.324579280117882365858796772348839, 
) 
Corresponding weight: 
10 ^ -3 x 4.30942396949340069685481460143672,

Generator: [ Fully symmetric ] 
( 0.114616540223995219683790300317451, 
0.114616540223995219683790300317451, 
0.114616540223995219683790300317451, 
) 
Corresponding weight: 
-0.0901847664812015251273931145424214,

Generator: [ Fully symmetric ] 
( 0.225489951911513918847703826550363, 
0.225489951911513918847703826550363, 
0.225489951911513918847703826550363, 
) 
Corresponding weight: 
0.0446725762025114446937850368088403,

Generator: [ Fully symmetric ] 
( 0.131627809246869809838976179197822, 
0.131627809246869809838976179197822, 
0.0836647016171849679665385541714202, 
) 
Corresponding weight: 
0.0347004058845507618227002576570271,

Generator: [ Fully symmetric ] 
( 0.433951461411406772411426685800524, 
0.433951461411406772411426685800524, 
0.107769859549428611825218473500043, 
) 
Corresponding weight: 
10 ^ -3 x 3.35258390266064697010720483063291,

Generator: [ Fully symmetric ] 
( 10 ^ -3 x -1.37627731813820071002030321419028, 
10 ^ -3 x -1.37627731813820071002030321419028, 
0.276553472636807342120192082186261, 
) 
Corresponding weight: 
10 ^ -4 x 4.31628875556996929641889902726182,


 */
template<class Rd>
void GQuadratureFormular<Rd>::Verification()
{
  QF_exact<GQuadratureFormular<Rd>,Rd::d+1>(exact,this);
  const int d=Rd::d;
  const double tol=1.e-12;
  R err=0;
  R a[d+1],h;
  
  for (int k=0;k<=exact;  k++)
    {

      //  formule magic: int_K  \lamda^k =   d! k! / ( d+k)! 
      
      R sa[d+1];
      for(int l=0;l<=d;++l)	      
	sa[l]=0.;
     
      for (int j=0;j<n;j++)
	{
	  h = p[j];
	  Rd P = p[j];
	  
	  for(int l=0;l<d;++l)	      
	    a[l]=p[j][l];
	  a[d] = 1.-p[j].sum();
	  
	  for(int l=0;l<=d;++l)	      	  
	    sa[l]+=h*pow(a[l],k);  
	}
      
      R se(1),see(1);
      for (int i=1;i<=k;i++) 
	se *= (R) i / (R) (i+d);
      see=se;
      
      for(int l=0;l<=d;++l)
	err = Max(err,Abs(se-sa[l]));
      
      if (err>tol)
	{
	  cerr << " d= " << d << "T Ordre= " << k << " d!k!/(d+k)!= " << se << " " ;
	  for(int l=0;l<=d;++l)
	    cerr << sa[l] << " ";
	  cerr << " err= " << err << endl;
	}
    }
  
  if(err>tol)
    {
      cerr << "Erreur dans la formule d'integration d=" <<d  << " exact = " << exact 
	   << " Nb Point = " << n << endl;
      assert(0);
    }
    
}
/*
  from:
http://www.cs.kuleuven.be/~nines/research/ecf/mtables.html
 */
typedef GQuadraturePoint<R3>  PQP3;
typedef GQuadratureFormular<R3>  PQF3;
// 2 4 (formule 2)
// A.H. Stroud, Approximate calculation of multiple integrals, Prentice-Hall, Englewood Cliffs, N.J., 1971.
// JavaScript:formule('t3-2-4b')
PQP3  QF_TET_2[]=
  { 
     PQP3(R3(0.58541019662496845446137605030968, 0.138196601125010515179541316563436, 0.138196601125010515179541316563436), 0.25),
     PQP3(R3(0.138196601125010515179541316563436,  0.58541019662496845446137605030968, 0.138196601125010515179541316563436), 0.25),
     PQP3(R3(0.138196601125010515179541316563436, 0.138196601125010515179541316563436,  0.58541019662496845446137605030968),0.25),
     PQP3(R3(0.138196601125010515179541316563436, 0.138196601125010515179541316563436, 0.138196601125010515179541316563436),0.25)
  };
  PQF3 const QuadratureFormular_Tet_2(2,4,QF_TET_2);
// 5  14  (formule 1)
/*
GM78
A. Grundmann and H.M. M鰈ler, Invariant integration formulas for the n-simplex by combinatorial methods, SIAM J. Numer. Anal. 15 (1978), 282--290.
 */
PQP3  QF_TET_5[]=
   { 
     PQP3(R3(0.7217942490673263207930282587889082 , 0.0927352503108912264023239137370306 , 0.0927352503108912264023239137370306) , 0.0122488405193936582572850342477212*6.),
     PQP3(R3(0.0927352503108912264023239137370306 , 0.7217942490673263207930282587889082 , 0.0927352503108912264023239137370306) , 0.0122488405193936582572850342477212*6.),
     PQP3(R3(0.0927352503108912264023239137370306 , 0.0927352503108912264023239137370306 , 0.7217942490673263207930282587889082) , 0.0122488405193936582572850342477212*6.),
     PQP3(R3(0.0927352503108912264023239137370306 , 0.0927352503108912264023239137370306 , 0.0927352503108912264023239137370306) , 0.0122488405193936582572850342477212*6),
     PQP3(R3(0.067342242210098170607962798709629 , 0.310885919263300609797345733763457 , 0.310885919263300609797345733763457) , 0.0187813209530026417998642753888810*6.),
     PQP3(R3(0.310885919263300609797345733763457 , 0.067342242210098170607962798709629 , 0.310885919263300609797345733763457) , 0.0187813209530026417998642753888810*6.),
     PQP3(R3(0.310885919263300609797345733763457 , 0.310885919263300609797345733763457 , 0.067342242210098170607962798709629) , 0.0187813209530026417998642753888810*6.),
     PQP3(R3(0.310885919263300609797345733763457 , 0.310885919263300609797345733763457 , 0.310885919263300609797345733763457) , 0.0187813209530026417998642753888810*6.),
     PQP3(R3(0.454496295874350350508119473720660,  0.454496295874350350508119473720660, 0.045503704125649649491880526279339),  7.09100346284691107301157135337624E-3*6.),
     PQP3(R3(00.454496295874350350508119473720660,  0.045503704125649649491880526279339, 0.454496295874350350508119473720660), 7.09100346284691107301157135337624E-3*6.),
     PQP3(R3(00.045503704125649649491880526279339,  0.454496295874350350508119473720660, 0.454496295874350350508119473720660), 7.09100346284691107301157135337624E-3*6.),
     PQP3(R3(00.045503704125649649491880526279339,  0.045503704125649649491880526279339, 0.454496295874350350508119473720660), 7.09100346284691107301157135337624E-3*6.),
     PQP3(R3(00.045503704125649649491880526279339,  0.454496295874350350508119473720660, 0.045503704125649649491880526279339), 7.09100346284691107301157135337624E-3*6.),
     PQP3(R3(00.454496295874350350508119473720660,  0.045503704125649649491880526279339, 0.045503704125649649491880526279339), 7.09100346284691107301157135337624E-3*6.)
   };

  PQF3 const QuadratureFormular_Tet_5(5,14,QF_TET_5);

  template<> const GQuadratureFormular<R1> * GQuadratureFormular<R1>::Default= & QF_GaussLegendre3;
  template<> const GQuadratureFormular<R2> * GQuadratureFormular<R2>::Default= & QuadratureFormular_T_5;
  template<> const GQuadratureFormular<R3> * GQuadratureFormular<R3>::Default= & QuadratureFormular_Tet_5;

} // end namespace Fem2D