|
rbmatlab 0.10.01
|
general/geometry/triaquadrature.m
00001 function res = triaquadrature(poldeg,func,varargin) 00002 %function res = triaquadrature(poldeg,func,varargin) 00003 % 00004 % integration of function func over reference triangle. 00005 % by Gaussian quadrature exactly integrating polynoms of poldeg. 00006 % func is a function getting a local coordinate vector and varargin 00007 % and giving a (vectorial or scalar) result 00008 % formulae from: 00009 % G.R. Cowper "Gaussian Quadrature Formulas for Triangles" 00010 % and dune-quadratures quadratures.cc 00011 % degrees 1-11 are copied, see the quadratures.cc included below 00012 % for more quadraturerules. 00013 00014 % Bernard Haasdonk 28.1.2009 00015 00016 switch poldeg 00017 case 0 00018 %npoints = 1; 00019 points = 1/3 * ones(1,2); 00020 weights = 0.5; 00021 case {1,2} 00022 %npoints = 3; 00023 points = [2/3, 1/6; 00024 1/6, 2/3; 00025 1/6, 1/6]; 00026 weights = 1/6 * ones(1,3); 00027 case {3,4} 00028 %npoints = 6; 00029 points = [ ... 00030 0.81684757298045851308085707319560, 0.091576213509770743459571463402202; 00031 0.091576213509770743459571463402202, 0.81684757298045851308085707319560; 00032 0.091576213509770743459571463402202, 0.091576213509770743459571463402202; 00033 0.10810301816807022736334149223390, 0.44594849091596488631832925388305; 00034 0.44594849091596488631832925388305, 0.10810301816807022736334149223390; 00035 0.44594849091596488631832925388305, 0.44594849091596488631832925388305; 00036 ]; 00037 weights = [ 0.5 * 0.10995174365532186763832632490021 * ones(3,1);... 00038 0.5 * 0.22338158967801146569500700843312 * ones(3,1)]; 00039 00040 % points = [0.816847572980459, 0.091576213509771, 0.091576213509771; ... 00041 % 0.091576213509771, 0.816847572980459, 0.091576213509771; ... 00042 % 0.091576213509771, 0.091576213509771, 0.816847572980459; ... 00043 % 0.108103018168070, 0.445948490915965, 0.445948490915965; ... 00044 % 0.445948490915965, 0.108103018168070, 0.445948490915965; ... 00045 % 0.445948490915965, 0.445948490915965, 0.108103018168070]'; 00046 % weights = 0.5 * [0.109951743655322, 0.109951743655322, 0.109951743655322, ... 00047 % 0.223381589678011, 0.223381589678011, ... 00048 % 0.223381589678011]; 00049 case {5,6,7} 00050 % npoints = 13; 00051 % points = [1/3, 1/3, 1/3;... 00052 % 0.479308067841923, 0.260345966079038, 0.260345966079038; ... 00053 % 0.260345966079038, 0.479308067841923, 0.260345966079038; ... 00054 % 0.260345966079038, 0.260345966079038, 0.479308067841923; ... 00055 % 0.869739794195568, 0.065130102902216, 0.065130102902216; ... 00056 % 0.065130102902216, 0.869739794195568, 0.065130102902216; ... 00057 % 0.065130102902216, 0.065130102902216, 0.869739794195568; ... 00058 % 0.638444188569809, 0.312865496004875, 0.048690315425316;... 00059 % 0.638444188569809, 0.048690315425316, 0.312865496004875;... 00060 % 0.048690315425316, 0.638444188569809, 0.312865496004875;... 00061 % 0.048690315425316, 0.312865496004875, 0.638444188569809;... 00062 % 0.312865496004875, 0.048690315425316, 0.638444188569809;... 00063 % 0.312865496004875, 0.638444188569809, 0.048690315425316 ... 00064 % ]'; 00065 % weights = 0.5 * [-0.149570044467670, ... 00066 % 0.175615257433204, 0.175615257433204, 0.175615257433204, ... 00067 % 0.053347235608839, 0.053347235608839, 0.053347235608839, ... 00068 % 0.077113760890257, 0.077113760890257, 0.077113760890257, ... 00069 % 0.077113760890257, 0.077113760890257, 0.077113760890257 ... 00070 % ]; 00071 %npoints = 12; 00072 points = [... 00073 0.0623822650944021181736830009963499, 0.0675178670739160854425571310508685; 00074 0.0675178670739160854425571310508685, 0.870099867831681796383759867952782; 00075 0.870099867831681796383759867952782, 0.0623822650944021181736830009963499; 00076 0.0552254566569266117374791902756449, 0.321502493851981822666307849199202; 00077 0.321502493851981822666307849199202, 0.623272049491091565596212960525153; 00078 0.623272049491091565596212960525153, 0.0552254566569266117374791902756449; 00079 0.0343243029450971464696306424839376,0.660949196186735657611980310197799; 00080 0.660949196186735657611980310197799, 0.304726500868167195918389047318263; 00081 0.304726500868167195918389047318263, 0.0343243029450971464696306424839376; 00082 0.515842334353591779257463386826430, 0.277716166976391782569581871393723; 00083 0.277716166976391782569581871393723, 0.20644149867001643817295474177985; 00084 0.20644149867001643817295474177985, 0.515842334353591779257463386826430]; 00085 00086 weights = 0.5 * [ 00087 0.053034056314872502857508360921478; 00088 0.053034056314872502857508360921478; 00089 0.053034056314872502857508360921478; 00090 0.087762817428892110073539806278575; 00091 0.087762817428892110073539806278575; 00092 0.087762817428892110073539806278575; 00093 0.057550085569963171476890993800437; 00094 0.057550085569963171476890993800437; 00095 0.057550085569963171476890993800437; 00096 0.13498637401960554892539417233284; 00097 0.13498637401960554892539417233284; 00098 0.13498637401960554892539417233284]; 00099 case {8,9,10,11} 00100 %npoints = 28; 00101 points = [ 00102 0.858870281282636704039173938058347, 0.141129718717363295960826061941652; 00103 0.858870281282636704039173938058347, 0.0; 00104 0.141129718717363295960826061941652, 0.858870281282636704039173938058347; 00105 0.141129718717363295960826061941652, 0.0; 00106 0.0, 0.858870281282636704039173938058347; 00107 0.0, 0.141129718717363295960826061941652; 00108 0.333333333333333333333333333333333, 0.333333333333333333333333333333333; 00109 0.025989140928287395260032485498841, 0.025989140928287395260032485498841; 00110 0.025989140928287395260032485498841, 0.94802171814342520947993502900232; 00111 0.94802171814342520947993502900232, 0.025989140928287395260032485498841; 00112 0.094287502647922495630569776275405, 0.094287502647922495630569776275405; 00113 0.094287502647922495630569776275405, 0.81142499470415500873886044744919; 00114 0.81142499470415500873886044744919, 0.094287502647922495630569776275405; 00115 0.49463677501721381374163260230644, 0.49463677501721381374163260230644; 00116 0.49463677501721381374163260230644, 0.01072644996557237251673479538713; 00117 0.01072644996557237251673479538713, 0.49463677501721381374163260230644; 00118 0.20734338261451133345293402411297, 0.20734338261451133345293402411297; 00119 0.20734338261451133345293402411297, 0.58531323477097733309413195177407; 00120 0.58531323477097733309413195177407, 0.20734338261451133345293402411297; 00121 0.43890780570049209506106538163613, 0.43890780570049209506106538163613; 00122 0.43890780570049209506106538163613, 0.12218438859901580987786923672775; 00123 0.12218438859901580987786923672775, 0.43890780570049209506106538163613; 00124 0.67793765488259040154212614118875, 0.044841677589130443309052391468801; 00125 0.67793765488259040154212614118875, 0.27722066752827915514882146734245; 00126 0.044841677589130443309052391468801, 0.67793765488259040154212614118875; 00127 0.044841677589130443309052391468801, 0.27722066752827915514882146734245; 00128 0.27722066752827915514882146734245, 0.67793765488259040154212614118875; 00129 0.27722066752827915514882146734245, 0.044841677589130443309052391468801]; 00130 00131 weights = 0.5 * [... 00132 0.0073623837833005542642588950473806; 00133 0.0073623837833005542642588950473806; 00134 0.0073623837833005542642588950473806; 00135 0.0073623837833005542642588950473806; 00136 0.0073623837833005542642588950473806; 00137 0.0073623837833005542642588950473806; 00138 0.087977301162232238798093169321456; 00139 0.0087443115537360230495164287998252; 00140 0.0087443115537360230495164287998252; 00141 0.0087443115537360230495164287998252; 00142 0.038081571993934937515024339435614; 00143 0.038081571993934937515024339435614; 00144 0.038081571993934937515024339435614; 00145 0.018855448056131292058476782591115; 00146 0.018855448056131292058476782591115; 00147 0.018855448056131292058476782591115; 00148 0.072159697544739526124029988586463; 00149 0.072159697544739526124029988586463; 00150 0.072159697544739526124029988586463; 00151 0.069329138705535899841765650903814; 00152 0.069329138705535899841765650903814; 00153 0.069329138705535899841765650903814; 00154 0.041056315429288566641652314907294; 00155 0.041056315429288566641652314907294; 00156 0.041056315429288566641652314907294; 00157 0.041056315429288566641652314907294; 00158 0.041056315429288566641652314907294; 00159 0.041056315429288566641652314907294]; 00160 otherwise 00161 error('quadrature of desired degree not available') 00162 end; 00163 00164 res = quadrature(weights,points,func,varargin{:}); 00165 00166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00167 % quadraturerules.cc from dune: 00168 % $$$ 00169 % $$$ #include "config.h" 00170 % $$$ #include "../quadraturerules.hh" 00171 % $$$ 00172 % $$$ namespace Dune { 00173 % $$$ 00174 % $$$ 00175 % $$$ template<int dim> 00176 % $$$ class SimplexQuadraturePoints; 00177 % $$$ 00178 % $$$ 00179 % $$$ template<int dim> 00180 % $$$ struct SimplexQuadraturePointsSingleton { 00181 % $$$ static SimplexQuadraturePoints<dim> sqp; 00182 % $$$ }; 00183 % $$$ 00184 % $$$ template<> 00185 % $$$ struct SimplexQuadraturePointsSingleton<2> { 00186 % $$$ static SimplexQuadraturePoints<2> sqp; 00187 % $$$ }; 00188 % $$$ 00189 % $$$ template<> 00190 % $$$ struct SimplexQuadraturePointsSingleton<3> { 00191 % $$$ static SimplexQuadraturePoints<3> sqp; 00192 % $$$ }; 00193 % $$$ 00194 % $$$ template<> 00195 % $$$ class SimplexQuadraturePoints<2> 00196 % $$$ { 00197 % $$$ public: 00198 % $$$ enum { MAXP=33}; 00199 % $$$ enum { highest_order=12 }; 00200 % $$$ 00201 % $$$ 00202 % $$$ SimplexQuadraturePoints () 00203 % $$$ { 00204 % $$$ int m = 0; 00205 % $$$ O[m] = 0; 00206 % $$$ 00207 % $$$ // polynom degree 1 00208 % $$$ m = 1; 00209 % $$$ G[m][0][0] = 0.333333333333333333333333333333333; 00210 % $$$ G[m][0][1] = 0.333333333333333333333333333333333; 00211 % $$$ W[m][0] = 0.5; 00212 % $$$ O[m] = 1; 00213 % $$$ 00214 % $$$ // polynom degree 2 00215 % $$$ // symmetric 00216 % $$$ m = 3; 00217 % $$$ G[m][0][0] = 4.0/6.0; 00218 % $$$ G[m][0][1] = 1.0/6.0; 00219 % $$$ G[m][1][0] = 1.0/6.0; 00220 % $$$ G[m][1][1] = 4.0/6.0; 00221 % $$$ G[m][2][0] = 1.0/6.0; 00222 % $$$ G[m][2][1] = 1.0/6.0; 00223 % $$$ W[m][0] = 0.5/3.0; 00224 % $$$ W[m][1] = 0.5/3.0; 00225 % $$$ W[m][2] = 0.5/3.0; 00226 % $$$ O[m] = 2; 00227 % $$$ 00228 % $$$ // polynom degree 3 00229 % $$$ // symmetric 00230 % $$$ m = 4; 00231 % $$$ G[m][0][0] = 10.0/30.0; 00232 % $$$ G[m][0][1] = 10.0/30.0; 00233 % $$$ G[m][1][0] = 18.0/30.0; 00234 % $$$ G[m][1][1] = 6.0/30.0; 00235 % $$$ G[m][2][0] = 6.0/30.0; 00236 % $$$ G[m][2][1] = 18.0/30.0; 00237 % $$$ G[m][3][0] = 6.0/30.0; 00238 % $$$ G[m][3][1] = 6.0/30.0; 00239 % $$$ 00240 % $$$ W[m][0] = 0.5 * -27.0/48.0; 00241 % $$$ W[m][1] = 0.5 * 25.0/48.0; 00242 % $$$ W[m][2] = 0.5 * 25.0/48.0; 00243 % $$$ W[m][3] = 0.5 * 25.0/48.0; 00244 % $$$ O[m] = 3; 00245 % $$$ 00246 % $$$ // polynomial degree 4 00247 % $$$ // symmetric points 00248 % $$$ m = 6; 00249 % $$$ G[m][0][0] = 0.81684757298045851308085707319560; 00250 % $$$ G[m][0][1] = 0.091576213509770743459571463402202; 00251 % $$$ G[m][1][0] = 0.091576213509770743459571463402202; 00252 % $$$ G[m][1][1] = 0.81684757298045851308085707319560; 00253 % $$$ G[m][2][0] = 0.091576213509770743459571463402202; 00254 % $$$ G[m][2][1] = 0.091576213509770743459571463402202; 00255 % $$$ G[m][3][0] = 0.10810301816807022736334149223390; 00256 % $$$ G[m][3][1] = 0.44594849091596488631832925388305; 00257 % $$$ G[m][4][0] = 0.44594849091596488631832925388305; 00258 % $$$ G[m][4][1] = 0.10810301816807022736334149223390; 00259 % $$$ G[m][5][0] = 0.44594849091596488631832925388305; 00260 % $$$ G[m][5][1] = 0.44594849091596488631832925388305; 00261 % $$$ 00262 % $$$ W[m][0] = 0.5 * 0.10995174365532186763832632490021; 00263 % $$$ W[m][1] = 0.5 * 0.10995174365532186763832632490021; 00264 % $$$ W[m][2] = 0.5 * 0.10995174365532186763832632490021; 00265 % $$$ W[m][3] = 0.5 * 0.22338158967801146569500700843312; 00266 % $$$ W[m][4] = 0.5 * 0.22338158967801146569500700843312; 00267 % $$$ W[m][5] = 0.5 * 0.22338158967801146569500700843312; 00268 % $$$ O[m] = 4; 00269 % $$$ 00270 % $$$ // polynomial degree 5 00271 % $$$ // symmetric points 00272 % $$$ 00273 % $$$ m = 7; 00274 % $$$ G[m][0][0] = 0.333333333333333333333333333333333; 00275 % $$$ G[m][0][1] = 0.333333333333333333333333333333333; 00276 % $$$ G[m][1][0] = 0.79742698535308732239802527616975; 00277 % $$$ G[m][1][1] = 0.1012865073234563388009873619151; 00278 % $$$ G[m][2][0] = 0.10128650732345633880098736191512; 00279 % $$$ G[m][2][1] = 0.79742698535308732239802527616975; 00280 % $$$ G[m][3][0] = 0.10128650732345633880098736191512; 00281 % $$$ G[m][3][1] = 0.10128650732345633880098736191512; 00282 % $$$ G[m][4][0] = 0.05971587178976982045911758097311; 00283 % $$$ G[m][4][1] = 0.47014206410511508977044120951345; 00284 % $$$ G[m][5][0] = 0.47014206410511508977044120951345; 00285 % $$$ G[m][5][1] = 0.05971587178976982045911758097311; 00286 % $$$ G[m][6][0] = 0.47014206410511508977044120951345; 00287 % $$$ G[m][6][1] = 0.47014206410511508977044120951345; 00288 % $$$ 00289 % $$$ W[m][0] = 0.5 * 0.225; 00290 % $$$ W[m][1] = 0.5 * 0.12593918054482715259568394550018; 00291 % $$$ W[m][2] = 0.5 * 0.12593918054482715259568394550018; 00292 % $$$ W[m][3] = 0.5 * 0.12593918054482715259568394550018; 00293 % $$$ W[m][4] = 0.5 * 0.13239415278850618073764938783315; 00294 % $$$ W[m][5] = 0.5 * 0.13239415278850618073764938783315; 00295 % $$$ W[m][6] = 0.5 * 0.13239415278850618073764938783315; 00296 % $$$ O[m] = 5; 00297 % $$$ 00298 % $$$ // polynomial degree 6 00299 % $$$ /* 12 inner Gauss points, positive weights */ 00300 % $$$ m=12; 00301 % $$$ G[m][0][0] = 0.063089014491502228340331602870819; 00302 % $$$ G[m][0][1] = 0.063089014491502228340331602870819; 00303 % $$$ G[m][1][0] = 0.063089014491502228340331602870819; 00304 % $$$ G[m][1][1] = 0.87382197101699554331933679425836; 00305 % $$$ G[m][2][0] = 0.87382197101699554331933679425836; 00306 % $$$ G[m][2][1] = 0.063089014491502228340331602870819; 00307 % $$$ G[m][3][0] = 0.24928674517091042129163855310702; 00308 % $$$ G[m][3][1] = 0.24928674517091042129163855310702; 00309 % $$$ G[m][4][0] = 0.24928674517091042129163855310702; 00310 % $$$ G[m][4][1] = 0.50142650965817915741672289378596; 00311 % $$$ G[m][5][0] = 0.50142650965817915741672289378596; 00312 % $$$ G[m][5][1] = 0.24928674517091042129163855310702; 00313 % $$$ G[m][6][0] = 0.053145049844816947353249671631398; 00314 % $$$ G[m][6][1] = 0.31035245103378440541660773395655; 00315 % $$$ G[m][7][0] = 0.053145049844816947353249671631398; 00316 % $$$ G[m][7][1] = 0.63650249912139864723014259441205; 00317 % $$$ G[m][8][0] = 0.31035245103378440541660773395655; 00318 % $$$ G[m][8][1] = 0.053145049844816947353249671631398; 00319 % $$$ G[m][9][0] = 0.31035245103378440541660773395655; 00320 % $$$ G[m][9][1] = 0.63650249912139864723014259441205; 00321 % $$$ G[m][10][0] = 0.63650249912139864723014259441205; 00322 % $$$ G[m][10][1] = 0.053145049844816947353249671631398; 00323 % $$$ G[m][11][0] = 0.63650249912139864723014259441205; 00324 % $$$ G[m][11][1] = 0.31035245103378440541660773395655; 00325 % $$$ 00326 % $$$ W[m][0] = 0.5 * 0.050844906370206816920936809106869; 00327 % $$$ W[m][1] = 0.5 * 0.050844906370206816920936809106869; 00328 % $$$ W[m][2] = 0.5 * 0.050844906370206816920936809106869; 00329 % $$$ W[m][3] = 0.5 * 0.11678627572637936602528961138558; 00330 % $$$ W[m][4] = 0.5 * 0.11678627572637936602528961138558; 00331 % $$$ W[m][5] = 0.5 * 0.11678627572637936602528961138558; 00332 % $$$ W[m][6] = 0.5 * 0.082851075618373575193553456420442; 00333 % $$$ W[m][7] = 0.5 * 0.082851075618373575193553456420442; 00334 % $$$ W[m][8] = 0.5 * 0.082851075618373575193553456420442; 00335 % $$$ W[m][9] = 0.5 * 0.082851075618373575193553456420442; 00336 % $$$ W[m][10] = 0.5 * 0.082851075618373575193553456420442; 00337 % $$$ W[m][11] = 0.5 * 0.082851075618373575193553456420442; 00338 % $$$ O[m] = 6; 00339 % $$$ 00340 % $$$ // polynomial degree 7 00341 % $$$ /* 12 inner Gauss points, positive weights */ 00342 % $$$ m=12; 00343 % $$$ G[m][0][0] = 0.0623822650944021181736830009963499; 00344 % $$$ G[m][0][1] = 0.0675178670739160854425571310508685; 00345 % $$$ G[m][1][0] = 0.0675178670739160854425571310508685; 00346 % $$$ G[m][1][1] = 0.870099867831681796383759867952782; 00347 % $$$ G[m][2][0] = 0.870099867831681796383759867952782; 00348 % $$$ G[m][2][1] = 0.0623822650944021181736830009963499; 00349 % $$$ G[m][3][0] = 0.0552254566569266117374791902756449; 00350 % $$$ G[m][3][1] = 0.321502493851981822666307849199202; 00351 % $$$ G[m][4][0] = 0.321502493851981822666307849199202; 00352 % $$$ G[m][4][1] = 0.623272049491091565596212960525153; 00353 % $$$ G[m][5][0] = 0.623272049491091565596212960525153; 00354 % $$$ G[m][5][1] = 0.0552254566569266117374791902756449; 00355 % $$$ G[m][6][0] = 0.0343243029450971464696306424839376; 00356 % $$$ G[m][6][1] = 0.660949196186735657611980310197799; 00357 % $$$ G[m][7][0] = 0.660949196186735657611980310197799; 00358 % $$$ G[m][7][1] = 0.304726500868167195918389047318263; 00359 % $$$ G[m][8][0] = 0.304726500868167195918389047318263; 00360 % $$$ G[m][8][1] = 0.0343243029450971464696306424839376; 00361 % $$$ G[m][9][0] = 0.515842334353591779257463386826430; 00362 % $$$ G[m][9][1] = 0.277716166976391782569581871393723; 00363 % $$$ G[m][10][0] = 0.277716166976391782569581871393723; 00364 % $$$ G[m][10][1] = 0.20644149867001643817295474177985; 00365 % $$$ G[m][11][0] = 0.20644149867001643817295474177985; 00366 % $$$ G[m][11][1] = 0.515842334353591779257463386826430; 00367 % $$$ 00368 % $$$ W[m][0] = 0.5 * 0.053034056314872502857508360921478; 00369 % $$$ W[m][1] = 0.5 * 0.053034056314872502857508360921478; 00370 % $$$ W[m][2] = 0.5 * 0.053034056314872502857508360921478; 00371 % $$$ W[m][3] = 0.5 * 0.087762817428892110073539806278575; 00372 % $$$ W[m][4] = 0.5 * 0.087762817428892110073539806278575; 00373 % $$$ W[m][5] = 0.5 * 0.087762817428892110073539806278575; 00374 % $$$ W[m][6] = 0.5 * 0.057550085569963171476890993800437; 00375 % $$$ W[m][7] = 0.5 * 0.057550085569963171476890993800437; 00376 % $$$ W[m][8] = 0.5 * 0.057550085569963171476890993800437; 00377 % $$$ W[m][9] = 0.5 * 0.13498637401960554892539417233284; 00378 % $$$ W[m][10] = 0.5 * 0.13498637401960554892539417233284; 00379 % $$$ W[m][11] = 0.5 * 0.13498637401960554892539417233284; 00380 % $$$ O[m] = 7; 00381 % $$$ 00382 % $$$ // polynomial degree 8 00383 % $$$ /* 16 inner Gauss points, positive weights */ 00384 % $$$ 00385 % $$$ m=16; 00386 % $$$ G[m][0][0] = 0.33333333333333333333333333333333; 00387 % $$$ G[m][0][1] = 0.33333333333333333333333333333333; 00388 % $$$ G[m][1][0] = 0.17056930775176020662229350149146; 00389 % $$$ G[m][1][1] = 0.17056930775176020662229350149146; 00390 % $$$ G[m][2][0] = 0.17056930775176020662229350149146; 00391 % $$$ G[m][2][1] = 0.65886138449647958675541299701707; 00392 % $$$ G[m][3][0] = 0.65886138449647958675541299701707; 00393 % $$$ G[m][3][1] = 0.17056930775176020662229350149146; 00394 % $$$ G[m][4][0] = 0.050547228317030975458423550596599; 00395 % $$$ G[m][4][1] = 0.050547228317030975458423550596599; 00396 % $$$ G[m][5][0] = 0.050547228317030975458423550596599; 00397 % $$$ G[m][5][1] = 0.89890554336593804908315289880680; 00398 % $$$ G[m][6][0] = 0.89890554336593804908315289880680; 00399 % $$$ G[m][6][1] = 0.050547228317030975458423550596599; 00400 % $$$ G[m][7][0] = 0.45929258829272315602881551449417; 00401 % $$$ G[m][7][1] = 0.45929258829272315602881551449417; 00402 % $$$ G[m][8][0] = 0.45929258829272315602881551449417; 00403 % $$$ G[m][8][1] = 0.08141482341455368794236897101166; 00404 % $$$ G[m][9][0] = 0.08141482341455368794236897101166; 00405 % $$$ G[m][9][1] = 0.45929258829272315602881551449417; 00406 % $$$ G[m][10][0] = 0.72849239295540428124100037917606; 00407 % $$$ G[m][10][1] = 0.26311282963463811342178578628464; 00408 % $$$ G[m][11][0] = 0.72849239295540428124100037917606; 00409 % $$$ G[m][11][1] = 0.00839477740995760533721383453930; 00410 % $$$ G[m][12][0] = 0.26311282963463811342178578628464; 00411 % $$$ G[m][12][1] = 0.72849239295540428124100037917606; 00412 % $$$ G[m][13][0] = 0.26311282963463811342178578628464; 00413 % $$$ G[m][13][1] = 0.00839477740995760533721383453930; 00414 % $$$ G[m][14][0] = 0.00839477740995760533721383453930; 00415 % $$$ G[m][14][1] = 0.72849239295540428124100037917606; 00416 % $$$ G[m][15][0] = 0.00839477740995760533721383453930; 00417 % $$$ G[m][15][1] = 0.26311282963463811342178578628464; 00418 % $$$ 00419 % $$$ W[m][0] = 0.5 * 0.14431560767778716825109111048906; 00420 % $$$ W[m][1] = 0.5 * 0.10321737053471825028179155029213; 00421 % $$$ W[m][2] = 0.5 * 0.10321737053471825028179155029213; 00422 % $$$ W[m][3] = 0.5 * 0.10321737053471825028179155029213; 00423 % $$$ W[m][4] = 0.5 * 0.032458497623198080310925928341780; 00424 % $$$ W[m][5] = 0.5 * 0.032458497623198080310925928341780; 00425 % $$$ W[m][6] = 0.5 * 0.032458497623198080310925928341780; 00426 % $$$ W[m][7] = 0.5 * 0.095091634267284624793896104388584; 00427 % $$$ W[m][8] = 0.5 * 0.095091634267284624793896104388584; 00428 % $$$ W[m][9] = 0.5 * 0.095091634267284624793896104388584; 00429 % $$$ W[m][10] = 0.5 * 0.027230314174434994264844690073909; 00430 % $$$ W[m][11] = 0.5 * 0.027230314174434994264844690073909; 00431 % $$$ W[m][12] = 0.5 * 0.027230314174434994264844690073909; 00432 % $$$ W[m][13] = 0.5 * 0.027230314174434994264844690073909; 00433 % $$$ W[m][14] = 0.5 * 0.027230314174434994264844690073909; 00434 % $$$ W[m][15] = 0.5 * 0.027230314174434994264844690073909; 00435 % $$$ O[m] = 8; 00436 % $$$ 00437 % $$$ // polynomial degree 9 00438 % $$$ /* 19 inner Gauss points, positive weights */ 00439 % $$$ 00440 % $$$ m=19; 00441 % $$$ G[m][0][0] = 0.333333333333333333333333333333333; 00442 % $$$ G[m][0][1] = 0.333333333333333333333333333333333; 00443 % $$$ G[m][1][0] = 0.48968251919873762778370692483619; 00444 % $$$ G[m][1][1] = 0.48968251919873762778370692483619; 00445 % $$$ G[m][2][0] = 0.48968251919873762778370692483619; 00446 % $$$ G[m][2][1] = 0.02063496160252474443258615032762; 00447 % $$$ G[m][3][0] = 0.02063496160252474443258615032762; 00448 % $$$ G[m][3][1] = 0.48968251919873762778370692483619; 00449 % $$$ G[m][4][0] = 0.43708959149293663726993036443535; 00450 % $$$ G[m][4][1] = 0.43708959149293663726993036443535; 00451 % $$$ G[m][5][0] = 0.43708959149293663726993036443535; 00452 % $$$ G[m][5][1] = 0.12582081701412672546013927112929; 00453 % $$$ G[m][6][0] = 0.12582081701412672546013927112929; 00454 % $$$ G[m][6][1] = 0.43708959149293663726993036443535; 00455 % $$$ G[m][7][0] = 0.18820353561903273024096128046733; 00456 % $$$ G[m][7][1] = 0.18820353561903273024096128046733; 00457 % $$$ G[m][8][0] = 0.18820353561903273024096128046733; 00458 % $$$ G[m][8][1] = 0.62359292876193453951807743906533; 00459 % $$$ G[m][9][0] = 0.62359292876193453951807743906533; 00460 % $$$ G[m][9][1] = 0.18820353561903273024096128046733; 00461 % $$$ G[m][10][0] = 0.044729513394452709865106589966276; 00462 % $$$ G[m][10][1] = 0.044729513394452709865106589966276; 00463 % $$$ G[m][11][0] = 0.044729513394452709865106589966276; 00464 % $$$ G[m][11][1] = 0.91054097321109458026978682006745; 00465 % $$$ G[m][12][0] = 0.91054097321109458026978682006745; 00466 % $$$ G[m][12][1] = 0.044729513394452709865106589966276; 00467 % $$$ G[m][13][0] = 0.74119859878449802069007987352342; 00468 % $$$ G[m][13][1] = 0.036838412054736283634817598783385; 00469 % $$$ G[m][14][0] = 0.74119859878449802069007987352342; 00470 % $$$ G[m][14][1] = 0.22196298916076569567510252769319; 00471 % $$$ G[m][15][0] = 0.036838412054736283634817598783385; 00472 % $$$ G[m][15][1] = 0.74119859878449802069007987352342; 00473 % $$$ G[m][16][0] = 0.036838412054736283634817598783385; 00474 % $$$ G[m][16][1] = 0.22196298916076569567510252769319; 00475 % $$$ G[m][17][0] = 0.22196298916076569567510252769319; 00476 % $$$ G[m][17][1] = 0.74119859878449802069007987352342; 00477 % $$$ G[m][18][0] = 0.22196298916076569567510252769319; 00478 % $$$ G[m][18][1] = 0.036838412054736283634817598783385; 00479 % $$$ 00480 % $$$ W[m][0] = 0.5 * 0.097135796282798833819241982507289; 00481 % $$$ W[m][1] = 0.5 * 0.031334700227139070536854831287209; 00482 % $$$ W[m][2] = 0.5 * 0.031334700227139070536854831287209; 00483 % $$$ W[m][3] = 0.5 * 0.031334700227139070536854831287209; 00484 % $$$ W[m][4] = 0.5 * 0.077827541004774279316739356299404; 00485 % $$$ W[m][5] = 0.5 * 0.077827541004774279316739356299404; 00486 % $$$ W[m][6] = 0.5 * 0.077827541004774279316739356299404; 00487 % $$$ W[m][7] = 0.5 * 0.079647738927210253032891774264045; 00488 % $$$ W[m][8] = 0.5 * 0.079647738927210253032891774264045; 00489 % $$$ W[m][9] = 0.5 * 0.079647738927210253032891774264045; 00490 % $$$ W[m][10] = 0.5 * 0.025577675658698031261678798559000; 00491 % $$$ W[m][11] = 0.5 * 0.025577675658698031261678798559000; 00492 % $$$ W[m][12] = 0.5 * 0.025577675658698031261678798559000; 00493 % $$$ W[m][13] = 0.5 * 0.043283539377289377289377289377289; 00494 % $$$ W[m][14] = 0.5 * 0.043283539377289377289377289377289; 00495 % $$$ W[m][15] = 0.5 * 0.043283539377289377289377289377289; 00496 % $$$ W[m][16] = 0.5 * 0.043283539377289377289377289377289; 00497 % $$$ W[m][17] = 0.5 * 0.043283539377289377289377289377289; 00498 % $$$ W[m][18] = 0.5 * 0.043283539377289377289377289377289; 00499 % $$$ O[m] = 9; 00500 % $$$ 00501 % $$$ // polynomial degree 10 00502 % $$$ /* 25 inner Gauss points, positive weights */ 00503 % $$$ m= 25; 00504 % $$$ G[m][0][0] = 0.333333333333333333333333333333333; 00505 % $$$ G[m][0][1] = 0.333333333333333333333333333333333; 00506 % $$$ G[m][1][0] = 0.42508621060209057296952951163804; 00507 % $$$ G[m][1][1] = 0.42508621060209057296952951163804; 00508 % $$$ G[m][2][0] = 0.42508621060209057296952951163804; 00509 % $$$ G[m][2][1] = 0.14982757879581885406094097672391; 00510 % $$$ G[m][3][0] = 0.14982757879581885406094097672391; 00511 % $$$ G[m][3][1] = 0.42508621060209057296952951163804; 00512 % $$$ G[m][4][0] = 0.023308867510000190714466386895980; 00513 % $$$ G[m][4][1] = 0.023308867510000190714466386895980; 00514 % $$$ G[m][5][0] = 0.023308867510000190714466386895980; 00515 % $$$ G[m][5][1] = 0.95338226497999961857106722620804; 00516 % $$$ G[m][6][0] = 0.95338226497999961857106722620804; 00517 % $$$ G[m][6][1] = 0.023308867510000190714466386895980; 00518 % $$$ G[m][7][0] = 0.62830740021349255642083766607883; 00519 % $$$ G[m][7][1] = 0.22376697357697300622568649026820; 00520 % $$$ G[m][8][0] = 0.62830740021349255642083766607883; 00521 % $$$ G[m][8][1] = 0.14792562620953443735347584365296; 00522 % $$$ G[m][9][0] = 0.22376697357697300622568649026820; 00523 % $$$ G[m][9][1] = 0.62830740021349255642083766607883; 00524 % $$$ G[m][10][0] = 0.22376697357697300622568649026820; 00525 % $$$ G[m][10][1] = 0.14792562620953443735347584365296; 00526 % $$$ G[m][11][0] = 0.14792562620953443735347584365296; 00527 % $$$ G[m][11][1] = 0.62830740021349255642083766607883; 00528 % $$$ G[m][12][0] = 0.14792562620953443735347584365296; 00529 % $$$ G[m][12][1] = 0.22376697357697300622568649026820; 00530 % $$$ G[m][13][0] = 0.61131382618139764891875500225390; 00531 % $$$ G[m][13][1] = 0.35874014186443146457815530072385; 00532 % $$$ G[m][14][0] = 0.61131382618139764891875500225390; 00533 % $$$ G[m][14][1] = 0.02994603195417088650308969702225; 00534 % $$$ G[m][15][0] = 0.35874014186443146457815530072385; 00535 % $$$ G[m][15][1] = 0.61131382618139764891875500225390; 00536 % $$$ G[m][16][0] = 0.35874014186443146457815530072385; 00537 % $$$ G[m][16][1] = 0.02994603195417088650308969702225; 00538 % $$$ G[m][17][0] = 0.02994603195417088650308969702225; 00539 % $$$ G[m][17][1] = 0.61131382618139764891875500225390; 00540 % $$$ G[m][18][0] = 0.02994603195417088650308969702225; 00541 % $$$ G[m][18][1] = 0.35874014186443146457815530072385; 00542 % $$$ G[m][19][0] = 0.82107206998562937337354441347218; 00543 % $$$ G[m][19][1] = 0.14329537042686714530585663061732; 00544 % $$$ G[m][20][0] = 0.82107206998562937337354441347218; 00545 % $$$ G[m][20][1] = 0.03563255958750348132059895591050; 00546 % $$$ G[m][21][0] = 0.14329537042686714530585663061732; 00547 % $$$ G[m][21][1] = 0.82107206998562937337354441347218; 00548 % $$$ G[m][22][0] = 0.14329537042686714530585663061732; 00549 % $$$ G[m][22][1] = 0.03563255958750348132059895591050; 00550 % $$$ G[m][23][0] = 0.03563255958750348132059895591050; 00551 % $$$ G[m][23][1] = 0.82107206998562937337354441347218; 00552 % $$$ G[m][24][0] = 0.03563255958750348132059895591050; 00553 % $$$ G[m][24][1] = 0.14329537042686714530585663061732; 00554 % $$$ 00555 % $$$ W[m][0] = 0.5 * 0.079894504741239707831247045213386; 00556 % $$$ W[m][1] = 0.5 * 0.071123802232377334639291287398658; 00557 % $$$ W[m][2] = 0.5 * 0.071123802232377334639291287398658; 00558 % $$$ W[m][3] = 0.5 * 0.071123802232377334639291287398658; 00559 % $$$ W[m][4] = 0.5 * 0.0082238186904641955186466203624719; 00560 % $$$ W[m][5] = 0.5 * 0.0082238186904641955186466203624719; 00561 % $$$ W[m][6] = 0.5 * 0.0082238186904641955186466203624719; 00562 % $$$ W[m][7] = 0.5 * 0.045430592296170018007073629243933; 00563 % $$$ W[m][8] = 0.5 * 0.045430592296170018007073629243933; 00564 % $$$ W[m][9] = 0.5 * 0.045430592296170018007073629243933; 00565 % $$$ W[m][10] = 0.5 * 0.045430592296170018007073629243933; 00566 % $$$ W[m][11] = 0.5 * 0.045430592296170018007073629243933; 00567 % $$$ W[m][12] = 0.5 * 0.045430592296170018007073629243933; 00568 % $$$ W[m][13] = 0.5 * 0.037359856234305276826236499001975; 00569 % $$$ W[m][14] = 0.5 * 0.037359856234305276826236499001975; 00570 % $$$ W[m][15] = 0.5 * 0.037359856234305276826236499001975; 00571 % $$$ W[m][16] = 0.5 * 0.037359856234305276826236499001975; 00572 % $$$ W[m][17] = 0.5 * 0.037359856234305276826236499001975; 00573 % $$$ W[m][18] = 0.5 * 0.037359856234305276826236499001975; 00574 % $$$ W[m][19] = 0.5 * 0.030886656884563988782513077004629; 00575 % $$$ W[m][20] = 0.5 * 0.030886656884563988782513077004629; 00576 % $$$ W[m][21] = 0.5 * 0.030886656884563988782513077004629; 00577 % $$$ W[m][22] = 0.5 * 0.030886656884563988782513077004629; 00578 % $$$ W[m][23] = 0.5 * 0.030886656884563988782513077004629; 00579 % $$$ W[m][24] = 0.5 * 0.030886656884563988782513077004629; 00580 % $$$ O[m] = 10; 00581 % $$$ 00582 % $$$ // polynomial degree 11 00583 % $$$ /* 28 inner Gauss points, positive weights */ 00584 % $$$ 00585 % $$$ m=28; 00586 % $$$ G[m][0][0] = 0.858870281282636704039173938058347; 00587 % $$$ G[m][0][1] = 0.141129718717363295960826061941652; 00588 % $$$ G[m][1][0] = 0.858870281282636704039173938058347; 00589 % $$$ G[m][1][1] = 0.0; 00590 % $$$ G[m][2][0] = 0.141129718717363295960826061941652; 00591 % $$$ G[m][2][1] = 0.858870281282636704039173938058347; 00592 % $$$ G[m][3][0] = 0.141129718717363295960826061941652; 00593 % $$$ G[m][3][1] = 0.0; 00594 % $$$ G[m][4][0] = 0.0; 00595 % $$$ G[m][4][1] = 0.858870281282636704039173938058347; 00596 % $$$ G[m][5][0] = 0.0; 00597 % $$$ G[m][5][1] = 0.141129718717363295960826061941652; 00598 % $$$ G[m][6][0] = 0.333333333333333333333333333333333; 00599 % $$$ G[m][6][1] = 0.333333333333333333333333333333333; 00600 % $$$ G[m][7][0] = 0.025989140928287395260032485498841; 00601 % $$$ G[m][7][1] = 0.025989140928287395260032485498841; 00602 % $$$ G[m][8][0] = 0.025989140928287395260032485498841; 00603 % $$$ G[m][8][1] = 0.94802171814342520947993502900232; 00604 % $$$ G[m][9][0] = 0.94802171814342520947993502900232; 00605 % $$$ G[m][9][1] = 0.025989140928287395260032485498841; 00606 % $$$ G[m][10][0] = 0.094287502647922495630569776275405; 00607 % $$$ G[m][10][1] = 0.094287502647922495630569776275405; 00608 % $$$ G[m][11][0] = 0.094287502647922495630569776275405; 00609 % $$$ G[m][11][1] = 0.81142499470415500873886044744919; 00610 % $$$ G[m][12][0] = 0.81142499470415500873886044744919; 00611 % $$$ G[m][12][1] = 0.094287502647922495630569776275405; 00612 % $$$ G[m][13][0] = 0.49463677501721381374163260230644; 00613 % $$$ G[m][13][1] = 0.49463677501721381374163260230644; 00614 % $$$ G[m][14][0] = 0.49463677501721381374163260230644; 00615 % $$$ G[m][14][1] = 0.01072644996557237251673479538713; 00616 % $$$ G[m][15][0] = 0.01072644996557237251673479538713; 00617 % $$$ G[m][15][1] = 0.49463677501721381374163260230644; 00618 % $$$ G[m][16][0] = 0.20734338261451133345293402411297; 00619 % $$$ G[m][16][1] = 0.20734338261451133345293402411297; 00620 % $$$ G[m][17][0] = 0.20734338261451133345293402411297; 00621 % $$$ G[m][17][1] = 0.58531323477097733309413195177407; 00622 % $$$ G[m][18][0] = 0.58531323477097733309413195177407; 00623 % $$$ G[m][18][1] = 0.20734338261451133345293402411297; 00624 % $$$ G[m][19][0] = 0.43890780570049209506106538163613; 00625 % $$$ G[m][19][1] = 0.43890780570049209506106538163613; 00626 % $$$ G[m][20][0] = 0.43890780570049209506106538163613; 00627 % $$$ G[m][20][1] = 0.12218438859901580987786923672775; 00628 % $$$ G[m][21][0] = 0.12218438859901580987786923672775; 00629 % $$$ G[m][21][1] = 0.43890780570049209506106538163613; 00630 % $$$ G[m][22][0] = 0.67793765488259040154212614118875; 00631 % $$$ G[m][22][1] = 0.044841677589130443309052391468801; 00632 % $$$ G[m][23][0] = 0.67793765488259040154212614118875; 00633 % $$$ G[m][23][1] = 0.27722066752827915514882146734245; 00634 % $$$ G[m][24][0] = 0.044841677589130443309052391468801; 00635 % $$$ G[m][24][1] = 0.67793765488259040154212614118875; 00636 % $$$ G[m][25][0] = 0.044841677589130443309052391468801; 00637 % $$$ G[m][25][1] = 0.27722066752827915514882146734245; 00638 % $$$ G[m][26][0] = 0.27722066752827915514882146734245; 00639 % $$$ G[m][26][1] = 0.67793765488259040154212614118875; 00640 % $$$ G[m][27][0] = 0.27722066752827915514882146734245; 00641 % $$$ G[m][27][1] = 0.044841677589130443309052391468801; 00642 % $$$ 00643 % $$$ W[m][0] = 0.5 * 0.0073623837833005542642588950473806; 00644 % $$$ W[m][1] = 0.5 * 0.0073623837833005542642588950473806; 00645 % $$$ W[m][2] = 0.5 * 0.0073623837833005542642588950473806; 00646 % $$$ 00647 % $$$ W[m][3] = 0.5 * 0.0073623837833005542642588950473806; 00648 % $$$ W[m][4] = 0.5 * 0.0073623837833005542642588950473806; 00649 % $$$ W[m][5] = 0.5 * 0.0073623837833005542642588950473806; 00650 % $$$ 00651 % $$$ W[m][6] = 0.5 * 0.087977301162232238798093169321456; 00652 % $$$ W[m][7] = 0.5 * 0.0087443115537360230495164287998252; 00653 % $$$ W[m][8] = 0.5 * 0.0087443115537360230495164287998252; 00654 % $$$ W[m][9] = 0.5 * 0.0087443115537360230495164287998252; 00655 % $$$ 00656 % $$$ W[m][10] = 0.5 * 0.038081571993934937515024339435614; 00657 % $$$ W[m][11] = 0.5 * 0.038081571993934937515024339435614; 00658 % $$$ W[m][12] = 0.5 * 0.038081571993934937515024339435614; 00659 % $$$ 00660 % $$$ W[m][13] = 0.5 * 0.018855448056131292058476782591115; 00661 % $$$ W[m][14] = 0.5 * 0.018855448056131292058476782591115; 00662 % $$$ W[m][15] = 0.5 * 0.018855448056131292058476782591115; 00663 % $$$ 00664 % $$$ W[m][16] = 0.5 * 0.072159697544739526124029988586463; 00665 % $$$ W[m][17] = 0.5 * 0.072159697544739526124029988586463; 00666 % $$$ W[m][18] = 0.5 * 0.072159697544739526124029988586463; 00667 % $$$ 00668 % $$$ W[m][19] = 0.5 * 0.069329138705535899841765650903814; 00669 % $$$ W[m][20] = 0.5 * 0.069329138705535899841765650903814; 00670 % $$$ W[m][21] = 0.5 * 0.069329138705535899841765650903814; 00671 % $$$ 00672 % $$$ W[m][22] = 0.5 * 0.041056315429288566641652314907294; 00673 % $$$ W[m][23] = 0.5 * 0.041056315429288566641652314907294; 00674 % $$$ W[m][24] = 0.5 * 0.041056315429288566641652314907294; 00675 % $$$ 00676 % $$$ W[m][25] = 0.5 * 0.041056315429288566641652314907294; 00677 % $$$ W[m][26] = 0.5 * 0.041056315429288566641652314907294; 00678 % $$$ W[m][27] = 0.5 * 0.041056315429288566641652314907294; 00679 % $$$ O[m] = 11; 00680 % $$$ 00681 % $$$ // polynomial degree 12 00682 % $$$ /* 33 inner Gauss points, positive weights */ 00683 % $$$ m=33; 00684 % $$$ G[m][0][0] = 0.02356522045239; 00685 % $$$ G[m][0][1] = 0.488217389773805; 00686 % $$$ G[m][1][0] = 0.488217389773805; 00687 % $$$ G[m][1][1] = 0.02356522045239; 00688 % $$$ G[m][2][0] = 0.488217389773805; 00689 % $$$ G[m][2][1] = 0.488217389773805; 00690 % $$$ G[m][3][0] = 0.43972439229446; 00691 % $$$ G[m][3][1] = 0.43972439229446; 00692 % $$$ G[m][4][0] = 0.43972439229446; 00693 % $$$ G[m][4][1] = 0.120551215411079; 00694 % $$$ G[m][5][0] = 0.120551215411079; 00695 % $$$ G[m][5][1] = 0.43972439229446; 00696 % $$$ G[m][6][0] = 0.271210385012116; 00697 % $$$ G[m][6][1] = 0.271210385012116; 00698 % $$$ G[m][7][0] = 0.271210385012116; 00699 % $$$ G[m][7][1] = 0.457579229975768; 00700 % $$$ G[m][8][0] = 0.457579229975768; 00701 % $$$ G[m][8][1] = 0.271210385012116; 00702 % $$$ G[m][9][0] = 0.127576145541586; 00703 % $$$ G[m][9][1] = 0.127576145541586; 00704 % $$$ G[m][10][0] = 0.127576145541586; 00705 % $$$ G[m][10][1] = 0.7448477089168279; 00706 % $$$ G[m][11][0] = 0.7448477089168279; 00707 % $$$ G[m][11][1] = 0.127576145541586; 00708 % $$$ G[m][12][0] = 0.02131735045321; 00709 % $$$ G[m][12][1] = 0.02131735045321; 00710 % $$$ G[m][13][0] = 0.02131735045321; 00711 % $$$ G[m][13][1] = 0.9573652990935799; 00712 % $$$ G[m][14][0] = 0.9573652990935799; 00713 % $$$ G[m][14][1] = 0.02131735045321; 00714 % $$$ G[m][15][0] = 0.115343494534698; 00715 % $$$ G[m][15][1] = 0.275713269685514; 00716 % $$$ G[m][16][0] = 0.115343494534698; 00717 % $$$ G[m][16][1] = 0.6089432357797879; 00718 % $$$ G[m][17][0] = 0.275713269685514; 00719 % $$$ G[m][17][1] = 0.115343494534698; 00720 % $$$ G[m][18][0] = 0.275713269685514; 00721 % $$$ G[m][18][1] = 0.6089432357797879; 00722 % $$$ G[m][19][0] = 0.6089432357797879; 00723 % $$$ G[m][19][1] = 0.115343494534698; 00724 % $$$ G[m][20][0] = 0.6089432357797879; 00725 % $$$ G[m][20][1] = 0.275713269685514; 00726 % $$$ G[m][21][0] = 0.022838332222257; 00727 % $$$ G[m][21][1] = 0.28132558098994; 00728 % $$$ G[m][22][0] = 0.022838332222257; 00729 % $$$ G[m][22][1] = 0.6958360867878031; 00730 % $$$ G[m][23][0] = 0.28132558098994; 00731 % $$$ G[m][23][1] = 0.022838332222257; 00732 % $$$ G[m][24][0] = 0.28132558098994; 00733 % $$$ G[m][24][1] = 0.6958360867878031; 00734 % $$$ G[m][25][0] = 0.6958360867878031; 00735 % $$$ G[m][25][1] = 0.022838332222257; 00736 % $$$ G[m][26][0] = 0.6958360867878031; 00737 % $$$ G[m][26][1] = 0.28132558098994; 00738 % $$$ G[m][27][0] = 0.02573405054833; 00739 % $$$ G[m][27][1] = 0.116251915907597; 00740 % $$$ G[m][28][0] = 0.02573405054833; 00741 % $$$ G[m][28][1] = 0.858014033544073; 00742 % $$$ G[m][29][0] = 0.116251915907597; 00743 % $$$ G[m][29][1] = 0.02573405054833; 00744 % $$$ G[m][30][0] = 0.116251915907597; 00745 % $$$ G[m][30][1] = 0.858014033544073; 00746 % $$$ G[m][31][0]= 0.858014033544073; 00747 % $$$ G[m][31][1] =0.02573405054833; 00748 % $$$ G[m][32][0] =0.858014033544073; 00749 % $$$ G[m][32][1]= 0.116251915907597; 00750 % $$$ 00751 % $$$ W[m][0] = 0.5 * 0.025731066440455; 00752 % $$$ W[m][1] = 0.5 * 0.025731066440455; 00753 % $$$ W[m][2] = 0.5 * 0.025731066440455; 00754 % $$$ W[m][3] = 0.5 * 0.043692544538038; 00755 % $$$ W[m][4] = 0.5 * 0.043692544538038; 00756 % $$$ W[m][5] = 0.5 * 0.043692544538038; 00757 % $$$ W[m][6] = 0.5 * 0.062858224217885; 00758 % $$$ W[m][7] = 0.5 * 0.062858224217885; 00759 % $$$ W[m][8] = 0.5 * 0.062858224217885; 00760 % $$$ W[m][9] = 0.5 * 0.034796112930709; 00761 % $$$ W[m][10] = 0.5 * 0.034796112930709; 00762 % $$$ W[m][11] = 0.5 * 0.034796112930709; 00763 % $$$ W[m][12] = 0.5 * 0.006166261051559; 00764 % $$$ W[m][13] = 0.5 * 0.006166261051559; 00765 % $$$ W[m][14] = 0.5 * 0.006166261051559; 00766 % $$$ W[m][15] = 0.5 * 0.040371557766381; 00767 % $$$ W[m][16] = 0.5 * 0.040371557766381; 00768 % $$$ W[m][17] = 0.5 * 0.040371557766381; 00769 % $$$ W[m][18] = 0.5 * 0.040371557766381; 00770 % $$$ W[m][19] = 0.5 * 0.040371557766381; 00771 % $$$ W[m][20] = 0.5 * 0.040371557766381; 00772 % $$$ W[m][21] = 0.5 * 0.022356773202303; 00773 % $$$ W[m][22] = 0.5 * 0.022356773202303; 00774 % $$$ W[m][23] = 0.5 * 0.022356773202303; 00775 % $$$ W[m][24] = 0.5 * 0.022356773202303; 00776 % $$$ W[m][25] = 0.5 * 0.022356773202303; 00777 % $$$ W[m][26] = 0.5 * 0.022356773202303; 00778 % $$$ W[m][27] = 0.5 * 0.017316231108659; 00779 % $$$ W[m][28] = 0.5 * 0.017316231108659; 00780 % $$$ W[m][29] = 0.5 * 0.017316231108659; 00781 % $$$ W[m][30] = 0.5 * 0.017316231108659; 00782 % $$$ W[m][31] = 0.5 * 0.017316231108659; 00783 % $$$ W[m][32] = 0.5 * 0.017316231108659; 00784 % $$$ O[m] = 12; 00785 % $$$ } 00786 % $$$ 00787 % $$$ FieldVector<double, 2> point(int m, int i) 00788 % $$$ { 00789 % $$$ return G[m][i]; 00790 % $$$ } 00791 % $$$ 00792 % $$$ double weight (int m, int i) 00793 % $$$ { 00794 % $$$ return W[m][i]; 00795 % $$$ } 00796 % $$$ 00797 % $$$ int order (int m) 00798 % $$$ { 00799 % $$$ return O[m]; 00800 % $$$ } 00801 % $$$ 00802 % $$$ private: 00803 % $$$ FieldVector<double, 2> G[MAXP+1][MAXP]; 00804 % $$$ 00805 % $$$ double W[MAXP+1][MAXP]; // weights associated with points 00806 % $$$ int O[MAXP+1]; // order of the rule 00807 % $$$ }; 00808 % $$$ 00809 % $$$ template<typename ct> 00810 % $$$ SimplexQuadratureRule<ct,2>.SimplexQuadratureRule(int p) : QuadratureRule<ct,2>(GeometryType(GeometryType.simplex, 2)) 00811 % $$$ { 00812 % $$$ int m; 00813 % $$$ if (p>highest_order) 00814 % $$$ DUNE_THROW(QuadratureOrderOutOfRange, 00815 % $$$ "QuadratureRule for order " << p << " and GeometryType " 00816 % $$$ << this->type() << " not available"); 00817 % $$$ 00818 % $$$ if (p>SimplexQuadraturePoints<2>.highest_order) 00819 % $$$ { 00820 % $$$ // Define the quadrature rules... 00821 % $$$ QuadratureRule<ct,1> gauss1D = 00822 % $$$ QuadratureRules<ct,1>.rule(GeometryType.cube, p, QuadratureType.Gauss); 00823 % $$$ QuadratureRule<ct,1> jac1D = 00824 % $$$ QuadratureRules<ct,1>.rule(GeometryType.cube, p, QuadratureType.Jacobian_1_0); 00825 % $$$ 00826 % $$$ // Compute the conical product 00827 % $$$ for (typename QuadratureRule<ct,1>.const_iterator 00828 % $$$ gp=gauss1D.begin(); gp!=gauss1D.end(); ++gp) 00829 % $$$ { 00830 % $$$ for (typename QuadratureRule<ct,1>.const_iterator 00831 % $$$ jp=jac1D.begin(); jp!=jac1D.end(); ++jp) 00832 % $$$ { 00833 % $$$ // compute coordinates and weight 00834 % $$$ double weight = 1.0; 00835 % $$$ FieldVector<ct,2> local(0.0); 00836 % $$$ local[0] = jp->position()[0]; 00837 % $$$ local[1] = gp->position()[0] * (1. - jp->position()[0]); 00838 % $$$ weight = gp->weight() * jp->weight() * (1. - jp->position()[0]); 00839 % $$$ // put in container 00840 % $$$ push_back(QuadraturePoint<ct,d>(local,weight)); 00841 % $$$ } 00842 % $$$ } 00843 % $$$ this->delivered_order = std.min(gauss1D.order(), jac1D.order()); 00844 % $$$ return; 00845 % $$$ } 00846 % $$$ 00847 % $$$ switch(p) 00848 % $$$ { 00849 % $$$ case 0: // to be verified 00850 % $$$ m=1; // to be verified 00851 % $$$ break; 00852 % $$$ case 1: 00853 % $$$ m=1; 00854 % $$$ break; 00855 % $$$ case 2: 00856 % $$$ m=3; 00857 % $$$ break; 00858 % $$$ case 3: 00859 % $$$ m=4; 00860 % $$$ break; 00861 % $$$ case 4: 00862 % $$$ m=6; 00863 % $$$ break; 00864 % $$$ case 5: 00865 % $$$ m=7; 00866 % $$$ break; 00867 % $$$ case 6: 00868 % $$$ m=12; 00869 % $$$ break; 00870 % $$$ case 7: 00871 % $$$ m=12; 00872 % $$$ break; 00873 % $$$ case 8: 00874 % $$$ m=16; 00875 % $$$ break; 00876 % $$$ case 9: 00877 % $$$ m=19; 00878 % $$$ break; 00879 % $$$ case 10: 00880 % $$$ m=25; 00881 % $$$ break; 00882 % $$$ case 11: 00883 % $$$ m=28; 00884 % $$$ break; 00885 % $$$ case 12: 00886 % $$$ m=33; 00887 % $$$ break; 00888 % $$$ default: m=33; 00889 % $$$ } 00890 % $$$ this->delivered_order = SimplexQuadraturePointsSingleton<2>.sqp.order(m); 00891 % $$$ FieldVector<ct, d> local; 00892 % $$$ double weight; 00893 % $$$ for(int i=0;i<m;++i) 00894 % $$$ { 00895 % $$$ for(int k=0;k<d;++k) 00896 % $$$ local[k]=SimplexQuadraturePointsSingleton<2>.sqp.point(m,i)[k]; 00897 % $$$ weight=SimplexQuadraturePointsSingleton<2>.sqp.weight(m,i); 00898 % $$$ // put in container 00899 % $$$ push_back(QuadraturePoint<ct,d>(local,weight)); 00900 % $$$ } 00901 % $$$ } 00902 % $$$ 00903 % $$$ 00904 % $$$ template<> 00905 % $$$ class SimplexQuadraturePoints<3> 00906 % $$$ { 00907 % $$$ public: 00908 % $$$ enum { MAXP=15}; 00909 % $$$ enum { highest_order=5 }; 00910 % $$$ 00911 % $$$ 00912 % $$$ SimplexQuadraturePoints() 00913 % $$$ { 00914 % $$$ int m = 0; 00915 % $$$ O[m] = 0; 00916 % $$$ 00917 % $$$ // polynom degree 1 00918 % $$$ m = 1; 00919 % $$$ G[m][0][0] = 0.25; 00920 % $$$ G[m][0][1] = 0.25; 00921 % $$$ G[m][0][2] = 0.25; 00922 % $$$ 00923 % $$$ W[m][0] = 1.0/6.0; 00924 % $$$ O[m] = 1; 00925 % $$$ 00926 % $$$ // polynom degree 2 00927 % $$$ // symmetric 00928 % $$$ m = 4; 00929 % $$$ static const double m_4_a = 0.585410196624968500; 00930 % $$$ static const double m_4_b = 0.138196601125010500; 00931 % $$$ G[m][0] = m_4_b; 00932 % $$$ G[m][1] = m_4_b; 00933 % $$$ G[m][2] = m_4_b; 00934 % $$$ G[m][3] = m_4_b; 00935 % $$$ G[m][0][0] = m_4_a; 00936 % $$$ G[m][1][1] = m_4_a; 00937 % $$$ G[m][2][2] = m_4_a; 00938 % $$$ 00939 % $$$ W[m][0] = 1.0/4.0/6.0; 00940 % $$$ W[m][1] = 1.0/4.0/6.0; 00941 % $$$ W[m][2] = 1.0/4.0/6.0; 00942 % $$$ W[m][3] = 1.0/4.0/6.0; 00943 % $$$ O[m] = 2; 00944 % $$$ 00945 % $$$ // polynom degree 3 00946 % $$$ // symmetric 00947 % $$$ m = 8; 00948 % $$$ G[m][0][0] = 0.0; 00949 % $$$ G[m][0][1] = 0.0; 00950 % $$$ G[m][0][2] = 0.0; 00951 % $$$ G[m][1][0] = 1.0; 00952 % $$$ G[m][1][1] = 0.0; 00953 % $$$ G[m][1][2] = 0.0; 00954 % $$$ G[m][2][0] = 0.0; 00955 % $$$ G[m][2][1] = 1.0; 00956 % $$$ G[m][2][2] = 0.0; 00957 % $$$ G[m][3][0] = 0.0; 00958 % $$$ G[m][3][1] = 0.0; 00959 % $$$ G[m][3][2] = 1.0; 00960 % $$$ G[m][4][0] = 1.0/3.0; 00961 % $$$ G[m][4][1] = 1.0/3.0; 00962 % $$$ G[m][4][2] = 0.0; 00963 % $$$ G[m][5][0] = 1.0/3.0; 00964 % $$$ G[m][5][1] = 0.0; 00965 % $$$ G[m][5][2] = 1.0/3.0; 00966 % $$$ G[m][6][0] = 0.0; 00967 % $$$ G[m][6][1] = 1.0/3.0; 00968 % $$$ G[m][6][2] = 1.0/3.0; 00969 % $$$ G[m][7][0] = 1.0/3.0; 00970 % $$$ G[m][7][1] = 1.0/3.0; 00971 % $$$ G[m][7][2] = 1.0/3.0; 00972 % $$$ W[m][0] = 0.025/6.0; 00973 % $$$ W[m][1] = 0.025/6.0; 00974 % $$$ W[m][2] = 0.025/6.0; 00975 % $$$ W[m][3] = 0.025/6.0; 00976 % $$$ W[m][4] = 0.225/6.0; 00977 % $$$ W[m][5] = 0.225/6.0; 00978 % $$$ W[m][6] = 0.225/6.0; 00979 % $$$ W[m][7] = 0.225/6.0; 00980 % $$$ O[m] = 3; 00981 % $$$ 00982 % $$$ 00983 % $$$ // polynomial degree 5 00984 % $$$ // symmetric points 00985 % $$$ static const double s_1=0.09197107805272303279; /* (7 - sqrt(15) ) / 34 */ 00986 % $$$ static const double s_2=0.31979362782962990839; /* (7 + sqrt(15) ) / 34 */ 00987 % $$$ static const double t_1=0.72408676584183090164; /* (13 + 3*sqrt(15) ) / 34 */ 00988 % $$$ static const double t_2=0.04061911651111027484; /* (13 - 3*sqrt(15) ) / 34 */ 00989 % $$$ static const double u =0.05635083268962915574; /* (10 - 2*sqrt(15) ) / 40 */ 00990 % $$$ static const double v =0.44364916731037084426; /* (10 + 2*sqrt(15) ) / 40 */ 00991 % $$$ static const double A =0.019753086419753086420; /* 16 / 135 / vol */ 00992 % $$$ static const double B_1=0.011989513963169770001; /* (2665 + 14*sqrt(15) ) / 37800 / vol */ 00993 % $$$ static const double B_2=0.011511367871045397547; /* (2665 - 14*sqrt(15) ) / 37800 / vol */ 00994 % $$$ static const double C =0.0088183421516754850088; /* 20 / 378 / vol */ 00995 % $$$ 00996 % $$$ m=15; 00997 % $$$ G[m][0][0] = 0.25; 00998 % $$$ G[m][0][1] = 0.25; 00999 % $$$ G[m][0][2] = 0.25; 01000 % $$$ G[m][1][0] = s_1; 01001 % $$$ G[m][1][1] = s_1; 01002 % $$$ G[m][1][2] = s_1; 01003 % $$$ G[m][2][0] = t_1; 01004 % $$$ G[m][2][1] = s_1; 01005 % $$$ G[m][2][2] = s_1; 01006 % $$$ G[m][3][0] = s_1; 01007 % $$$ G[m][3][1] = t_1; 01008 % $$$ G[m][3][2] = s_1; 01009 % $$$ G[m][4][0] = s_1; 01010 % $$$ G[m][4][1] = s_1; 01011 % $$$ G[m][4][2] = t_1; 01012 % $$$ G[m][5][0] = s_2; 01013 % $$$ G[m][5][1] = s_2; 01014 % $$$ G[m][5][2] = s_2; 01015 % $$$ G[m][6][0] = t_2; 01016 % $$$ G[m][6][1] = s_2; 01017 % $$$ G[m][6][2] = s_2; 01018 % $$$ G[m][7][0] = s_2; 01019 % $$$ G[m][7][1] = t_2; 01020 % $$$ G[m][7][2] = s_2; 01021 % $$$ G[m][8][0] = s_2; 01022 % $$$ G[m][8][1] = s_2; 01023 % $$$ G[m][8][2] = t_2; 01024 % $$$ G[m][9][0] = v; 01025 % $$$ G[m][9][1] = u; 01026 % $$$ G[m][9][2] = u; 01027 % $$$ G[m][10][0] = u; 01028 % $$$ G[m][10][1] = v; 01029 % $$$ G[m][10][2] = u; 01030 % $$$ G[m][11][0] = u; 01031 % $$$ G[m][11][1] = u; 01032 % $$$ G[m][11][2] = v; 01033 % $$$ G[m][12][0] = v; 01034 % $$$ G[m][12][1] = v; 01035 % $$$ G[m][12][2] = u; 01036 % $$$ G[m][13][0] = v; 01037 % $$$ G[m][13][1] = u; 01038 % $$$ G[m][13][2] = v; 01039 % $$$ G[m][14][0] = u; 01040 % $$$ G[m][14][1] = v; 01041 % $$$ G[m][14][2] = v; 01042 % $$$ 01043 % $$$ W[m][0] = A; 01044 % $$$ W[m][1] = B_1; 01045 % $$$ W[m][2] = B_1; 01046 % $$$ W[m][3] = B_1; 01047 % $$$ W[m][4] = B_1; 01048 % $$$ W[m][5] = B_2; 01049 % $$$ W[m][6] = B_2; 01050 % $$$ W[m][7] = B_2; 01051 % $$$ W[m][8] = B_2; 01052 % $$$ W[m][9] = C; 01053 % $$$ W[m][10] = C; 01054 % $$$ W[m][11] = C; 01055 % $$$ W[m][12] = C; 01056 % $$$ W[m][13] = C; 01057 % $$$ W[m][14] = C; 01058 % $$$ O[m] = 5; 01059 % $$$ 01060 % $$$ } 01061 % $$$ 01062 % $$$ FieldVector<double, 3> point(int m, int i) 01063 % $$$ { 01064 % $$$ return G[m][i]; 01065 % $$$ } 01066 % $$$ 01067 % $$$ double weight (int m, int i) 01068 % $$$ { 01069 % $$$ return W[m][i]; 01070 % $$$ } 01071 % $$$ 01072 % $$$ int order (int m) 01073 % $$$ { 01074 % $$$ return O[m]; 01075 % $$$ } 01076 % $$$ 01077 % $$$ private: 01078 % $$$ FieldVector<double, 3> G[MAXP+1][MAXP]; 01079 % $$$ double W[MAXP+1][MAXP]; // weights associated with points 01080 % $$$ int O[MAXP+1]; // order of the rule 01081 % $$$ }; 01082 % $$$ 01083 % $$$ template<typename ct> 01084 % $$$ SimplexQuadratureRule<ct,3>.SimplexQuadratureRule(int p) : QuadratureRule<ct,3>(GeometryType(GeometryType.simplex, 3)) 01085 % $$$ { 01086 % $$$ int m; 01087 % $$$ if (p>highest_order) 01088 % $$$ DUNE_THROW(QuadratureOrderOutOfRange, 01089 % $$$ "QuadratureRule for order " << p << " and GeometryType " 01090 % $$$ << this->type() << " not available"); 01091 % $$$ 01092 % $$$ // quadrature rules >= 3 are known to work. 01093 % $$$ if (p>SimplexQuadraturePoints<3>.highest_order) 01094 % $$$ { 01095 % $$$ // Define the quadrature rules... 01096 % $$$ QuadratureRule<ct,1> gauss1D = 01097 % $$$ QuadratureRules<ct,1>.rule(GeometryType.cube, p, QuadratureType.Gauss); 01098 % $$$ QuadratureRule<ct,1> jacA1D = 01099 % $$$ QuadratureRules<ct,1>.rule(GeometryType.cube, p, QuadratureType.Jacobian_1_0); 01100 % $$$ QuadratureRule<ct,1> jacB1D = 01101 % $$$ QuadratureRules<ct,1>.rule(GeometryType.cube, p, QuadratureType.Jacobian_2_0); 01102 % $$$ 01103 % $$$ // Compute the tensor product 01104 % $$$ 01105 % $$$ // All rules should be of the same order 01106 % $$$ assert(gauss1D.size() == jacA1D.size()); 01107 % $$$ assert(gauss1D.size() == jacB1D.size()); 01108 % $$$ 01109 % $$$ // Compute the conical product 01110 % $$$ for (typename QuadratureRule<ct,1>.const_iterator 01111 % $$$ gp=gauss1D.begin(); gp!=gauss1D.end(); ++gp) 01112 % $$$ { 01113 % $$$ for (typename QuadratureRule<ct,1>.const_iterator 01114 % $$$ j1p=jacA1D.begin(); j1p!=jacA1D.end(); ++j1p) 01115 % $$$ { 01116 % $$$ for (typename QuadratureRule<ct,1>.const_iterator 01117 % $$$ j2p=jacB1D.begin(); j2p!=jacB1D.end(); ++j2p) 01118 % $$$ { 01119 % $$$ // compute coordinates and weight 01120 % $$$ double weight = 1.0; 01121 % $$$ FieldVector<ct, d> local; 01122 % $$$ local[0] = j2p->position()[0]; 01123 % $$$ local[1] = j1p->position()[0] * (1.0-j2p->position()[0]); 01124 % $$$ local[2] = gp->position()[0] * (1.-j1p->position()[0]) * (1.0-j2p->position()[0]); 01125 % $$$ weight = (1.0-j2p->position()[0]) * (1.-j1p->position()[0]) * (1.0-j2p->position()[0]) 01126 % $$$ * gp->weight() * j1p->weight() * j2p->weight(); 01127 % $$$ // put in container 01128 % $$$ push_back(QuadraturePoint<ct,d>(local,weight)); 01129 % $$$ } 01130 % $$$ } 01131 % $$$ } 01132 % $$$ this->delivered_order = std.min(gauss1D.order(), std.min(jacA1D.order(), jacB1D.order())); 01133 % $$$ return; 01134 % $$$ } 01135 % $$$ 01136 % $$$ switch(p) 01137 % $$$ { 01138 % $$$ case 0: // to be verified 01139 % $$$ m=1; // to be verified 01140 % $$$ break; 01141 % $$$ case 1: 01142 % $$$ m=1; 01143 % $$$ break; 01144 % $$$ case 2: 01145 % $$$ m=4; 01146 % $$$ break; 01147 % $$$ case 3: 01148 % $$$ m=8; 01149 % $$$ break; 01150 % $$$ case 4: 01151 % $$$ case 5: 01152 % $$$ m=15; 01153 % $$$ break; 01154 % $$$ default: m=15; 01155 % $$$ } 01156 % $$$ this->delivered_order = SimplexQuadraturePointsSingleton<3>.sqp.order(m); 01157 % $$$ FieldVector<ct, d> local; 01158 % $$$ double weight; 01159 % $$$ for(int i=0;i<m;++i) 01160 % $$$ { 01161 % $$$ for(int k=0;k<d;++k) 01162 % $$$ local[k]=SimplexQuadraturePointsSingleton<3>.sqp.point(m,i)[k]; 01163 % $$$ weight=SimplexQuadraturePointsSingleton<3>.sqp.weight(m,i); 01164 % $$$ // put in container 01165 % $$$ push_back(QuadraturePoint<ct,d>(local,weight)); 01166 % $$$ 01167 % $$$ } 01168 % $$$ 01169 % $$$ } 01170 % $$$ 01171 % $$$ 01172 % $$$ SimplexQuadraturePoints<2> SimplexQuadraturePointsSingleton<2>.sqp; 01173 % $$$ 01174 % $$$ 01175 % $$$ SimplexQuadraturePoints<3> SimplexQuadraturePointsSingleton<3>.sqp; 01176 % $$$ 01177 % $$$ 01178 % $$$ PrismQuadraturePoints<3> PrismQuadraturePointsSingleton<3>.prqp; 01179 % $$$ 01180 % $$$ 01181 % $$$ PyramidQuadraturePoints PyramidQuadraturePointsSingleton<3>.pyqp; 01182 % $$$ 01183 % $$$ template SimplexQuadratureRule<float, 2>.SimplexQuadratureRule(int); 01184 % $$$ template SimplexQuadratureRule<double, 2>.SimplexQuadratureRule(int); 01185 % $$$ template SimplexQuadratureRule<float, 3>.SimplexQuadratureRule(int); 01186 % $$$ template SimplexQuadratureRule<double, 3>.SimplexQuadratureRule(int); 01187 % $$$ 01188 % $$$ } // namespace 01189 %| \docupdate
