quadrature_gauss_2d.c

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

		{0.099520061958437e+00, 0.450239969020782e+00,                    0., 0.024350878353672e+00 / 2.0}, // 3-perm
		{0.199467521245206e+00, 0.400266239377397e+00,                    0., 0.031107550868969e+00 / 2.0}, // 3-perm
		{0.495717464058095e+00, 0.252141267970953e+00,                    0., 0.031257111218620e+00 / 2.0}, // 3-perm
		{0.675905990683077e+00, 0.162047004658461e+00,                    0., 0.024815654339665e+00 / 2.0}, // 3-perm
		{0.848248235478508e+00, 0.075875882260746e+00,                    0., 0.014056073070557e+00 / 2.0}, // 3-perm
		{0.968690546064356e+00, 0.015654726967822e+00,                    0., 0.003194676173779e+00 / 2.0}, // 3-perm
		{0.010186928826919e+00, 0.334319867363658e+00, 0.655493203809423e+00, 0.008119655318993e+00 / 2.0}, // 6-perm
		{0.135440871671036e+00, 0.292221537796944e+00, 0.572337590532020e+00, 0.026805742283163e+00 / 2.0}, // 6-perm
		{0.054423924290583e+00, 0.319574885423190e+00, 0.626001190286228e+00, 0.018459993210822e+00 / 2.0}, // 6-perm 
		{0.012868560833637e+00, 0.190704224192292e+00, 0.796427214974071e+00, 0.008476868534328e+00 / 2.0}, // 6-perm
		{0.067165782413524e+00, 0.180483211648746e+00, 0.752351005937729e+00, 0.018292796770025e+00 / 2.0}, // 6-perm 
		{0.014663182224828e+00, 0.080711313679564e+00, 0.904625504095608e+00, 0.006665632004165e+00 / 2.0}  // 6-perm
	      };

	      
	      // Now call the dunavant routine to generate _points and _weights
	      dunavant_rule(p, 15);

	      return;
	    }


	    
	    // Dunavant's 18th-order rule contains points outside the region and is therefore unsuitable
	    // for our FEM calculations.  His 19th-order rule has 73 points, compared with 100 points for
	    // a comparable-order conical product rule.
	    //
	    // It was copied 23rd June 2008 from:
	    // http://people.scs.fsu.edu/~burkardt/f_src/dunavant/dunavant.f90
	  case EIGHTTEENTH:  
	  case NINTEENTH:
	    {
	      _points.resize (73);
	      _weights.resize(73);

	      // The raw data for the quadrature rule.
	      const Real p[17][4] = {
		{                1./3.,                    0.,                    0., 0.032906331388919e+00 / 2.0}, // 1-perm
		{0.020780025853987e+00, 0.489609987073006e+00,                    0., 0.010330731891272e+00 / 2.0}, // 3-perm
		{0.090926214604215e+00, 0.454536892697893e+00,                    0., 0.022387247263016e+00 / 2.0}, // 3-perm
		{0.197166638701138e+00, 0.401416680649431e+00,                    0., 0.030266125869468e+00 / 2.0}, // 3-perm
		{0.488896691193805e+00, 0.255551654403098e+00,                    0., 0.030490967802198e+00 / 2.0}, // 3-perm
		{0.645844115695741e+00, 0.177077942152130e+00,                    0., 0.024159212741641e+00 / 2.0}, // 3-perm
		{0.779877893544096e+00, 0.110061053227952e+00,                    0., 0.016050803586801e+00 / 2.0}, // 3-perm
		{0.888942751496321e+00, 0.055528624251840e+00,                    0., 0.008084580261784e+00 / 2.0}, // 3-perm
		{0.974756272445543e+00, 0.012621863777229e+00,                    0., 0.002079362027485e+00 / 2.0}, // 3-perm
		{0.003611417848412e+00, 0.395754787356943e+00, 0.600633794794645e+00, 0.003884876904981e+00 / 2.0}, // 6-perm
		{0.134466754530780e+00, 0.307929983880436e+00, 0.557603261588784e+00, 0.025574160612022e+00 / 2.0}, // 6-perm
		{0.014446025776115e+00, 0.264566948406520e+00, 0.720987025817365e+00, 0.008880903573338e+00 / 2.0}, // 6-perm 
		{0.046933578838178e+00, 0.358539352205951e+00, 0.594527068955871e+00, 0.016124546761731e+00 / 2.0}, // 6-perm
		{0.002861120350567e+00, 0.157807405968595e+00, 0.839331473680839e+00, 0.002491941817491e+00 / 2.0}, // 6-perm 
		{0.223861424097916e+00, 0.075050596975911e+00, 0.701087978926173e+00, 0.018242840118951e+00 / 2.0}, // 6-perm
		{0.034647074816760e+00, 0.142421601113383e+00, 0.822931324069857e+00, 0.010258563736199e+00 / 2.0}, // 6-perm 
		{0.010161119296278e+00, 0.065494628082938e+00, 0.924344252620784e+00, 0.003799928855302e+00 / 2.0}  // 6-perm
	      };

	      
	      // Now call the dunavant routine to generate _points and _weights
	      dunavant_rule(p, 17);

	      return;
	    }

	    
	    // 20th-order rule by Wandzura.  
	    //
	    // Stephen Wandzura, Hong Xiao,
	    // Symmetric Quadrature Rules on a Triangle,
	    // Computers and Mathematics with Applications,
	    // Volume 45, Number 12, June 2003, pages 1829-1840.
	    //
	    // Wandzura's work extends the work of Dunavant by providing degree
	    // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
	    //
	    // Copied on 3rd July 2008 from:
	    // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
	  case TWENTIETH:
	    {
	      // The equivalent concial product rule would have 121 points
	      _points.resize (85);
	      _weights.resize(85);

	      // The raw data for the quadrature rule.
	      const Real p[19][4] = {
		{0.33333333333333e+00,                  0.0,                  0.0, 0.2761042699769952e-01 / 2.0}, // 1-perm
		{0.00150064932443e+00, 0.49924967533779e+00,                  0.0, 0.1779029547326740e-02 / 2.0}, // 3-perm 
		{0.09413975193895e+00, 0.45293012403052e+00,                  0.0, 0.2011239811396117e-01 / 2.0}, // 3-perm
		{0.20447212408953e+00, 0.39776393795524e+00,                  0.0, 0.2681784725933157e-01 / 2.0}, // 3-perm
		{0.47099959493443e+00, 0.26450020253279e+00,                  0.0, 0.2452313380150201e-01 / 2.0}, // 3-perm
		{0.57796207181585e+00, 0.21101896409208e+00,                  0.0, 0.1639457841069539e-01 / 2.0}, // 3-perm
		{0.78452878565746e+00, 0.10773560717127e+00,                  0.0, 0.1479590739864960e-01 / 2.0}, // 3-perm
		{0.92186182432439e+00, 0.03906908783780e+00,                  0.0, 0.4579282277704251e-02 / 2.0}, // 3-perm
		{0.97765124054134e+00, 0.01117437972933e+00,                  0.0, 0.1651826515576217e-02 / 2.0}, // 3-perm
		{0.00534961818734e+00, 0.06354966590835e+00, 0.93110071590431e+00, 0.2349170908575584e-02 / 2.0}, // 6-perm
		{0.00795481706620e+00, 0.15710691894071e+00, 0.83493826399309e+00, 0.4465925754181793e-02 / 2.0}, // 6-perm
		{0.01042239828126e+00, 0.39564211436437e+00, 0.59393548735436e+00, 0.6099566807907972e-02 / 2.0}, // 6-perm
		{0.01096441479612e+00, 0.27316757071291e+00, 0.71586801449097e+00, 0.6891081327188203e-02 / 2.0}, // 6-perm
		{0.03856671208546e+00, 0.10178538248502e+00, 0.85964790542952e+00, 0.7997475072478163e-02 / 2.0}, // 6-perm
		{0.03558050781722e+00, 0.44665854917641e+00, 0.51776094300637e+00, 0.7386134285336024e-02 / 2.0}, // 6-perm
		{0.04967081636276e+00, 0.19901079414950e+00, 0.75131838948773e+00, 0.1279933187864826e-01 / 2.0}, // 6-perm
		{0.05851972508433e+00, 0.32426118369228e+00, 0.61721909122339e+00, 0.1725807117569655e-01 / 2.0}, // 6-perm
		{0.12149778700439e+00, 0.20853136321013e+00, 0.66997084978547e+00, 0.1867294590293547e-01 / 2.0}, // 6-perm
		{0.14071084494394e+00, 0.32317056653626e+00, 0.53611858851980e+00, 0.2281822405839526e-01 / 2.0}  // 6-perm
	      };

	      
	      // Now call the dunavant routine to generate _points and _weights
	      dunavant_rule(p, 19);

	      return;
	    }



	    // 25th-order rule by Wandzura.  
	    //
	    // Stephen Wandzura, Hong Xiao,
	    // Symmetric Quadrature Rules on a Triangle,
	    // Computers and Mathematics with Applications,
	    // Volume 45, Number 12, June 2003, pages 1829-1840.
	    //
	    // Wandzura's work extends the work of Dunavant by providing degree
	    // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
	    //
	    // Copied on 3rd July 2008 from:
	    // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
	    // case TWENTYFIRST: // fall through to 121 point conical product rule below
	  case TWENTYSECOND:  
	  case TWENTYTHIRD:   
	  case TWENTYFOURTH:  	    
	  case TWENTYFIFTH:
	    {
	      // The equivalent concial product rule would have 169 points
	      _points.resize (126);
	      _weights.resize(126);

	      // The raw data for the quadrature rule.
	      const Real p[26][4] = {
		{0.02794648307317e+00, 0.48602675846341e+00,                  0.0, 0.8005581880020417e-02 / 2.0},  // 3-perm
		{0.13117860132765e+00, 0.43441069933617e+00,                  0.0, 0.1594707683239050e-01 / 2.0},  // 3-perm
		{0.22022172951207e+00, 0.38988913524396e+00,                  0.0, 0.1310914123079553e-01 / 2.0},  // 3-perm
		{0.40311353196039e+00, 0.29844323401980e+00,                  0.0, 0.1958300096563562e-01 / 2.0},  // 3-perm
		{0.53191165532526e+00, 0.23404417233737e+00,                  0.0, 0.1647088544153727e-01 / 2.0},  // 3-perm
		{0.69706333078196e+00, 0.15146833460902e+00,                  0.0, 0.8547279074092100e-02 / 2.0},  // 3-perm
		{0.77453221290801e+00, 0.11273389354599e+00,                  0.0, 0.8161885857226492e-02 / 2.0},  // 3-perm
		{0.84456861581695e+00, 0.07771569209153e+00,                  0.0, 0.6121146539983779e-02 / 2.0},  // 3-perm
		{0.93021381277141e+00, 0.03489309361430e+00,                  0.0, 0.2908498264936665e-02 / 2.0},  // 3-perm
		{0.98548363075813e+00, 0.00725818462093e+00,                  0.0, 0.6922752456619963e-03 / 2.0},  // 3-perm
		{0.00129235270444e+00, 0.22721445215336e+00, 0.77149319514219e+00, 0.1248289199277397e-02 / 2.0},  // 6-perm
		{0.00539970127212e+00, 0.43501055485357e+00, 0.55958974387431e+00, 0.3404752908803022e-02 / 2.0},  // 6-perm
		{0.00638400303398e+00, 0.32030959927220e+00, 0.67330639769382e+00, 0.3359654326064051e-02 / 2.0},  // 6-perm
		{0.00502821150199e+00, 0.09175032228001e+00, 0.90322146621800e+00, 0.1716156539496754e-02 / 2.0},  // 6-perm
		{0.00682675862178e+00, 0.03801083585872e+00, 0.95516240551949e+00, 0.1480856316715606e-02 / 2.0},  // 6-perm
		{0.01001619963993e+00, 0.15742521848531e+00, 0.83255858187476e+00, 0.3511312610728685e-02 / 2.0},  // 6-perm
		{0.02575781317339e+00, 0.23988965977853e+00, 0.73435252704808e+00, 0.7393550149706484e-02 / 2.0},  // 6-perm
		{0.03022789811992e+00, 0.36194311812606e+00, 0.60782898375402e+00, 0.7983087477376558e-02 / 2.0},  // 6-perm
		{0.03050499010716e+00, 0.08355196095483e+00, 0.88594304893801e+00, 0.4355962613158041e-02 / 2.0},  // 6-perm
		{0.04595654736257e+00, 0.14844322073242e+00, 0.80560023190501e+00, 0.7365056701417832e-02 / 2.0},  // 6-perm
		{0.06744280054028e+00, 0.28373970872753e+00, 0.64881749073219e+00, 0.1096357284641955e-01 / 2.0},  // 6-perm
		{0.07004509141591e+00, 0.40689937511879e+00, 0.52305553346530e+00, 0.1174996174354112e-01 / 2.0},  // 6-perm
		{0.08391152464012e+00, 0.19411398702489e+00, 0.72197448833499e+00, 0.1001560071379857e-01 / 2.0},  // 6-perm
		{0.12037553567715e+00, 0.32413434700070e+00, 0.55549011732214e+00, 0.1330964078762868e-01 / 2.0},  // 6-perm
		{0.14806689915737e+00, 0.22927748355598e+00, 0.62265561728665e+00, 0.1415444650522614e-01 / 2.0},  // 6-perm
		{0.19177186586733e+00, 0.32561812259598e+00, 0.48261001153669e+00, 0.1488137956116801e-01 / 2.0}   // 6-perm
	      };

	      
	      // Now call the dunavant routine to generate _points and _weights
	      dunavant_rule(p, 26);

	      return;
	    }



	    // 30th-order rule by Wandzura.  
	    //
	    // Stephen Wandzura, Hong Xiao,
	    // Symmetric Quadrature Rules on a Triangle,
	    // Computers and Mathematics with Applications,
	    // Volume 45, Number 12, June 2003, pages 1829-1840.
	    //
	    // Wandzura's work extends the work of Dunavant by providing degree
	    // 5,10,15,20,25, and 30 rules with positive weights for the triangle.
	    //
	    // Copied on 3rd July 2008 from:
	    // http://people.scs.fsu.edu/~burkardt/f_src/wandzura/wandzura.f90
	  case TWENTYSIXTH:   
	  case TWENTYSEVENTH: 
	  case TWENTYEIGHTH:  
	  case TWENTYNINTH:   	    
	  case THIRTIETH:
	    {
	      // The equivalent concial product rule would have 256 points
	      _points.resize (175);
	      _weights.resize(175);

	      // The raw data for the quadrature rule.
	      const Real p[36][4] = {
		{0.33333333333333e+00,                  0.0,                  0.0, 0.1557996020289920e-01 / 2.0}, // 1-perm 
		{0.00733011643277e+00, 0.49633494178362e+00,                  0.0, 0.3177233700534134e-02 / 2.0}, // 3-perm  
		{0.08299567580296e+00, 0.45850216209852e+00,                  0.0, 0.1048342663573077e-01 / 2.0}, // 3-perm  
		{0.15098095612541e+00, 0.42450952193729e+00,                  0.0, 0.1320945957774363e-01 / 2.0}, // 3-perm  
		{0.23590585989217e+00, 0.38204707005392e+00,                  0.0, 0.1497500696627150e-01 / 2.0}, // 3-perm  
		{0.43802430840785e+00, 0.28098784579608e+00,                  0.0, 0.1498790444338419e-01 / 2.0}, // 3-perm  
		{0.54530204829193e+00, 0.22734897585403e+00,                  0.0, 0.1333886474102166e-01 / 2.0}, // 3-perm  
		{0.65088177698254e+00, 0.17455911150873e+00,                  0.0, 0.1088917111390201e-01 / 2.0}, // 3-perm  
		{0.75348314559713e+00, 0.12325842720144e+00,                  0.0, 0.8189440660893461e-02 / 2.0}, // 3-perm  
		{0.83983154221561e+00, 0.08008422889220e+00,                  0.0, 0.5575387588607785e-02 / 2.0}, // 3-perm  
		{0.90445106518420e+00, 0.04777446740790e+00,                  0.0, 0.3191216473411976e-02 / 2.0}, // 3-perm  
		{0.95655897063972e+00, 0.02172051468014e+00,                  0.0, 0.1296715144327045e-02 / 2.0}, // 3-perm  
		{0.99047064476913e+00, 0.00476467761544e+00,                  0.0, 0.2982628261349172e-03 / 2.0}, // 3-perm  
		{0.00092537119335e+00, 0.41529527091331e+00, 0.58377935789334e+00, 0.9989056850788964e-03 / 2.0}, // 6-perm  
		{0.00138592585556e+00, 0.06118990978535e+00, 0.93742416435909e+00, 0.4628508491732533e-03 / 2.0}, // 6-perm  
		{0.00368241545591e+00, 0.16490869013691e+00, 0.83140889440718e+00, 0.1234451336382413e-02 / 2.0}, // 6-perm  
		{0.00390322342416e+00, 0.02503506223200e+00, 0.97106171434384e+00, 0.5707198522432062e-03 / 2.0}, // 6-perm  
		{0.00323324815501e+00, 0.30606446515110e+00, 0.69070228669389e+00, 0.1126946125877624e-02 / 2.0}, // 6-perm  
		{0.00646743211224e+00, 0.10707328373022e+00, 0.88645928415754e+00, 0.1747866949407337e-02 / 2.0}, // 6-perm  
		{0.00324747549133e+00, 0.22995754934558e+00, 0.76679497516308e+00, 0.1182818815031657e-02 / 2.0}, // 6-perm  
		{0.00867509080675e+00, 0.33703663330578e+00, 0.65428827588746e+00, 0.1990839294675034e-02 / 2.0}, // 6-perm  
		{0.01559702646731e+00, 0.05625657618206e+00, 0.92814639735063e+00, 0.1900412795035980e-02 / 2.0}, // 6-perm  
		{0.01797672125369e+00, 0.40245137521240e+00, 0.57957190353391e+00, 0.4498365808817451e-02 / 2.0}, // 6-perm  
		{0.01712424535389e+00, 0.24365470201083e+00, 0.73922105263528e+00, 0.3478719460274719e-02 / 2.0}, // 6-perm  
		{0.02288340534658e+00, 0.16538958561453e+00, 0.81172700903888e+00, 0.4102399036723953e-02 / 2.0}, // 6-perm  
		{0.03273759728777e+00, 0.09930187449585e+00, 0.86796052821639e+00, 0.4021761549744162e-02 / 2.0}, // 6-perm  
		{0.03382101234234e+00, 0.30847833306905e+00, 0.65770065458860e+00, 0.6033164660795066e-02 / 2.0}, // 6-perm  
		{0.03554761446002e+00, 0.46066831859211e+00, 0.50378406694787e+00, 0.3946290302129598e-02 / 2.0}, // 6-perm  
		{0.05053979030687e+00, 0.21881529945393e+00, 0.73064491023920e+00, 0.6644044537680268e-02 / 2.0}, // 6-perm  
		{0.05701471491573e+00, 0.37920955156027e+00, 0.56377573352399e+00, 0.8254305856078458e-02 / 2.0}, // 6-perm  
		{0.06415280642120e+00, 0.14296081941819e+00, 0.79288637416061e+00, 0.6496056633406411e-02 / 2.0}, // 6-perm  
		{0.08050114828763e+00, 0.28373128210592e+00, 0.63576756960645e+00, 0.9252778144146602e-02 / 2.0}, // 6-perm  
		{0.10436706813453e+00, 0.19673744100444e+00, 0.69889549086103e+00, 0.9164920726294280e-02 / 2.0}, // 6-perm  
		{0.11384489442875e+00, 0.35588914121166e+00, 0.53026596435959e+00, 0.1156952462809767e-01 / 2.0}, // 6-perm  
		{0.14536348771552e+00, 0.25981868535191e+00, 0.59481782693256e+00, 0.1176111646760917e-01 / 2.0}, // 6-perm  
		{0.18994565282198e+00, 0.32192318123130e+00, 0.48813116594672e+00, 0.1382470218216540e-01 / 2.0}  // 6-perm   
	      };

	      
	      // Now call the dunavant routine to generate _points and _weights
	      dunavant_rule(p, 36);

	      return;
	    }
	    
	    
	    // By default, we fall back on the conical product rules.  If the user
	    // requests an order higher than what is currently available in the 1D
	    // rules, an error will be thrown from the respective 1D code.
	  default:
	    {
	      // 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 a
	      // 1D Jacobi-Gauss rule on [0,1].

	      // Define the quadrature rules...
              // FIXME - why can't we init() these explicitly? [RHS]
	      QGauss  gauss1D(1,static_cast<Order>(_order+2*p));
	      QJacobi jac1D(1,static_cast<Order>(_order+2*p),1,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_tri(gauss1D, jac1D);
	      return;
	    }
	    
	  }
      }

	    
      //---------------------------------------------
      // Unsupported type
    default:
      {
	std::cerr << "Element type not supported!:" << _type << std::endl;
	libmesh_error();
      }
    }

  libmesh_error();

  return;

#endif
}