00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include <cmath>
00030
00031 #include "GaussianQuadratureRule.hpp"
00032 #include "Exception.hpp"
00033 #include "UblasCustomFunctions.hpp"
00034
00035 template<unsigned ELEMENT_DIM>
00036 const ChastePoint<ELEMENT_DIM>& GaussianQuadratureRule<ELEMENT_DIM>::rGetQuadPoint(unsigned index) const
00037 {
00038 assert(index < mNumQuadPoints);
00039 return mPoints[index];
00040 }
00041
00042 template<unsigned ELEMENT_DIM>
00043 double GaussianQuadratureRule<ELEMENT_DIM>::GetWeight(unsigned index) const
00044 {
00045 assert(index < mNumQuadPoints);
00046 return mWeights[index];
00047 }
00048
00049 template<unsigned ELEMENT_DIM>
00050 unsigned GaussianQuadratureRule<ELEMENT_DIM>::GetNumQuadPoints() const
00051 {
00052 return mNumQuadPoints;
00053 }
00054
00059 template<>
00060 GaussianQuadratureRule<0>::GaussianQuadratureRule(unsigned numPointsInEachDimension)
00061 {
00062 mNumQuadPoints = 1;
00063 mWeights.reserve(mNumQuadPoints);
00064 mPoints.reserve(mNumQuadPoints);
00065 mWeights.push_back(1);
00066 mPoints.push_back(ChastePoint<0>());
00067 }
00068
00073 template<>
00074 GaussianQuadratureRule<1>::GaussianQuadratureRule(unsigned numPointsInEachDimension)
00075 {
00076 mNumQuadPoints = numPointsInEachDimension;
00077
00078 mWeights.reserve(mNumQuadPoints);
00079 mPoints.reserve(mNumQuadPoints);
00080 switch (numPointsInEachDimension)
00081 {
00082 case 1:
00083 mWeights.push_back(1);
00084 mPoints.push_back(ChastePoint<1>(0.5));
00085 break;
00086
00087 case 2:
00088 mWeights.push_back(0.5);
00089 mWeights.push_back(0.5);
00090
00091 mPoints.push_back(ChastePoint<1>(0.21132486540519));
00092 mPoints.push_back(ChastePoint<1>(0.78867513459481));
00093 break;
00094
00095 case 3:
00096 mWeights.push_back(5.0/18.0);
00097 mWeights.push_back(4.0/9.0);
00098 mWeights.push_back(5.0/18.0);
00099
00100 mPoints.push_back(ChastePoint<1>(0.1127016654));
00101 mPoints.push_back(ChastePoint<1>(0.5));
00102 mPoints.push_back(ChastePoint<1>(0.8872983346));
00103 break;
00104
00105 default:
00106 EXCEPTION("Number of gauss points per dimension not supported.");
00107 }
00108 }
00109
00114 template<>
00115 GaussianQuadratureRule<2>::GaussianQuadratureRule(unsigned numPointsInEachDimension)
00116 {
00117 mNumQuadPoints = numPointsInEachDimension * numPointsInEachDimension;
00118
00119 mWeights.reserve(mNumQuadPoints);
00120 mPoints.reserve(mNumQuadPoints);
00121
00122 switch (numPointsInEachDimension)
00123 {
00124 case 1:
00125 mWeights.push_back(0.5);
00126 mPoints.push_back(ChastePoint<2>(0.25,0.5));
00127 break;
00128
00129 case 2:
00130 mWeights.push_back(0.19716878364870);
00131 mWeights.push_back(0.19716878364870);
00132 mWeights.push_back(0.05283121635130);
00133 mWeights.push_back(0.05283121635130);
00134
00135 mPoints.push_back(ChastePoint<2>(0.16666666666667,0.21132486540519));
00136 mPoints.push_back(ChastePoint<2>(0.62200846792815,0.21132486540519));
00137 mPoints.push_back(ChastePoint<2>(0.04465819873852,0.78867513459481));
00138 mPoints.push_back(ChastePoint<2>(0.16666666666667,0.78867513459481));
00139 break;
00140
00141 case 3:
00142 mWeights.push_back(0.06846437766975);
00143 mWeights.push_back(0.10954300427160);
00144 mWeights.push_back(0.06846437766975);
00145 mWeights.push_back(0.06172839506173);
00146 mWeights.push_back(0.09876543209877);
00147 mWeights.push_back(0.06172839506173);
00148 mWeights.push_back(0.00869611615741);
00149 mWeights.push_back(0.01391378585185);
00150 mWeights.push_back(0.00869611615741);
00151
00152 mPoints.push_back(ChastePoint<2>(0.10000000001607,0.11270166540000));
00153 mPoints.push_back(ChastePoint<2>(0.44364916730000,0.11270166540000));
00154 mPoints.push_back(ChastePoint<2>(0.78729833458393,0.11270166540000));
00155 mPoints.push_back(ChastePoint<2>(0.05635083270000,0.50000000000000));
00156 mPoints.push_back(ChastePoint<2>(0.25000000000000,0.50000000000000));
00157 mPoints.push_back(ChastePoint<2>(0.44364916730000,0.50000000000000));
00158 mPoints.push_back(ChastePoint<2>(0.01270166538393,0.88729833460000));
00159 mPoints.push_back(ChastePoint<2>(0.05635083270000,0.88729833460000));
00160 mPoints.push_back(ChastePoint<2>(0.10000000001607,0.88729833460000));
00161 break;
00162
00163 default:
00164 EXCEPTION("Number of gauss points per dimension not supported.");
00165 }
00166 }
00167
00172 template<>
00173 GaussianQuadratureRule<3>::GaussianQuadratureRule(unsigned numPointsInEachDimension)
00174 {
00175 mNumQuadPoints = numPointsInEachDimension * numPointsInEachDimension * numPointsInEachDimension;
00176
00177 mWeights.reserve(mNumQuadPoints);
00178 mPoints.reserve(mNumQuadPoints);
00179
00180 switch (numPointsInEachDimension)
00181 {
00182 case 1:
00183 mWeights.push_back(0.12500000000000);
00184 mPoints.push_back(ChastePoint<3>(0.25000000000000,0.50000000000000,0.12500000000000));
00185 break;
00186
00187 case 2:
00188 mWeights.push_back(0.06132032652029);
00189 mWeights.push_back(0.01643073197073);
00190 mWeights.push_back(0.00440260136261);
00191 mWeights.push_back(0.00117967347971);
00192 mWeights.push_back(0.06132032652029);
00193 mWeights.push_back(0.01643073197073);
00194 mWeights.push_back(0.00440260136261);
00195 mWeights.push_back(0.00117967347971);
00196
00197 mPoints.push_back(ChastePoint<3>(0.16666666666667, 0.21132486540519, 0.13144585576580));
00198 mPoints.push_back(ChastePoint<3>(0.62200846792815, 0.21132486540519, 0.03522081090086));
00199 mPoints.push_back(ChastePoint<3>(0.04465819873852, 0.78867513459481, 0.03522081090086));
00200 mPoints.push_back(ChastePoint<3>(0.16666666666667, 0.78867513459481, 0.00943738783766));
00201 mPoints.push_back(ChastePoint<3>(0.16666666666667, 0.21132486540519, 0.49056261216234));
00202 mPoints.push_back(ChastePoint<3>(0.62200846792815, 0.21132486540519, 0.13144585576580));
00203 mPoints.push_back(ChastePoint<3>(0.04465819873852, 0.78867513459481, 0.13144585576580));
00204 mPoints.push_back(ChastePoint<3>(0.16666666666667, 0.78867513459481, 0.03522081090086));
00205 break;
00206
00207 case 3:
00208 mWeights.push_back(0.01497274736603);
00209 mWeights.push_back(0.01349962850795);
00210 mWeights.push_back(0.00190178826891);
00211 mWeights.push_back(0.00760715307442);
00212 mWeights.push_back(0.00685871056241);
00213 mWeights.push_back(0.00096623512860);
00214 mWeights.push_back(0.00024155878219);
00215 mWeights.push_back(0.00021779261632);
00216 mWeights.push_back(0.00003068198821);
00217 mWeights.push_back(0.02395639578565);
00218 mWeights.push_back(0.02159940561273);
00219 mWeights.push_back(0.00304286123026);
00220 mWeights.push_back(0.01217144491907);
00221 mWeights.push_back(0.01097393689986);
00222 mWeights.push_back(0.00154597620576);
00223 mWeights.push_back(0.00038649405150);
00224 mWeights.push_back(0.00034846818612);
00225 mWeights.push_back(0.00004909118114);
00226 mWeights.push_back(0.01497274736603);
00227 mWeights.push_back(0.01349962850795);
00228 mWeights.push_back(0.00190178826891);
00229 mWeights.push_back(0.00760715307442);
00230 mWeights.push_back(0.00685871056241);
00231 mWeights.push_back(0.00096623512860);
00232 mWeights.push_back(0.00024155878219);
00233 mWeights.push_back(0.00021779261632);
00234 mWeights.push_back(0.00003068198821);
00235
00236 mPoints.push_back(ChastePoint<3>(0.10000000001607, 0.11270166540000, 0.08872983347426));
00237 mPoints.push_back(ChastePoint<3>(0.44364916730000, 0.11270166540000, 0.05000000000803));
00238 mPoints.push_back(ChastePoint<3>(0.78729833458393, 0.11270166540000, 0.01127016654181));
00239 mPoints.push_back(ChastePoint<3>(0.05635083270000, 0.50000000000000, 0.05000000000803));
00240 mPoints.push_back(ChastePoint<3>(0.25000000000000, 0.50000000000000, 0.02817541635000));
00241 mPoints.push_back(ChastePoint<3>(0.44364916730000, 0.50000000000000, 0.00635083269197));
00242 mPoints.push_back(ChastePoint<3>(0.01270166538393, 0.88729833460000, 0.01127016654181));
00243 mPoints.push_back(ChastePoint<3>(0.05635083270000, 0.88729833460000, 0.00635083269197));
00244 mPoints.push_back(ChastePoint<3>(0.10000000001607, 0.88729833460000, 0.00143149884212));
00245 mPoints.push_back(ChastePoint<3>(0.10000000001607, 0.11270166540000, 0.39364916729197));
00246 mPoints.push_back(ChastePoint<3>(0.44364916730000, 0.11270166540000, 0.22182458365000));
00247 mPoints.push_back(ChastePoint<3>(0.78729833458393, 0.11270166540000, 0.05000000000803));
00248 mPoints.push_back(ChastePoint<3>(0.05635083270000, 0.50000000000000, 0.22182458365000));
00249 mPoints.push_back(ChastePoint<3>(0.25000000000000, 0.50000000000000, 0.12500000000000));
00250 mPoints.push_back(ChastePoint<3>(0.44364916730000, 0.50000000000000, 0.02817541635000));
00251 mPoints.push_back(ChastePoint<3>(0.01270166538393, 0.88729833460000, 0.05000000000803));
00252 mPoints.push_back(ChastePoint<3>(0.05635083270000, 0.88729833460000, 0.02817541635000));
00253 mPoints.push_back(ChastePoint<3>(0.10000000001607, 0.88729833460000, 0.00635083269197));
00254 mPoints.push_back(ChastePoint<3>(0.10000000001607, 0.11270166540000, 0.69856850110968));
00255 mPoints.push_back(ChastePoint<3>(0.44364916730000, 0.11270166540000, 0.39364916729197));
00256 mPoints.push_back(ChastePoint<3>(0.78729833458393, 0.11270166540000, 0.08872983347426));
00257 mPoints.push_back(ChastePoint<3>(0.05635083270000, 0.50000000000000, 0.39364916729197));
00258 mPoints.push_back(ChastePoint<3>(0.25000000000000, 0.50000000000000, 0.22182458365000));
00259 mPoints.push_back(ChastePoint<3>(0.44364916730000, 0.50000000000000, 0.05000000000803));
00260 mPoints.push_back(ChastePoint<3>(0.01270166538393, 0.88729833460000, 0.08872983347426));
00261 mPoints.push_back(ChastePoint<3>(0.05635083270000, 0.88729833460000, 0.05000000000803));
00262 mPoints.push_back(ChastePoint<3>(0.10000000001607, 0.88729833460000, 0.01127016654181));
00263 break;
00264
00265 case 4:
00266 mWeights.push_back(0.00423982561968);
00267 mWeights.push_back(0.00572288385156);
00268 mWeights.push_back(0.00281885467361);
00269 mWeights.push_back(0.00031634320391);
00270 mWeights.push_back(0.00412036229051);
00271 mWeights.push_back(0.00556163317318);
00272 mWeights.push_back(0.00273942929295);
00273 mWeights.push_back(0.00030742976838);
00274 mWeights.push_back(0.00099965677330);
00275 mWeights.push_back(0.00134932898618);
00276 mWeights.push_back(0.00066462336430);
00277 mWeights.push_back(0.00007458670588);
00278 mWeights.push_back(0.00002360309872);
00279 mWeights.push_back(0.00003185928022);
00280 mWeights.push_back(0.00001569255698);
00281 mWeights.push_back(0.00000176108183);
00282 mWeights.push_back(0.00794866986669);
00283 mWeights.push_back(0.01072905315027);
00284 mWeights.push_back(0.00528468555374);
00285 mWeights.push_back(0.00059306865848);
00286 mWeights.push_back(0.00772470439029);
00287 mWeights.push_back(0.01042674628127);
00288 mWeights.push_back(0.00513578175757);
00289 mWeights.push_back(0.00057635807584);
00290 mWeights.push_back(0.00187411992466);
00291 mWeights.push_back(0.00252967258912);
00292 mWeights.push_back(0.00124601155388);
00293 mWeights.push_back(0.00013983242583);
00294 mWeights.push_back(0.00004425022545);
00295 mWeights.push_back(0.00005972861231);
00296 mWeights.push_back(0.00002941983138);
00297 mWeights.push_back(0.00000330161175);
00298 mWeights.push_back(0.00794866986669);
00299 mWeights.push_back(0.01072905315027);
00300 mWeights.push_back(0.00528468555374);
00301 mWeights.push_back(0.00059306865848);
00302 mWeights.push_back(0.00772470439029);
00303 mWeights.push_back(0.01042674628127);
00304 mWeights.push_back(0.00513578175757);
00305 mWeights.push_back(0.00057635807584);
00306 mWeights.push_back(0.00187411992466);
00307 mWeights.push_back(0.00252967258912);
00308 mWeights.push_back(0.00124601155388);
00309 mWeights.push_back(0.00013983242583);
00310 mWeights.push_back(0.00004425022545);
00311 mWeights.push_back(0.00005972861231);
00312 mWeights.push_back(0.00002941983138);
00313 mWeights.push_back(0.00000330161175);
00314 mWeights.push_back(0.00423982561968);
00315 mWeights.push_back(0.00572288385156);
00316 mWeights.push_back(0.00281885467361);
00317 mWeights.push_back(0.00031634320391);
00318 mWeights.push_back(0.00412036229051);
00319 mWeights.push_back(0.00556163317318);
00320 mWeights.push_back(0.00273942929295);
00321 mWeights.push_back(0.00030742976838);
00322 mWeights.push_back(0.00099965677330);
00323 mWeights.push_back(0.00134932898618);
00324 mWeights.push_back(0.00066462336430);
00325 mWeights.push_back(0.00007458670588);
00326 mWeights.push_back(0.00002360309872);
00327 mWeights.push_back(0.00003185928022);
00328 mWeights.push_back(0.00001569255698);
00329 mWeights.push_back(0.00000176108183);
00330
00331 mPoints.push_back(ChastePoint<3>(0.06461106321099, 0.06943184420000, 0.06012499793653));
00332 mPoints.push_back(ChastePoint<3>(0.30709631152509, 0.06943184420000, 0.04328879995478));
00333 mPoints.push_back(ChastePoint<3>(0.62347184427491, 0.06943184420000, 0.02132226325621));
00334 mPoints.push_back(ChastePoint<3>(0.86595709258901, 0.06943184420000, 0.00448606527446));
00335 mPoints.push_back(ChastePoint<3>(0.04651867752509, 0.33000947820000, 0.04328879995478));
00336 mPoints.push_back(ChastePoint<3>(0.22110322249816, 0.33000947820000, 0.03116707302848));
00337 mPoints.push_back(ChastePoint<3>(0.44888729930184, 0.33000947820000, 0.01535160449661));
00338 mPoints.push_back(ChastePoint<3>(0.62347184427491, 0.33000947820000, 0.00322987757031));
00339 mPoints.push_back(ChastePoint<3>(0.02291316667491, 0.66999052180000, 0.02132226325621));
00340 mPoints.push_back(ChastePoint<3>(0.10890625570184, 0.66999052180000, 0.01535160449661));
00341 mPoints.push_back(ChastePoint<3>(0.22110322249816, 0.66999052180000, 0.00756156217830));
00342 mPoints.push_back(ChastePoint<3>(0.30709631152509, 0.66999052180000, 0.00159090341870));
00343 mPoints.push_back(ChastePoint<3>(0.00482078098901, 0.93056815580000, 0.00448606527446));
00344 mPoints.push_back(ChastePoint<3>(0.02291316667491, 0.93056815580000, 0.00322987757031));
00345 mPoints.push_back(ChastePoint<3>(0.04651867752509, 0.93056815580000, 0.00159090341870));
00346 mPoints.push_back(ChastePoint<3>(0.06461106321099, 0.93056815580000, 0.00033471571455));
00347 mPoints.push_back(ChastePoint<3>(0.06461106321099, 0.06943184420000, 0.28577404826889));
00348 mPoints.push_back(ChastePoint<3>(0.30709631152509, 0.06943184420000, 0.20575161800155));
00349 mPoints.push_back(ChastePoint<3>(0.62347184427491, 0.06943184420000, 0.10134469352354));
00350 mPoints.push_back(ChastePoint<3>(0.86595709258901, 0.06943184420000, 0.02132226325621));
00351 mPoints.push_back(ChastePoint<3>(0.04651867752509, 0.33000947820000, 0.20575161800155));
00352 mPoints.push_back(ChastePoint<3>(0.22110322249816, 0.33000947820000, 0.14813706341321));
00353 mPoints.push_back(ChastePoint<3>(0.44888729930184, 0.33000947820000, 0.07296615908496));
00354 mPoints.push_back(ChastePoint<3>(0.62347184427491, 0.33000947820000, 0.01535160449661));
00355 mPoints.push_back(ChastePoint<3>(0.02291316667491, 0.66999052180000, 0.10134469352354));
00356 mPoints.push_back(ChastePoint<3>(0.10890625570184, 0.66999052180000, 0.07296615908496));
00357 mPoints.push_back(ChastePoint<3>(0.22110322249816, 0.66999052180000, 0.03594009661688));
00358 mPoints.push_back(ChastePoint<3>(0.30709631152509, 0.66999052180000, 0.00756156217830));
00359 mPoints.push_back(ChastePoint<3>(0.00482078098901, 0.93056815580000, 0.02132226325621));
00360 mPoints.push_back(ChastePoint<3>(0.02291316667491, 0.93056815580000, 0.01535160449661));
00361 mPoints.push_back(ChastePoint<3>(0.04651867752509, 0.93056815580000, 0.00756156217830));
00362 mPoints.push_back(ChastePoint<3>(0.06461106321099, 0.93056815580000, 0.00159090341870));
00363 mPoints.push_back(ChastePoint<3>(0.06461106321099, 0.06943184420000, 0.58018304432012));
00364 mPoints.push_back(ChastePoint<3>(0.30709631152509, 0.06943184420000, 0.41772022627335));
00365 mPoints.push_back(ChastePoint<3>(0.62347184427491, 0.06943184420000, 0.20575161800155));
00366 mPoints.push_back(ChastePoint<3>(0.86595709258901, 0.06943184420000, 0.04328879995478));
00367 mPoints.push_back(ChastePoint<3>(0.04651867752509, 0.33000947820000, 0.41772022627335));
00368 mPoints.push_back(ChastePoint<3>(0.22110322249816, 0.33000947820000, 0.30075023588863));
00369 mPoints.push_back(ChastePoint<3>(0.44888729930184, 0.33000947820000, 0.14813706341321));
00370 mPoints.push_back(ChastePoint<3>(0.62347184427491, 0.33000947820000, 0.03116707302848));
00371 mPoints.push_back(ChastePoint<3>(0.02291316667491, 0.66999052180000, 0.20575161800155));
00372 mPoints.push_back(ChastePoint<3>(0.10890625570184, 0.66999052180000, 0.14813706341321));
00373 mPoints.push_back(ChastePoint<3>(0.22110322249816, 0.66999052180000, 0.07296615908496));
00374 mPoints.push_back(ChastePoint<3>(0.30709631152509, 0.66999052180000, 0.01535160449661));
00375 mPoints.push_back(ChastePoint<3>(0.00482078098901, 0.93056815580000, 0.04328879995478));
00376 mPoints.push_back(ChastePoint<3>(0.02291316667491, 0.93056815580000, 0.03116707302848));
00377 mPoints.push_back(ChastePoint<3>(0.04651867752509, 0.93056815580000, 0.01535160449661));
00378 mPoints.push_back(ChastePoint<3>(0.06461106321099, 0.93056815580000, 0.00322987757031));
00379 mPoints.push_back(ChastePoint<3>(0.06461106321099, 0.06943184420000, 0.80583209465249));
00380 mPoints.push_back(ChastePoint<3>(0.30709631152509, 0.06943184420000, 0.58018304432012));
00381 mPoints.push_back(ChastePoint<3>(0.62347184427491, 0.06943184420000, 0.28577404826889));
00382 mPoints.push_back(ChastePoint<3>(0.86595709258901, 0.06943184420000, 0.06012499793653));
00383 mPoints.push_back(ChastePoint<3>(0.04651867752509, 0.33000947820000, 0.58018304432012));
00384 mPoints.push_back(ChastePoint<3>(0.22110322249816, 0.33000947820000, 0.41772022627335));
00385 mPoints.push_back(ChastePoint<3>(0.44888729930184, 0.33000947820000, 0.20575161800155));
00386 mPoints.push_back(ChastePoint<3>(0.62347184427491, 0.33000947820000, 0.04328879995478));
00387 mPoints.push_back(ChastePoint<3>(0.02291316667491, 0.66999052180000, 0.28577404826889));
00388 mPoints.push_back(ChastePoint<3>(0.10890625570184, 0.66999052180000, 0.20575161800155));
00389 mPoints.push_back(ChastePoint<3>(0.22110322249816, 0.66999052180000, 0.10134469352354));
00390 mPoints.push_back(ChastePoint<3>(0.30709631152509, 0.66999052180000, 0.02132226325621));
00391 mPoints.push_back(ChastePoint<3>(0.00482078098901, 0.93056815580000, 0.06012499793653));
00392 mPoints.push_back(ChastePoint<3>(0.02291316667491, 0.93056815580000, 0.04328879995478));
00393 mPoints.push_back(ChastePoint<3>(0.04651867752509, 0.93056815580000, 0.02132226325621));
00394 mPoints.push_back(ChastePoint<3>(0.06461106321099, 0.93056815580000, 0.00448606527446));
00395 break;
00396
00397 default:
00398 EXCEPTION("Number of gauss points per dimension not supported.");
00399 }
00400 }
00401
00402 template<unsigned ELEMENT_DIM>
00403 GaussianQuadratureRule<ELEMENT_DIM>::GaussianQuadratureRule(unsigned numPointsInEachDimension)
00404 {
00405 EXCEPTION("Gauss points not available for this dimension.");
00406 }
00407
00408
00410
00412
00413 template class GaussianQuadratureRule<0>;
00414 template class GaussianQuadratureRule<1>;
00415 template class GaussianQuadratureRule<2>;
00416 template class GaussianQuadratureRule<3>;
00417 template class GaussianQuadratureRule<4>;