1!/*****************************************************************************/ 2! * 3! * Elmer, A Finite Element Software for Multiphysical Problems 4! * 5! * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland 6! * 7! * This library is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU Lesser General Public 9! * License as published by the Free Software Foundation; either 10! * version 2.1 of the License, or (at your option) any later version. 11! * 12! * This library is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15! * Lesser General Public License for more details. 16! * 17! * You should have received a copy of the GNU Lesser General Public 18! * License along with this library (in file ../LGPL-2.1); if not, write 19! * to the Free Software Foundation, Inc., 51 Franklin Street, 20! * Fifth Floor, Boston, MA 02110-1301 USA 21! * 22! *****************************************************************************/ 23! 24!/****************************************************************************** 25! * 26! * Authors: Juha Ruokolainen 27! * Email: Juha.Ruokolainen@csc.fi 28! * Web: http://www.csc.fi/elmer 29! * Address: CSC - IT Center for Science Ltd. 30! * Keilaranta 14 31! * 02101 Espoo, Finland 32! * 33! * Original Date: 01 Oct 1996 34! * 35! ****************************************************************************/ 36 37!> \ingroup ElmerLib 38!> \{ 39 40!----------------------------------------------------------------------------- 41!> Module defining Gauss integration points and 42!> containing various integration routines. 43!----------------------------------------------------------------------------- 44MODULE Integration 45 USE Types 46 USE LoadMod 47 USE LinearAlgebra 48 49 IMPLICIT NONE 50 51 INTEGER, PARAMETER, PRIVATE :: MAXN = 13, MAXNPAD = 16 ! Padded to 64-byte alignment 52 INTEGER, PARAMETER, PRIVATE :: MAX_INTEGRATION_POINTS = MAXN**3 53 54 LOGICAL, PRIVATE :: GInit = .FALSE. 55 !$OMP THREADPRIVATE(GInit) 56 57!------------------------------------------------------------------------------ 58 TYPE GaussIntegrationPoints_t 59 INTEGER :: N 60 REAL(KIND=dp), POINTER CONTIG :: u(:),v(:),w(:),s(:) 61!DIR$ ATTRIBUTES ALIGN:64 :: u, v, w, s 62 END TYPE GaussIntegrationPoints_t 63 64 TYPE(GaussIntegrationPoints_t), TARGET, PRIVATE, SAVE :: IntegStuff 65 !$OMP THREADPRIVATE(IntegStuff) 66!------------------------------------------------------------------------------ 67 68!------------------------------------------------------------------------------ 69! Storage for 1d Gauss points, and weights. The values are computed on the 70! fly (see ComputeGaussPoints1D below). These values are used for quads and 71! bricks as well. To avoid NUMA issues, Points and Weights are private for each 72! thread. 73!------------------------------------------------------------------------------ 74 REAL(KIND=dp), PRIVATE, SAVE :: Points(MAXNPAD,MAXN),Weights(MAXNPAD,MAXN) 75 !$OMP THREADPRIVATE(Points, Weights) 76!DIR$ ATTRIBUTES ALIGN:64::Points, Weights 77!------------------------------------------------------------------------------ 78 79!------------------------------------------------------------------------------ 80! Triangle - 1 point rule; exact integration of x^py^q, p+q<=1 81!------------------------------------------------------------------------------ 82 REAL(KIND=dp),DIMENSION(1),PRIVATE :: UTriangle1P=(/ 0.3333333333333333D0 /) 83 REAL(KIND=dp),DIMENSION(1),PRIVATE :: VTriangle1P=(/ 0.3333333333333333D0 /) 84 REAL(KIND=dp),DIMENSION(1),PRIVATE :: STriangle1P=(/ 1.0000000000000000D0 /) 85 86!------------------------------------------------------------------------------ 87! Triangle - 3 point rule; exact integration of x^py^q, p+q<=2 88!------------------------------------------------------------------------------ 89 REAL(KIND=dp), DIMENSION(3), PRIVATE :: UTriangle3P = & 90 (/ 0.16666666666666667D0, 0.66666666666666667D0, 0.16666666666666667D0 /) 91 92 REAL(KIND=dp), DIMENSION(3), PRIVATE :: VTriangle3P = & 93 (/ 0.16666666666666667D0, 0.16666666666666667D0, 0.66666666666666667D0 /) 94 95 REAL(KIND=dp), DIMENSION(3), PRIVATE :: STriangle3P = & 96 (/ 0.33333333333333333D0, 0.33333333333333333D0, 0.33333333333333333D0 /) 97 98!------------------------------------------------------------------------------ 99! Triangle - 4 point rule; exact integration of x^py^q, p+q<=3 100!------------------------------------------------------------------------------ 101 REAL(KIND=dp), DIMENSION(4), PRIVATE :: UTriangle4P = & 102 (/ 0.33333333333333333D0, 0.2000000000000000D0, & 103 0.60000000000000000D0, 0.20000000000000000D0 /) 104 105 REAL(KIND=dp), DIMENSION(4), PRIVATE :: VTriangle4P = & 106 (/ 0.33333333333333333D0, 0.2000000000000000D0, & 107 0.20000000000000000D0, 0.60000000000000000D0 /) 108 109 REAL(KIND=dp), DIMENSION(4), PRIVATE :: STriangle4P = & 110 (/ -0.56250000000000000D0, 0.52083333333333333D0, & 111 0.52083333333333333D0, 0.52083333333333333D0 /) 112 113!------------------------------------------------------------------------------ 114! Triangle - 6 point rule; exact integration of x^py^q, p+q<=4 115!------------------------------------------------------------------------------ 116 REAL(KIND=dp), DIMENSION(6), PRIVATE :: UTriangle6P = & 117 (/ 0.091576213509771D0, 0.816847572980459D0, 0.091576213509771D0, & 118 0.445948490915965D0, 0.108103018168070D0, 0.445948490915965D0 /) 119 120 REAL(KIND=dp), DIMENSION(6), PRIVATE :: VTriangle6P = & 121 (/ 0.091576213509771D0, 0.091576213509771D0, 0.816847572980459D0, & 122 0.445948490915965D0, 0.445948490915965D0, 0.108103018168070D0 /) 123 124 REAL(KIND=dp), DIMENSION(6), PRIVATE :: STriangle6P = & 125 (/ 0.109951743655322D0, 0.109951743655322D0, 0.109951743655322D0, & 126 0.223381589678011D0, 0.223381589678011D0, 0.223381589678011D0 /) 127 128!------------------------------------------------------------------------------ 129! Triangle - 7 point rule; exact integration of x^py^q, p+q<=5 130!------------------------------------------------------------------------------ 131 REAL(KIND=dp), DIMENSION(7), PRIVATE :: UTriangle7P = & 132 (/ 0.333333333333333D0, 0.101286507323456D0, 0.797426985353087D0, & 133 0.101286507323456D0, 0.470142064105115D0, 0.059715871789770D0, & 134 0.470142064105115D0 /) 135 136 REAL(KIND=dp), DIMENSION(7), PRIVATE :: VTriangle7P = & 137 (/ 0.333333333333333D0, 0.101286507323456D0, 0.101286507323456D0, & 138 0.797426985353087D0, 0.470142064105115D0, 0.470142064105115D0, & 139 0.059715871789770D0 /) 140 141 REAL(KIND=dp), DIMENSION(7), PRIVATE :: STriangle7P = & 142 (/ 0.225000000000000D0, 0.125939180544827D0, 0.125939180544827D0, & 143 0.125939180544827D0, 0.132394152788506D0, 0.132394152788506D0, & 144 0.132394152788506D0 /) 145 146!------------------------------------------------------------------------------ 147! Triangle - 11 point rule; exact integration of x^py^q, p+q<=6 148!------------------------------------------------------------------------------ 149 REAL(KIND=dp), DIMENSION(11), PRIVATE :: UTriangle11P = & 150 (/ 0.3019427231413448D-01, 0.5298143569082113D-01, & 151 0.4972454892773975D-01, 0.7697772693248785D-01, & 152 0.7008117469890058D+00, 0.5597774797709894D+00, & 153 0.5428972301980696D+00, 0.3437947421925572D+00, & 154 0.2356669356664465D+00, 0.8672623210691472D+00, & 155 0.2151020995173866D+00 /) 156 157 REAL(KIND=dp), DIMENSION(11), PRIVATE :: VTriangle11P = & 158 (/ 0.2559891985673773D+00, 0.1748087863744473D-01, & 159 0.6330812033358987D+00, 0.8588528075577063D+00, & 160 0.2708520519075563D+00, 0.1870602768957014D-01, & 161 0.2027008533579804D+00, 0.5718576583152437D+00, & 162 0.1000777578531811D+00, 0.4654861310422605D-01, & 163 0.3929681357810497D+00 /) 164 165 REAL(KIND=dp), DIMENSION(11), PRIVATE :: STriangle11P = & 166 (/ 0.3375321205342688D-01, 0.1148426034648707D-01, & 167 0.4197958777582435D-01, 0.3098130358202468D-01, & 168 0.2925899761167147D-01, 0.2778515729102349D-01, & 169 0.8323049608963519D-01, 0.6825761580824108D-01, & 170 0.6357334991651026D-01, 0.2649352562792455D-01, & 171 0.8320249389723097D-01 /) 172 173!------------------------------------------------------------------------------ 174! Triangle - 12 point rule; exact integration of x^py^q, p+q<=7 175!------------------------------------------------------------------------------ 176 REAL(KIND=dp), DIMENSION(12), PRIVATE :: UTriangle12P = & 177 (/ 0.6232720494911090D+00, 0.3215024938520235D+00, & 178 0.5522545665686063D-01, 0.2777161669760336D+00, & 179 0.5158423343535236D+00, 0.2064414986704435D+00, & 180 0.3432430294535058D-01, 0.3047265008682535D+00, & 181 0.6609491961864082D+00, 0.6238226509441210D-01, & 182 0.8700998678316921D+00, 0.6751786707389329D-01 /) 183 184 REAL(KIND=dp), DIMENSION(12), PRIVATE :: VTriangle12P = & 185 (/ 0.3215024938520235D+00, 0.5522545665686063D-01, & 186 0.6232720494911090D+00, 0.5158423343535236D+00, & 187 0.2064414986704435D+00, 0.2777161669760336D+00, & 188 0.3047265008682535D+00, 0.6609491961864082D+00, & 189 0.3432430294535058D-01, 0.8700998678316921D+00, & 190 0.6751786707389329D-01, 0.6238226509441210D-01 /) 191 192 REAL(KIND=dp), DIMENSION(12), PRIVATE :: STriangle12P = & 193 (/ 0.4388140871440586D-01, 0.4388140871440586D-01, & 194 0.4388140871440587D-01, 0.6749318700971417D-01, & 195 0.6749318700971417D-01, 0.6749318700971417D-01, & 196 0.2877504278510970D-01, 0.2877504278510970D-01, & 197 0.2877504278510969D-01, 0.2651702815743698D-01, & 198 0.2651702815743698D-01, 0.2651702815743698D-01 /) 199 200!------------------------------------------------------------------------------ 201! Triangle - 17 point rule; exact integration of x^py^q, p+q<=8 202!------------------------------------------------------------------------------ 203 REAL(KIND=dp), DIMENSION(17), PRIVATE :: UTriangle17P = & 204 (/ 0.2292423642627924D+00, 0.4951220175479885D-01, & 205 0.3655948407066446D+00, 0.4364350639589269D+00, & 206 0.1596405673569602D+00, 0.9336507149305228D+00, & 207 0.5219569066777245D+00, 0.7110782758797098D+00, & 208 0.5288509041694864D+00, 0.1396967677642513D-01, & 209 0.4205421906708996D-01, 0.4651359156686354D-01, & 210 0.1975981349257204D+00, 0.7836841874017514D+00, & 211 0.4232808751402256D-01, 0.4557097415216423D+00, & 212 0.2358934246935281D+00 /) 213 214 REAL(KIND=dp), DIMENSION(17), PRIVATE :: VTriangle17P = & 215 (/ 0.5117407211006358D+00, 0.7589103637479163D+00, & 216 0.1529647481767193D+00, 0.3151398735074337D-01, & 217 0.5117868393288316D-01, 0.1964516824106966D-01, & 218 0.2347490459725670D+00, 0.4908682577187765D-01, & 219 0.4382237537321878D+00, 0.3300210677033395D-01, & 220 0.2088758614636060D+00, 0.9208929246654702D+00, & 221 0.2742740954674795D+00, 0.1654585179097472D+00, & 222 0.4930011699833554D+00, 0.4080804967846944D+00, & 223 0.7127872162741824D+00 /) 224 225 REAL(KIND=dp), DIMENSION(17), PRIVATE :: STriangle17P = & 226 (/ 0.5956595662857148D-01, 0.2813390230006461D-01, & 227 0.3500735477096827D-01, 0.2438077450393263D-01, & 228 0.2843374448051010D-01, 0.7822856634218779D-02, & 229 0.5179111341004783D-01, 0.3134229539096806D-01, & 230 0.2454951584925144D-01, 0.5371382557647114D-02, & 231 0.2571565514768072D-01, 0.1045933340802507D-01, & 232 0.4937780841212319D-01, 0.2824772362317584D-01, & 233 0.3218881684015661D-01, 0.2522089247693226D-01, & 234 0.3239087356572598D-01 /) 235 236!------------------------------------------------------------------------------ 237! Triangle - 20 point rule; exact integration of x^py^q, p+q<=9 238!------------------------------------------------------------------------------ 239 REAL(KIND=dp), DIMENSION(20), PRIVATE :: UTriangle20P = & 240 (/ 0.2469118866487856D-01, 0.3348782965514246D-00, & 241 0.4162560937597861D-00, 0.1832492889417431D-00, & 242 0.2183952668281443D-00, 0.4523362527443628D-01, & 243 0.4872975112073226D-00, 0.7470127381316580D-00, & 244 0.7390287107520658D-00, 0.3452260444515281D-01, & 245 0.4946745467572288D-00, 0.3747439678780460D-01, & 246 0.2257524791528391D-00, 0.9107964437563798D-00, & 247 0.4254445629399445D-00, 0.1332215072275240D-00, & 248 0.5002480151788234D-00, 0.4411517722238269D-01, & 249 0.1858526744057914D-00, 0.6300024376672695D-00 /) 250 251 REAL(KIND=dp), DIMENSION(20), PRIVATE :: VTriangle20P = & 252 (/ 0.4783451248176442D-00, 0.3373844236168988D-00, & 253 0.1244378463254732D-00, 0.6365569723648120D-00, & 254 0.3899759363237886D-01, 0.9093437140096456D-00, & 255 0.1968266037596590D-01, 0.2191311129347709D-00, & 256 0.3833588560240875D-01, 0.7389795063475102D-00, & 257 0.4800989285800525D-00, 0.2175137165318176D-00, & 258 0.7404716820879975D-00, 0.4413531509926682D-01, & 259 0.4431292142978816D-00, 0.4440953593652837D-00, & 260 0.1430831401051367D-00, 0.4392970158411878D-01, & 261 0.1973209364545017D-00, 0.1979381059170009D-00 /) 262 263 REAL(KIND=dp), DIMENSION(20), PRIVATE :: STriangle20P = & 264 (/ 0.1776913091122958D-01, 0.4667544936904065D-01, & 265 0.2965283331432967D-01, 0.3880447634997608D-01, & 266 0.2251511457011248D-01, 0.1314162394636178D-01, & 267 0.1560341736610505D-01, 0.1967065434689744D-01, & 268 0.2247962849501080D-01, 0.2087108394969067D-01, & 269 0.1787661200700672D-01, 0.2147695865607915D-01, & 270 0.2040998247303970D-01, 0.1270342300533680D-01, & 271 0.3688713099356314D-01, 0.3813199811535777D-01, & 272 0.1508642325812160D-01, 0.1238422287692121D-01, & 273 0.3995072336992735D-01, 0.3790911262589247D-01 /) 274 275!------------------------------------------------------------------------------ 276! Tetrahedron - 1 point rule; exact integration of x^py^qz^r, p+q+r<=1 277!------------------------------------------------------------------------------ 278 REAL(KIND=dp), DIMENSION(1), PRIVATE :: UTetra1P = (/ 0.25D0 /) 279 REAL(KIND=dp), DIMENSION(1), PRIVATE :: VTetra1P = (/ 0.25D0 /) 280 REAL(KIND=dp), DIMENSION(1), PRIVATE :: WTetra1P = (/ 0.25D0 /) 281 REAL(KIND=dp), DIMENSION(1), PRIVATE :: STetra1P = (/ 1.00D0 /) 282 283!------------------------------------------------------------------------------ 284! Tetrahedron - 4 point rule; exact integration of x^py^qz^r, p+q+r<=2 285!------------------------------------------------------------------------------ 286 REAL(KIND=dp), DIMENSION(4), PRIVATE :: UTetra4P = & 287 (/ 0.1757281246520584D0, 0.2445310270213291D0, & 288 0.5556470949048655D0, 0.0240937534217468D0 /) 289 290 REAL(KIND=dp), DIMENSION(4), PRIVATE :: VTetra4P = & 291 (/ 0.5656137776620919D0, 0.0501800797762026D0, & 292 0.1487681308666864D0, 0.2354380116950194D0 /) 293 294 REAL(KIND=dp), DIMENSION(4), PRIVATE :: WTetra4P = & 295 (/ 0.2180665126782654D0, 0.5635595064952189D0, & 296 0.0350112499848832D0, 0.1833627308416330D0 /) 297 298 REAL(KIND=dp), DIMENSION(4), PRIVATE :: STetra4P = & 299 (/ 0.2500000000000000D0, 0.2500000000000000D0, & 300 0.2500000000000000D0, 0.2500000000000000D0 /) 301!------------------------------------------------------------------------------ 302! Tetrahedron - 5 point rule; exact integration of x^py^qz^r, p+q+r<=3 303!------------------------------------------------------------------------------ 304 REAL(KIND=dp), DIMENSION(5), PRIVATE :: UTetra5P = & 305 (/ 0.25000000000000000D0, 0.50000000000000000D0, & 306 0.16666666666666667D0, 0.16666666666666667D0, & 307 0.16666666666666667D0 /) 308 309 REAL(KIND=dp), DIMENSION(5), PRIVATE :: VTetra5P = & 310 (/ 0.25000000000000000D0, 0.16666666666666667D0, & 311 0.50000000000000000D0, 0.16666666666666667D0, & 312 0.16666666666666667D0 /) 313 314 REAL(KIND=dp), DIMENSION(5), PRIVATE :: WTetra5P = & 315 (/ 0.25000000000000000D0, 0.16666666666666667D0, & 316 0.16666666666666667D0, 0.50000000000000000D0, & 317 0.16666666666666667D0 /) 318 319 REAL(KIND=dp), DIMENSION(5), PRIVATE :: STetra5P = & 320 (/-0.80000000000000000D0, 0.45000000000000000D0, & 321 0.45000000000000000D0, 0.45000000000000000D0, & 322 0.45000000000000000D0 /) 323!------------------------------------------------------------------------------ 324! Tetrahedron - 11 point rule; exact integration of x^py^qz^r, p+q+r<=4 325!------------------------------------------------------------------------------ 326 REAL(KIND=dp), DIMENSION(11), PRIVATE :: UTetra11P = & 327 (/ 0.3247902050850455D+00, 0.4381969657060433D+00, & 328 0.8992592373310454D-01, 0.1092714936292849D+00, & 329 0.3389119319942253D-01, 0.5332363613904868D-01, & 330 0.1935618747806815D+00, 0.4016250624424964D-01, & 331 0.3878132182319405D+00, 0.7321489692875428D+00, & 332 0.8066342495294049D-01 /) 333 334 REAL(KIND=dp), DIMENSION(11), PRIVATE :: VTetra11P = & 335 (/ 0.4573830181783998D+00, 0.9635325047480842D-01, & 336 0.3499588148445295D+00, 0.1228957438582778D+00, & 337 0.4736224692062527D-01, 0.4450376952468180D+00, & 338 0.2165626476982170D+00, 0.8033385922433729D+00, & 339 0.7030897281814283D-01, 0.1097836536360084D+00, & 340 0.1018859284267242D-01 /) 341 342 REAL(KIND=dp), DIMENSION(11), PRIVATE :: WTetra11P = & 343 (/ 0.1116787541193331D+00, 0.6966288385119494D-01, & 344 0.5810783971325720D-01, 0.3424607753785182D+00, & 345 0.7831772466208499D+00, 0.3688112094344830D+00, & 346 0.5872345323698884D+00, 0.6178518963560731D-01, & 347 0.4077342860913465D+00, 0.9607290317342082D-01, & 348 0.8343823045787845D-01 /) 349 350 REAL(KIND=dp), DIMENSION(11), PRIVATE :: STetra11P = & 351 (/ 0.1677896627448221D+00, 0.1128697325878004D+00, & 352 0.1026246621329828D+00, 0.1583002576888426D+00, & 353 0.3847841737508437D-01, 0.1061709382037234D+00, & 354 0.5458124994014422D-01, 0.3684475128738168D-01, & 355 0.1239234851349682D+00, 0.6832098141300300D-01, & 356 0.3009586149124714D-01 /) 357!------------------------------------------------------------------------------ 358 359!------------------------------------------------------------------------------ 360! A SELECTION OF QUADRATURE RULES UP TO A HIGH ORDER DESIGNED FOR ECONOMIC 361! INTEGRATION OF COMPLETE POLYNOMIALS 362! 363! These rules have been reproduced from the source "P. Solin, K. Segeth 364! and I. Dolezel: Higher-Order Finite Element Methods", 365! Chapman & Hall/CRC Press, 2003. 366!------------------------------------------------------------------------------ 367 368!------------------------------------------------------------------------------ 369! Quadrilateral - 8-point rule for complete polynomials of order p<=5 370!------------------------------------------------------------------------------ 371 REAL(KIND=dp), DIMENSION(8), PRIVATE :: UPQuad8 = & 372 (/ 0.683130051063973d0,-0.683130051063973d0, 0.000000000000000d0, & 373 0.000000000000000d0, 0.881917103688197d0, 0.881917103688197d0, & 374 -0.881917103688197d0,-0.881917103688197d0 /) 375 376 REAL(KIND=dp), DIMENSION(8), PRIVATE :: VPQuad8 = & 377 (/ 0.000000000000000d0, 0.000000000000000d0, 0.683130051063973d0, & 378 -0.683130051063973d0, 0.881917103688197d0,-0.881917103688197d0, & 379 0.881917103688197d0,-0.881917103688197d0 /) 380 381 REAL(KIND=dp), DIMENSION(8), PRIVATE :: SPQuad8 = & 382 (/ 0.816326530612245d0, 0.816326530612245d0, 0.816326530612245d0, & 383 0.816326530612245d0, 0.183673469387755d0, 0.183673469387755d0, & 384 0.183673469387755d0, 0.183673469387755d0 /) 385!------------------------------------------------------------------------------ 386! Quadrilateral - 12-point rule for complete polynomials of order p<=7 387!------------------------------------------------------------------------------ 388 REAL(KIND=dp), DIMENSION(12), PRIVATE :: UPQuad12 = & 389 (/ 0.925820099772551d0, -0.925820099772551d0, 0.000000000000000d0, & 390 0.000000000000000d0, 0.805979782918599d0, 0.805979782918599d0, & 391 -0.805979782918599d0, -0.805979782918599d0, 0.380554433208316d0, & 392 0.380554433208316d0, -0.380554433208316d0, -0.380554433208316d0 /) 393 394 REAL(KIND=dp), DIMENSION(12), PRIVATE :: VPQuad12 = & 395 (/ 0.000000000000000d0, 0.000000000000000d0, 0.925820099772551d0, & 396 -0.925820099772551d0, 0.805979782918599d0, -0.805979782918599d0, & 397 0.805979782918599d0, -0.805979782918599d0, 0.380554433208316d0, & 398 -0.380554433208316d0, 0.380554433208316d0, -0.380554433208316d0 /) 399 400 REAL(KIND=dp), DIMENSION(12), PRIVATE :: SPQuad12 = & 401 (/ 0.241975308641975d0, 0.241975308641975d0, 0.241975308641975d0, & 402 0.241975308641975d0, 0.237431774690630d0, 0.237431774690630d0, & 403 0.237431774690630d0, 0.237431774690630d0, 0.520592916667394d0, & 404 0.520592916667394d0, 0.520592916667394d0, 0.520592916667394d0 /) 405!------------------------------------------------------------------------------ 406! Quadrilateral - 20-point rule for complete polynomials of order p<=9 407! Note: Some points lie outside the reference element 408!------------------------------------------------------------------------------ 409 REAL(KIND=dp), DIMENSION(20), PRIVATE :: UPQuad20 = & 410 (/ 0.112122576386656d1, -0.112122576386656d1, 0.000000000000000d0, & 411 0.000000000000000d0, 0.451773049920657d0, -0.451773049920657d0, & 412 0.000000000000000d0, 0.000000000000000d0, 0.891849420851512d0, & 413 0.891849420851512d0, -0.891849420851512d0, -0.891849420851512d0, & 414 0.824396370749276d0, 0.824396370749276d0, -0.824396370749276d0, & 415 -0.824396370749276d0, 0.411623426336542d0, 0.411623426336542d0, & 416 -0.411623426336542d0, -0.411623426336542d0 /) 417 418 REAL(KIND=dp), DIMENSION(20), PRIVATE :: VPQuad20 = & 419 (/ 0.000000000000000d0, 0.000000000000000d0, 0.112122576386656d1, & 420 -0.112122576386656d1, 0.000000000000000d0, 0.000000000000000d0, & 421 0.451773049920657d0, -0.451773049920657d0, 0.891849420851512d0, & 422 -0.891849420851512d0, 0.891849420851512d0, -0.891849420851512d0, & 423 0.411623426336542d0, -0.411623426336542d0, 0.411623426336542d0, & 424 -0.411623426336542d0, 0.824396370749276d0, -0.824396370749276d0, & 425 0.824396370749276d0, -0.824396370749276d0 /) 426 427 REAL(KIND=dp), DIMENSION(20), PRIVATE :: SPQuad20 = & 428 (/ 0.184758425074910d-1, 0.184758425074910d-1, 0.184758425074910d-1, & 429 0.184758425074910d-1, 0.390052939160735d0, 0.390052939160735d0, & 430 0.390052939160735d0, 0.390052939160735d0, 0.830951780264820d-1, & 431 0.830951780264820d-1, 0.830951780264820d-1, 0.830951780264820d-1, & 432 0.254188020152646d0, 0.254188020152646d0, 0.254188020152646d0, & 433 0.254188020152646d0, 0.254188020152646d0, 0.254188020152646d0, & 434 0.254188020152646d0, 0.254188020152646d0 /) 435!------------------------------------------------------------------------------ 436! Quadrilateral - 25-point rule for complete polynomials of order p<=11 437! Note: Some points lie outside the reference element 438!------------------------------------------------------------------------------ 439 REAL(KIND=dp), DIMENSION(25), PRIVATE :: UPQuad25 = & 440 (/0.000000000000000d0, 0.104440291540981d1, -0.104440291540981d1, & 441 0.000000000000000d0, 0.000000000000000d0, 0.769799068396649d0, & 442 -0.769799068396649d0, 0.000000000000000d0, 0.000000000000000d0, & 443 0.935787012440540d0, 0.935787012440540d0, -0.935787012440540d0, & 444 -0.935787012440540d0, 0.413491953449114d0, 0.413491953449114d0, & 445 -0.413491953449114d0, -0.413491953449114d0, 0.883025508525690d0, & 446 0.883025508525690d0, -0.883025508525690d0, -0.883025508525690d0, & 447 0.575653595840465d0, 0.575653595840465d0, -0.575653595840465d0, & 448 -0.575653595840465d0 /) 449 450 REAL(KIND=dp), DIMENSION(25), PRIVATE :: VPQuad25 = & 451 (/ 0.000000000000000d0, 0.000000000000000d0, 0.000000000000000d0, & 452 0.104440291540981d1, -0.104440291540981d1, 0.000000000000000d0, & 453 0.000000000000000d0, 0.769799068396649d0, -0.769799068396649d0, & 454 0.935787012440540d0, -0.935787012440540d0, 0.935787012440540d0, & 455 -0.935787012440540d0, 0.413491953449114d0, -0.413491953449114d0, & 456 0.413491953449114d0, -0.413491953449114d0, 0.575653595840465d0, & 457 -0.575653595840465d0, 0.575653595840465d0, -0.575653595840465d0, & 458 0.883025508525690d0, -0.883025508525690d0, 0.883025508525690d0, & 459 -0.883025508525690d0 /) 460 461 REAL(KIND=dp), DIMENSION(25), PRIVATE :: SPQuad25 = & 462 (/ 0.365379525585903d0, 0.277561655642040d-1, 0.277561655642040d-1, & 463 0.277561655642040d-1, 0.277561655642040d-1, 0.244272057751754d0, & 464 0.244272057751754d0, 0.244272057751754d0, 0.244272057751754d0, & 465 0.342651038512290d-1, 0.342651038512290d-1, 0.342651038512290d-1, & 466 0.342651038512290d-1, 0.308993036133713d0, 0.308993036133713d0, & 467 0.308993036133713d0, 0.308993036133713d0, 0.146684377651312d0, & 468 0.146684377651312d0, 0.146684377651312d0, 0.146684377651312d0, & 469 0.146684377651312d0, 0.146684377651312d0, 0.146684377651312d0, & 470 0.146684377651312d0 /) 471!------------------------------------------------------------------------------ 472! Quadrilateral - 36-point rule for complete polynomials of order p<=13 473! Note 1: Some points lie outside the reference element 474! Note 2: Weights somewhat inaccurate? 475!------------------------------------------------------------------------------ 476 REAL(KIND=dp), DIMENSION(36), PRIVATE :: UPQuad36 = & 477 (/ 0.108605615857397d1, -0.108605615857397d1, 0.000000000000000d0, & 478 0.000000000000000d0, 0.658208197042585d0, -0.658208197042585d0, & 479 0.000000000000000d0, 0.000000000000000d0, 0.100130060299173d1, & 480 0.100130060299173d1, -0.100130060299173d1, -0.100130060299173d1, & 481 0.584636168775946d0, 0.584636168775946d0, -0.584636168775946d0, & 482 -0.584636168775946d0, 0.246795612720261d0, 0.246795612720261d0, & 483 -0.246795612720261d0, -0.246795612720261d0, 0.900258815287201d0, & 484 0.900258815287201d0, -0.900258815287201d0, -0.900258815287201d0, & 485 0.304720678579870d0, 0.304720678579870d0, -0.304720678579870d0, & 486 -0.304720678579870d0, 0.929866705560780d0, 0.929866705560780d0, & 487 -0.929866705560780d0, -0.929866705560780d0, 0.745052720131169d0, & 488 0.745052720131169d0, -0.745052720131169d0, -0.745052720131169d0 /) 489 490 REAL(KIND=dp), DIMENSION(36), PRIVATE :: VPQuad36 = & 491 (/ 0.000000000000000d0, 0.000000000000000d0, 0.108605615857397d1, & 492 -0.108605615857397d1, 0.000000000000000d0, 0.000000000000000d0, & 493 0.658208197042585d0, -0.658208197042585d0, 0.100130060299173d1, & 494 -0.100130060299173d1, 0.100130060299173d1, -0.100130060299173d1, & 495 0.584636168775946d0, -0.584636168775946d0, 0.584636168775946d0, & 496 -0.584636168775946d0, 0.246795612720261d0, -0.246795612720261d0, & 497 0.246795612720261d0, -0.246795612720261d0, 0.304720678579870d0, & 498 -0.304720678579870d0, 0.304720678579870d0, -0.304720678579870d0, & 499 0.900258815287201d0, -0.900258815287201d0, 0.900258815287201d0, & 500 -0.900258815287201d0, 0.745052720131169d0, -0.745052720131169d0, & 501 0.745052720131169d0, -0.745052720131169d0, 0.929866705560780d0, & 502 -0.929866705560780d0, 0.929866705560780d0, -0.929866705560780d0 /) 503 504 REAL(KIND=dp), DIMENSION(36), PRIVATE :: SPQuad36 = & 505 (/ 0.565616969376400d-2, 0.565616969376400d-2, 0.565616969376400d-2, & 506 0.565616969376400d-2, 0.192443967470396d0, 0.192443967470396d0, & 507 0.192443967470396d0, 0.192443967470396d0, 0.516683297977300d-2, & 508 0.516683297977300d-2, 0.516683297977300d-2, 0.516683297977300d-2, & 509 0.200302559622138d0, 0.200302559622138d0, 0.200302559622138d0, & 510 0.200302559622138d0, 0.228125175912536d0, 0.228125175912536d0, & 511 0.228125175912536d0, 0.228125175912536d0, 0.117496926974491d0, & 512 0.117496926974491d0, 0.117496926974491d0, 0.117496926974491d0, & 513 0.117496926974491d0, 0.117496926974491d0, 0.117496926974491d0, & 514 0.117496926974491d0, 0.666557701862050d-1, 0.666557701862050d-1, & 515 0.666557701862050d-1, 0.666557701862050d-1, 0.666557701862050d-1, & 516 0.666557701862050d-1, 0.666557701862050d-1, 0.666557701862050d-1 /) 517!------------------------------------------------------------------------------ 518! Quadrilateral - 45-point rule for complete polynomials of order p<=15 519! Note: Some points lie outside the reference element 520!------------------------------------------------------------------------------ 521 REAL(KIND=dp), DIMENSION(45), PRIVATE :: UPQuad45 = & 522 (/ 0.000000000000000d0, 0.102731435771937d1, -0.102731435771937d1, & 523 0.000000000000000d0, 0.000000000000000d0, 0.856766776147643d0, & 524 -0.856766776147643d0, 0.000000000000000d0, 0.000000000000000d0, & 525 0.327332998189723d0, -0.327332998189723d0, 0.000000000000000d0, & 526 0.000000000000000d0, 0.967223740028505d0, 0.967223740028505d0, & 527 -0.967223740028505d0, -0.967223740028505d0, 0.732168901749711d0, & 528 0.732168901749711d0, -0.732168901749711d0, -0.732168901749711d0, & 529 0.621974427996805d0, 0.621974427996805d0, -0.621974427996805d0, & 530 -0.621974427996805d0, 0.321696694921009d0, 0.321696694921009d0, & 531 -0.321696694921009d0, -0.321696694921009d0, 0.928618480068352d0, & 532 0.928618480068352d0, -0.928618480068352d0, -0.928618480068352d0, & 533 0.455124178121179d0, 0.455124178121179d0, -0.455124178121179d0, & 534 -0.455124178121179d0, 0.960457474887516d0, 0.960457474887516d0, & 535 -0.960457474887516d0, -0.960457474887516d0, 0.809863684081217d0, & 536 0.809863684081217d0, -0.809863684081217d0, -0.809863684081217d0 /) 537 538 REAL(KIND=dp), DIMENSION(45), PRIVATE :: VPQuad45 = & 539 (/ 0.000000000000000d0, 0.000000000000000d0, 0.000000000000000d0, & 540 0.102731435771937d1, -0.102731435771937d1, 0.000000000000000d0, & 541 0.000000000000000d0, 0.856766776147643d0, -0.856766776147643d0, & 542 0.000000000000000d0, 0.000000000000000d0, 0.327332998189723d0, & 543 -0.327332998189723d0, 0.967223740028505d0, -0.967223740028505d0, & 544 0.967223740028505d0, -0.967223740028505d0, 0.732168901749711d0, & 545 -0.732168901749711d0, 0.732168901749711d0, -0.732168901749711d0, & 546 0.321696694921009d0, -0.321696694921009d0, 0.321696694921009d0, & 547 -0.321696694921009d0, 0.621974427996805d0, -0.621974427996805d0, & 548 0.621974427996805d0, -0.621974427996805d0, 0.455124178121179d0, & 549 -0.455124178121179d0, 0.455124178121179d0, -0.455124178121179d0, & 550 0.928618480068352d0, -0.928618480068352d0, 0.928618480068352d0, & 551 -0.928618480068352d0, 0.809863684081217d0, -0.809863684081217d0, & 552 0.809863684081217d0, -0.809863684081217d0, 0.960457474887516d0, & 553 -0.960457474887516d0, 0.960457474887516d0, -0.960457474887516d0 /) 554 555 REAL(KIND=dp), DIMENSION(45), PRIVATE :: SPQuad45 = & 556 (/ -0.176897982720700d-2, 0.128167266175120d-1, 0.128167266175120d-1, & 557 0.128167266175120d-1, 0.128167266175120d-1, 0.119897873101347d0, & 558 0.119897873101347d0, 0.119897873101347d0, 0.119897873101347d0, & 559 0.210885452208801d0, 0.210885452208801d0, 0.210885452208801d0, & 560 0.210885452208801d0, 0.639272012821500d-2, 0.639272012821500d-2, & 561 0.639272012821500d-2, 0.639272012821500d-2, 0.104415680788580d0, & 562 0.104415680788580d0, 0.104415680788580d0, 0.104415680788580d0, & 563 0.168053047203816d0, 0.168053047203816d0, 0.168053047203816d0, & 564 0.168053047203816d0, 0.168053047203816d0, 0.168053047203816d0, & 565 0.168053047203816d0, 0.168053047203816d0, 0.761696944522940d-1, & 566 0.761696944522940d-1, 0.761696944522940d-1, 0.761696944522940d-1, & 567 0.761696944522940d-1, 0.761696944522940d-1, 0.761696944522940d-1, & 568 0.761696944522940d-1, 0.287941544000640d-1, 0.287941544000640d-1, & 569 0.287941544000640d-1, 0.287941544000640d-1, 0.287941544000640d-1, & 570 0.287941544000640d-1, 0.287941544000640d-1, 0.287941544000640d-1 /) 571!------------------------------------------------------------------------------ 572! Quadrilateral - 60-point rule for complete polynomials of order p<=17 573!------------------------------------------------------------------------------ 574 REAL(KIND=dp), DIMENSION(60), PRIVATE :: UPQuad60 = & 575 (/ 0.989353074512600d0, -0.989353074512600d0, 0.000000000000000d0, & 576 0.000000000000000d0, 0.376285207157973d0, -0.376285207157973d0, & 577 0.000000000000000d0, 0.000000000000000d0, 0.978848279262233d0, & 578 0.978848279262233d0, -0.978848279262233d0, -0.978848279262233d0, & 579 0.885794729164116d0, 0.885794729164116d0, -0.885794729164116d0, & 580 -0.885794729164116d0, 0.171756123838348d0, 0.171756123838348d0, & 581 -0.171756123838348d0, -0.171756123838348d0, 0.590499273806002d0, & 582 0.590499273806002d0, -0.590499273806002d0, -0.590499273806002d0, & 583 0.319505036634574d0, 0.319505036634574d0, -0.319505036634574d0, & 584 -0.319505036634574d0, 0.799079131916863d0, 0.799079131916863d0, & 585 -0.799079131916863d0, -0.799079131916863d0, 0.597972451929457d0, & 586 0.597972451929457d0, -0.597972451929457d0, -0.597972451929457d0, & 587 0.803743962958745d0, 0.803743962958745d0, -0.803743962958745d0, & 588 -0.803743962958745d0, 0.583444817765510d-1, 0.583444817765510d-1, & 589 -0.583444817765510d-1, -0.583444817765510d-1, 0.936506276127495d0, & 590 0.936506276127495d0, -0.936506276127495d0, -0.936506276127495d0, & 591 0.347386316166203d0, 0.347386316166203d0, -0.347386316166203d0, & 592 -0.347386316166203d0, 0.981321179805452d0, 0.981321179805452d0, & 593 -0.981321179805452d0, -0.981321179805452d0, 0.706000287798646d0, & 594 0.706000287798646d0, -0.706000287798646d0, -0.706000287798646d0 /) 595 596 REAL(KIND=dp), DIMENSION(60), PRIVATE :: VPQuad60 = & 597 (/0.000000000000000d0, 0.000000000000000d0, 0.989353074512600d0, & 598 -0.989353074512600d0, 0.000000000000000d0, 0.000000000000000d0, & 599 0.376285207157973d0, -0.376285207157973d0, 0.978848279262233d0, & 600 -0.978848279262233d0, 0.978848279262233d0, -0.978848279262233d0, & 601 0.885794729164116d0, -0.885794729164116d0, 0.885794729164116d0, & 602 -0.885794729164116d0, 0.171756123838348d0, -0.171756123838348d0, & 603 0.171756123838348d0, -0.171756123838348d0, 0.319505036634574d0, & 604 -0.319505036634574d0, 0.319505036634574d0, -0.319505036634574d0, & 605 0.590499273806002d0, -0.590499273806002d0, 0.590499273806002d0, & 606 -0.590499273806002d0, 0.597972451929457d0, -0.597972451929457d0, & 607 0.597972451929457d0, -0.597972451929457d0, 0.799079131916863d0, & 608 -0.799079131916863d0, 0.799079131916863d0, -0.799079131916863d0, & 609 0.583444817765510d-1, -0.583444817765510d-1, 0.583444817765510d-1, & 610 -0.583444817765510d-1, 0.803743962958745d0, -0.803743962958745d0, & 611 0.803743962958745d0, -0.803743962958745d0, 0.347386316166203d0, & 612 -0.347386316166203d0, 0.347386316166203d0, -0.347386316166203d0, & 613 0.936506276127495d0, -0.936506276127495d0, 0.936506276127495d0, & 614 -0.936506276127495d0, 0.706000287798646d0, -0.706000287798646d0, & 615 0.706000287798646d0, -0.706000287798646d0, 0.981321179805452d0, & 616 -0.981321179805452d0, 0.981321179805452d0, -0.981321179805452d0 /) 617 618 REAL(KIND=dp), DIMENSION(60), PRIVATE :: SPQuad60 = & 619 (/ 0.206149159199910d-1, 0.206149159199910d-1, 0.206149159199910d-1, & 620 0.206149159199910d-1, 0.128025716179910d0, 0.128025716179910d0, & 621 0.128025716179910d0, 0.128025716179910d0, 0.551173953403200d-2, & 622 0.551173953403200d-2, 0.551173953403200d-2, 0.551173953403200d-2, & 623 0.392077124571420d-1, 0.392077124571420d-1, 0.392077124571420d-1, & 624 0.392077124571420d-1, 0.763969450798630d-1, 0.763969450798630d-1, & 625 0.763969450798630d-1, 0.763969450798630d-1, 0.141513729949972d0, & 626 0.141513729949972d0, 0.141513729949972d0, 0.141513729949972d0, & 627 0.141513729949972d0, 0.141513729949972d0, 0.141513729949972d0, & 628 0.141513729949972d0, 0.839032793637980d-1, 0.839032793637980d-1, & 629 0.839032793637980d-1, 0.839032793637980d-1, 0.839032793637980d-1, & 630 0.839032793637980d-1, 0.839032793637980d-1, 0.839032793637980d-1, & 631 0.603941636496850d-1, 0.603941636496850d-1, 0.603941636496850d-1, & 632 0.603941636496850d-1, 0.603941636496850d-1, 0.603941636496850d-1, & 633 0.603941636496850d-1, 0.603941636496850d-1, 0.573877529692130d-1, & 634 0.573877529692130d-1, 0.573877529692130d-1, 0.573877529692130d-1, & 635 0.573877529692130d-1, 0.573877529692130d-1, 0.573877529692130d-1, & 636 0.573877529692130d-1, 0.219225594818640d-1, 0.219225594818640d-1, & 637 0.219225594818640d-1, 0.219225594818640d-1, 0.219225594818640d-1, & 638 0.219225594818640d-1, 0.219225594818640d-1, 0.219225594818640d-1 /) 639 640 641!------------------------------------------------------------------------------ 642! These rules have been reproduced from the paper 643! Ethan J. Kubatko, Benjamin A. Yeager, Ashley L. Magg and I. Dolezel: 644! "New computationally efficient quadrature formulas for triangular 645! prism elements.", Computers & Fluids 73 (2013) 187–201. 646!------------------------------------------------------------------------------ 647 648!------------------------------------------------------------------------------ 649! Wedge - 4 point rule, 2nd order polynoms 650!------------------------------------------------------------------------------ 651 REAL(KIND=dp), DIMENSION(4), PRIVATE :: SWedge4P = [ & 652 1.000000000000000d0, & 653 1.000000000000000d0, & 654 1.000000000000000d0, & 655 1.000000000000000d0] 656 REAL(KIND=dp), DIMENSION(4), PRIVATE :: UWedge4P = [ & 657 0.483163247594393d0, & 658 -0.605498860309242d0, & 659 -0.605498860309242d0, & 660 -0.605498860309242d0] 661 REAL(KIND=dp), DIMENSION(4), PRIVATE :: VWedge4P = [ & 662 -0.741581623797196d0, & 663 0.469416096821288d0, & 664 -0.530583903178712d0, & 665 -0.530583903178712d0] 666 REAL(KIND=dp), DIMENSION(4), PRIVATE :: WWedge4P = [ & 667 0.000000000000000d0, & 668 0.000000000000000d0, & 669 0.816496580927726d0, & 670 -0.816496580927726d0 ] 671 672!------------------------------------------------------------------------------ 673! Wedge - n=5, p=2 674!------------------------------------------------------------------------------ 675 REAL(KIND=dp), DIMENSION(5), PRIVATE :: SWedge5P = [ & 676 0.333333333333333d0, & 677 0.333333333333333d0, & 678 0.333333333333333d0, & 679 1.500000000000000d0, & 680 1.500000000000000d0 ] 681 REAL(KIND=dp), DIMENSION(5), PRIVATE :: UWedge5P = [ & 682 -1.000000000000000d0, & 683 -1.000000000000000d0, & 684 1.000000000000000d0, & 685 -0.333333333333333d0, & 686 -0.333333333333333d0 ] 687 REAL(KIND=dp), DIMENSION(5), PRIVATE :: VWedge5P = [ & 688 -1.000000000000000d0, & 689 1.000000000000000d0, & 690 -1.000000000000000d0, & 691 -0.333333333333333d0, & 692 -0.333333333333333d0 ] 693 REAL(KIND=dp), DIMENSION(5), PRIVATE :: WWedge5P = [ & 694 0.000000000000000d0, & 695 0.000000000000000d0, & 696 0.000000000000000d0, & 697 0.666666666666667d0, & 698 -0.666666666666667d0 ] 699 700!------------------------------------------------------------------------------ 701! Wedge - n=6, p=3 702!------------------------------------------------------------------------------ 703 REAL(KIND=dp), DIMENSION(6), PRIVATE :: SWedge6P = [ & 704 0.742534890852309d0, & 705 0.375143463443327d0, & 706 0.495419047908462d0, & 707 0.523999970843238d0, & 708 0.980905839025611d0, & 709 0.881996787927053d0 ] 710 REAL(KIND=dp), DIMENSION(6), PRIVATE :: UWedge6P = [ & 711 0.240692796349625d0, & 712 -0.968326281451138d0, & 713 0.467917833640195d0, & 714 -0.786144119530819d0, & 715 -0.484844897886675d0, & 716 -0.559053711932125d0 ] 717 REAL(KIND=dp), DIMENSION(6), PRIVATE :: VWedge6P = [ & 718 -0.771991660873412d0, & 719 -0.568046512457875d0, & 720 -0.549342790168347d0, & 721 0.362655041695561d0, & 722 -0.707931130717342d0, & 723 0.260243325246813d0 ] 724 REAL(KIND=dp), DIMENSION(6), PRIVATE :: WWedge6P = [ & 725 0.614747128207527d0, & 726 0.676689529541421d0, & 727 -0.599905857322635d0, & 728 -0.677609795694786d0, & 729 -0.502482717716373d0, & 730 0.493010512161538d0 ] 731 732!------------------------------------------------------------------------------ 733! Wedge - n=7, p=3 734!------------------------------------------------------------------------------ 735 REAL(KIND=dp), DIMENSION(7), PRIVATE :: SWedge7P = [ & 736 -2.250000000000000d0, & 737 1.041666666666667d0, & 738 1.041666666666667d0, & 739 1.041666666666667d0, & 740 1.041666666666667d0, & 741 1.041666666666667d0, & 742 1.041666666666667d0 ] 743 REAL(KIND=dp), DIMENSION(7), PRIVATE :: UWedge7P = [ & 744 -0.333333333333333d0, & 745 -0.600000000000000d0, & 746 -0.600000000000000d0, & 747 0.200000000000000d0, & 748 -0.600000000000000d0, & 749 -0.600000000000000d0, & 750 0.200000000000000d0 ] 751 REAL(KIND=dp), DIMENSION(7), PRIVATE :: VWedge7P = [ & 752 -0.333333333333333d0, & 753 -0.600000000000000d0, & 754 0.200000000000000d0, & 755 -0.600000000000000d0, & 756 -0.600000000000000d0, & 757 0.200000000000000d0, & 758 -0.600000000000000d0 ] 759 REAL(KIND=dp), DIMENSION(7), PRIVATE :: WWedge7P = [ & 760 0.000000000000000d0, & 761 0.461880215351701d0, & 762 0.461880215351701d0, & 763 0.461880215351701d0, & 764 -0.461880215351701d0, & 765 -0.461880215351701d0, & 766 -0.461880215351701d0 ] 767 768!------------------------------------------------------------------------------ 769! Wedge - n=10, p=4 770!------------------------------------------------------------------------------ 771 REAL(KIND=dp), DIMENSION(10), PRIVATE :: SWedge10P = [ & 772 0.111155943811228d0, & 773 0.309060899887509d0, & 774 0.516646862442958d0, & 775 0.567975205132714d0, & 776 0.382742555939017d0, & 777 0.355960928492268d0, & 778 0.108183228294342d0, & 779 0.126355242780924d0, & 780 0.587370828592853d0, & 781 0.934548304626188d0 ] 782 REAL(KIND=dp), DIMENSION(10), PRIVATE :: UWedge10P = [ & 783 0.812075900047562d0, & 784 -0.792166223585545d0, & 785 -0.756726179789306d0, & 786 -0.552495167978340d0, & 787 -0.357230019521233d0, & 788 -0.987225392999058d0, & 789 -0.816603728785918d0, & 790 0.423489172633859d0, & 791 0.363041084609230d0, & 792 -0.175780343149613d0 ] 793 REAL(KIND=dp), DIMENSION(10), PRIVATE :: VWedge10P = [ & 794 -0.986242751499303d0, & 795 0.687201105597868d0, & 796 -0.731311840596107d0, & 797 0.015073398439985d0, & 798 0.126888850505978d0, & 799 0.082647545710800d0, & 800 -0.915066171481315d0, & 801 -1.112801167237130d0, & 802 -0.499011410082669d0, & 803 -0.654971142379686d0 ] 804 REAL(KIND=dp), DIMENSION(10), PRIVATE :: WWedge10P = [ & 805 0.850716248413834d0, & 806 -0.115214772515700d0, & 807 -0.451491675441927d0, & 808 -0.824457000064439d0, & 809 0.855349689995606d0, & 810 0.452976444667786d0, & 811 0.997939285245240d0, & 812 -0.963298774205756d0, & 813 -0.299892769705443d0, & 814 0.367947041936472d0 ] 815 816!------------------------------------------------------------------------------ 817! Wedge - n=11, p=4 818!------------------------------------------------------------------------------ 819 REAL(KIND=dp), DIMENSION(11), PRIVATE :: SWedge11P = [ & 820 0.545658450421913d0, & 821 0.545658450421913d0, & 822 0.545658450421913d0, & 823 0.431647899262139d0, & 824 0.249954808368331d0, & 825 0.249954808368331d0, & 826 0.249954808368331d0, & 827 0.431647899262139d0, & 828 0.249954808368331d0, & 829 0.249954808368331d0, & 830 0.249954808368331d0 ] 831 REAL(KIND=dp), DIMENSION(11), PRIVATE :: UWedge11P = [ & 832 -0.062688380276010d0, & 833 -0.062688380276010d0, & 834 -0.874623239447980d0, & 835 -0.333333333333333d0, & 836 -0.798519188402179d0, & 837 -0.798519188402179d0, & 838 0.597038376804357d0, & 839 -0.333333333333333d0, & 840 -0.798519188402179d0, & 841 -0.798519188402179d0, & 842 0.597038376804357d0 ] 843 REAL(KIND=dp), DIMENSION(11), PRIVATE :: VWedge11P = [ & 844 -0.062688380276010d0, & 845 -0.874623239447980d0, & 846 -0.062688380276010d0, & 847 -0.333333333333333d0, & 848 -0.798519188402179d0, & 849 0.597038376804357d0, & 850 -0.798519188402179d0, & 851 -0.333333333333333d0, & 852 -0.798519188402179d0, & 853 0.597038376804357d0, & 854 -0.798519188402179d0 ] 855 REAL(KIND=dp), DIMENSION(11), PRIVATE :: WWedge11P = [ & 856 0.000000000000000d0, & 857 0.000000000000000d0, & 858 0.000000000000000d0, & 859 0.866861974009030d0, & 860 0.675639823682265d0, & 861 0.675639823682265d0, & 862 0.675639823682265d0, & 863 -0.866861974009030d0, & 864 -0.675639823682265d0, & 865 -0.675639823682265d0, & 866 -0.675639823682265d0 ] 867 868!------------------------------------------------------------------------------ 869! Wedge - n=14, p=5 870!------------------------------------------------------------------------------ 871 REAL(KIND=dp), DIMENSION(14), PRIVATE :: SWedge14P = [ & 872 0.087576186678730d0, & 873 0.229447629454892d0, & 874 0.229447629454891d0, & 875 0.833056798542985d0, & 876 0.166443428304729d0, & 877 0.166443428304729d0, & 878 0.376993270712316d0, & 879 0.170410864470884d0, & 880 0.298194157223163d0, & 881 0.298194157223162d0, & 882 0.376993270712316d0, & 883 0.170410864470884d0, & 884 0.298194157223163d0, & 885 0.298194157223162d0 ] 886 REAL(KIND=dp), DIMENSION(14), PRIVATE :: UWedge14P = [ & 887 -0.955901336867574d0, & 888 -0.051621305926029d0, & 889 -1.017063354038640d0, & 890 -0.297388746460523d0, & 891 0.774849157169622d0, & 892 -0.849021640062097d0, & 893 -0.665292008657551d0, & 894 -0.012171148087346d0, & 895 -0.734122164680096d0, & 896 0.334122164680066d0, & 897 -0.665292008657551d0, & 898 -0.012171148087346d0, & 899 -0.734122164680096d0, & 900 0.334122164680066d0 ] 901 REAL(KIND=dp), DIMENSION(14), PRIVATE :: VWedge14P = [ & 902 -0.955901336867577d0, & 903 -1.017063354038640d0, & 904 -0.051621305926032d0, & 905 -0.297388746460521d0, & 906 -0.849021640062096d0, & 907 0.774849157169620d0, & 908 -0.665292008657551d0, & 909 -0.012171148087349d0, & 910 0.334122164680066d0, & 911 -0.734122164680096d0, & 912 -0.665292008657551d0, & 913 -0.012171148087349d0, & 914 0.334122164680066d0, & 915 -0.734122164680096d0 ] 916 REAL(KIND=dp), DIMENSION(14), PRIVATE :: WWedge14P = [ & 917 0.000000000000286d0, & 918 0.000000000000038d0, & 919 0.000000000000038d0, & 920 0.000000000000008d0, & 921 0.000000000000080d0, & 922 0.000000000000080d0, & 923 0.756615409654429d0, & 924 0.600149379161583d0, & 925 0.808115770496521d0, & 926 0.808115770496521d0, & 927 -0.756615409654429d0, & 928 -0.600149379161583d0, & 929 -0.808115770496521d0, & 930 -0.808115770496521d0 ] 931 932!------------------------------------------------------------------------------ 933! Wedge - n=15, p=5 934!------------------------------------------------------------------------------ 935 REAL(KIND=dp), DIMENSION(15), PRIVATE :: SWedge15P = [ & 936 0.213895020288765d0, & 937 0.141917375616806d0, & 938 0.295568859378071d0, & 939 0.256991945593379d0, & 940 0.122121979248381d0, & 941 0.175194917962627d0, & 942 0.284969106392719d0, & 943 0.323336131783334d0, & 944 0.159056110329063d0, & 945 0.748067388709644d0, & 946 0.280551223607231d0, & 947 0.147734016552639d0, & 948 0.259874920383688d0, & 949 0.235144061421191d0, & 950 0.355576942732463d0 ] 951 REAL(KIND=dp), DIMENSION(15), PRIVATE :: UWedge15P = [ & 952 -0.820754415297359d0, & 953 0.611831616907812d0, & 954 -0.951379065092975d0, & 955 0.200535109198601d0, & 956 -0.909622841605196d0, & 957 0.411514133287729d0, & 958 -0.127534496411879d0, & 959 -0.555217727817199d0, & 960 0.706942532529193d0, & 961 -0.278092963133809d0, & 962 -0.057824844208300d0, & 963 -0.043308436222116d0, & 964 -0.774478920734726d0, & 965 -0.765638443571458d0, & 966 -0.732830649614460d0 ] 967 REAL(KIND=dp), DIMENSION(15), PRIVATE :: VWedge15P = [ & 968 0.701020947925133d0, & 969 -0.869995576950389d0, & 970 -0.087091980055873d0, & 971 -0.783721735474016d0, & 972 -0.890218158063352d0, & 973 -0.725374126531787d0, & 974 -0.953467283619037d0, & 975 -0.530472194607007d0, & 976 -0.782248553944540d0, & 977 -0.291936530517119d0, & 978 -0.056757587543798d0, & 979 0.012890722780611d0, & 980 0.476188541042454d0, & 981 0.177195164202219d0, & 982 -0.737447982744191d0 ] 983 REAL(KIND=dp), DIMENSION(15), PRIVATE :: WWedge15P = [ & 984 -0.300763696502910d0, & 985 0.348546607420888d0, & 986 0.150656042323906d0, & 987 -0.844285153176719d0, & 988 0.477120081549168d0, & 989 0.864653509536562d0, & 990 0.216019762875977d0, & 991 0.873409672725819d0, & 992 -0.390653804976705d0, & 993 -0.126030507204870d0, & 994 0.539907869785112d0, & 995 -0.776314479909204d0, & 996 0.745875967497062d0, & 997 -0.888355356215127d0, & 998 -0.651653242952189d0 ] 999 1000!------------------------------------------------------------------------------ 1001! Wedge - n=16, p=5 1002!------------------------------------------------------------------------------ 1003 REAL(KIND=dp), DIMENSION(16), PRIVATE :: SWedge16P = [ & 1004 0.711455555931488d0, & 1005 0.224710067228267d0, & 1006 0.224710067228267d0, & 1007 0.224710067228267d0, & 1008 0.185661421316158d0, & 1009 0.185661421316158d0, & 1010 0.185661421316158d0, & 1011 0.250074285747794d0, & 1012 0.250074285747794d0, & 1013 0.250074285747794d0, & 1014 0.185661421316158d0, & 1015 0.185661421316158d0, & 1016 0.185661421316158d0, & 1017 0.250074285747794d0, & 1018 0.250074285747794d0, & 1019 0.250074285747794d0 ] 1020 REAL(KIND=dp), DIMENSION(16), PRIVATE :: UWedge16P = [ & 1021 -0.333333333333333d0, & 1022 -0.025400070899509d0, & 1023 -0.025400070899509d0, & 1024 -0.949199858200983d0, & 1025 -0.108803790659256d0, & 1026 -0.108803790659256d0, & 1027 -0.782392418681488d0, & 1028 -0.798282108034583d0, & 1029 -0.798282108034583d0, & 1030 0.596564216069166d0, & 1031 -0.108803790659256d0, & 1032 -0.108803790659256d0, & 1033 -0.782392418681488d0, & 1034 -0.798282108034583d0, & 1035 -0.798282108034583d0, & 1036 0.596564216069166d0 ] 1037 REAL(KIND=dp), DIMENSION(16), PRIVATE :: VWedge16P = [ & 1038 -0.333333333333333d0, & 1039 -0.025400070899509d0, & 1040 -0.949199858200983d0, & 1041 -0.025400070899509d0, & 1042 -0.108803790659256d0, & 1043 -0.782392418681488d0, & 1044 -0.108803790659256d0, & 1045 -0.798282108034583d0, & 1046 0.596564216069166d0, & 1047 -0.798282108034583d0, & 1048 -0.108803790659256d0, & 1049 -0.782392418681488d0, & 1050 -0.108803790659256d0, & 1051 -0.798282108034583d0, & 1052 0.596564216069166d0, & 1053 -0.798282108034583d0 ] 1054 REAL(KIND=dp), DIMENSION(16), PRIVATE :: WWedge16P = [ & 1055 0.000000000000000d0, & 1056 0.000000000000000d0, & 1057 0.000000000000000d0, & 1058 0.000000000000000d0, & 1059 0.871002934865444d0, & 1060 0.871002934865444d0, & 1061 0.871002934865444d0, & 1062 0.570426980705159d0, & 1063 0.570426980705159d0, & 1064 0.570426980705159d0, & 1065 -0.871002934865444d0, & 1066 -0.871002934865444d0, & 1067 -0.871002934865444d0, & 1068 -0.570426980705159d0, & 1069 -0.570426980705159d0, & 1070 -0.570426980705159d0 ] 1071 1072!------------------------------------------------------------------------------ 1073! Wedge - n=24, p=6 1074!------------------------------------------------------------------------------ 1075 REAL(KIND=dp), DIMENSION(24), PRIVATE :: SWedge24P = [ & 1076 0.168480079940386d0, & 1077 0.168480079940386d0, & 1078 0.168480079940386d0, & 1079 0.168480079940386d0, & 1080 0.168480079940386d0, & 1081 0.168480079940386d0, & 1082 0.000079282874851d0, & 1083 0.000079282874851d0, & 1084 0.000079282874851d0, & 1085 0.544286440652304d0, & 1086 0.026293733850586d0, & 1087 0.283472344926041d0, & 1088 0.283472344926041d0, & 1089 0.283472344926041d0, & 1090 0.115195615637235d0, & 1091 0.115195615637235d0, & 1092 0.115195615637235d0, & 1093 0.026293733850586d0, & 1094 0.283472344926041d0, & 1095 0.283472344926041d0, & 1096 0.283472344926041d0, & 1097 0.115195615637235d0, & 1098 0.115195615637235d0, & 1099 0.115195615637235d0 ] 1100 REAL(KIND=dp), DIMENSION(24), PRIVATE :: UWedge24P = [ & 1101 0.513019949700545d0, & 1102 -0.582925557762337d0, & 1103 -0.582925557762337d0, & 1104 0.513019949700545d0, & 1105 -0.930094391938207d0, & 1106 -0.930094391938207d0, & 1107 -1.830988812620400d0, & 1108 -1.830988812620400d0, & 1109 2.661977625240810d0, & 1110 -0.333333333333333d0, & 1111 -0.333333333333333d0, & 1112 -0.098283514203544d0, & 1113 -0.098283514203544d0, & 1114 -0.803432971592912d0, & 1115 -0.812603471654584d0, & 1116 -0.812603471654584d0, & 1117 0.625206943309168d0, & 1118 -0.333333333333333d0, & 1119 -0.098283514203544d0, & 1120 -0.098283514203544d0, & 1121 -0.803432971592912d0, & 1122 -0.812603471654584d0, & 1123 -0.812603471654584d0, & 1124 0.625206943309168d0 ] 1125 REAL(KIND=dp), DIMENSION(24), PRIVATE :: VWedge24P = [ & 1126 -0.930094391938207d0, & 1127 -0.930094391938207d0, & 1128 0.513019949700545d0, & 1129 -0.582925557762337d0, & 1130 -0.582925557762337d0, & 1131 0.513019949700545d0, & 1132 -1.830988812620400d0, & 1133 2.661977625240810d0, & 1134 -1.830988812620400d0, & 1135 -0.333333333333333d0, & 1136 -0.333333333333333d0, & 1137 -0.098283514203544d0, & 1138 -0.803432971592912d0, & 1139 -0.098283514203544d0, & 1140 -0.812603471654584d0, & 1141 0.625206943309168d0, & 1142 -0.812603471654584d0, & 1143 -0.333333333333333d0, & 1144 -0.098283514203544d0, & 1145 -0.803432971592912d0, & 1146 -0.098283514203544d0, & 1147 -0.812603471654584d0, & 1148 0.625206943309168d0, & 1149 -0.812603471654584d0 ] 1150 REAL(KIND=dp), DIMENSION(24), PRIVATE :: WWedge24P = [ & 1151 0.000000000000000d0, & 1152 0.000000000000000d0, & 1153 0.000000000000000d0, & 1154 0.000000000000000d0, & 1155 0.000000000000000d0, & 1156 0.000000000000000d0, & 1157 0.000000000000000d0, & 1158 0.000000000000000d0, & 1159 0.000000000000000d0, & 1160 0.000000000000000d0, & 1161 1.250521622121900d0, & 1162 0.685008566774710d0, & 1163 0.685008566774710d0, & 1164 0.685008566774710d0, & 1165 0.809574716992997d0, & 1166 0.809574716992997d0, & 1167 0.809574716992997d0, & 1168 -1.250521622121900d0, & 1169 -0.685008566774710d0, & 1170 -0.685008566774710d0, & 1171 -0.685008566774710d0, & 1172 -0.809574716992997d0, & 1173 -0.809574716992997d0, & 1174 -0.809574716992997d0 ] 1175 1176! continue wedge rules here 1177!------------------------------------------------------------------------------ 1178! Wedge - n=27, p=6 1179!------------------------------------------------------------------------------ 1180! REAL(KIND=dp), DIMENSION(27), PRIVATE :: SWedge27P = [ & 1181! REAL(KIND=dp), DIMENSION(27), PRIVATE :: UWedge27P = [ & 1182! REAL(KIND=dp), DIMENSION(27), PRIVATE :: VWedge27P = [ & 1183! REAL(KIND=dp), DIMENSION(27), PRIVATE :: WWedge27P = [ & 1184 1185 1186 1187 CONTAINS 1188 1189 1190!------------------------------------------------------------------------------ 1191!> Subroutine to compute Fejer integration points by FFT. 1192!------------------------------------------------------------------------------ 1193 SUBROUTINE ComputeFejerPoints1D( Points,Weights,n ) 1194!------------------------------------------------------------------------------ 1195 INTEGER :: n 1196 REAL(KIND=dp) :: Points(n),Weights(n) 1197!------------------------------------------------------------------------------ 1198 INTEGER :: i,j,k,l,m 1199 TYPE(C_FFTCMPLX) :: W(n+1) 1200 REAL(KIND=dp) :: arr((n+1)/2+1), V(n+2) 1201 1202 DO i=1,(n+1)/2 1203 Points(i) = COS(i*PI/(n+1._dp)) 1204 Points(n-i+1) = -Points(i) 1205 END DO 1206 1207 k = 0 1208 DO i=1,n+1,2 1209 k = k + 1 1210 arr(k) = i 1211 END DO 1212 1213 V = 0._dp 1214 DO i=1,k 1215 V(i) = 2._dp/(arr(i)*(arr(i)-2)) 1216 END DO 1217 V(k+1) = 1._dp/arr(k) 1218 1219 DO i=1,n+1 1220 W(i) % rval = -(V(i)+V(n-i+3)) 1221 END DO 1222 CALL frfftb(n+1,W,V) 1223 1224 Weights(1:n) = V(2:n+1)/(n+1) 1225 DO i=1,n/2 1226 Weights(i) = (Weights(i) + Weights(n-i+1))/2 1227 Weights(n-i+1) = Weights(i) 1228 END DO 1229 Weights(1:n) = 2 * Weights(1:n) / SUM(Weights(1:n)) 1230!------------------------------------------------------------------------------ 1231 END SUBROUTINE ComputeFejerPoints1D 1232!------------------------------------------------------------------------------ 1233 1234 1235!------------------------------------------------------------------------------ 1236!> Subroutine to compute Gaussian integration points and weights in [-1,1] 1237!> as roots of Legendre polynomials. 1238!------------------------------------------------------------------------------ 1239 SUBROUTINE ComputeGaussPoints1D( Points,Weights,n ) 1240!------------------------------------------------------------------------------ 1241 INTEGER :: n 1242 REAL(KIND=dp) :: Points(n),Weights(n) 1243!------------------------------------------------------------------------------ 1244 REAL(KIND=dp) :: A(n/2,n/2),s,x,Work(8*n) 1245 COMPLEX(KIND=dp) :: Eigs(n/2) 1246 REAL(KIND=dp) :: P(n+1),Q(n),P0(n),P1(n+1) 1247 INTEGER :: i,j,k,np,info 1248!------------------------------------------------------------------------------ 1249! One point is trivial 1250!------------------------------------------------------------------------------ 1251 IF ( n <= 1 ) THEN 1252 Points(1) = 0.0d0 1253 Weights(1) = 2.0d0 1254 RETURN 1255 END IF 1256!------------------------------------------------------------------------------ 1257! Compute coefficients of n:th Legendre polynomial from the recurrence: 1258! 1259! (i+1)P_{i+1}(x) = (2i+1)*x*P_i(x) - i*P_{i-1}(x), P_{0} = 1; P_{1} = x; 1260! 1261! CAVEAT: Computed coefficients inaccurate for n > ~15 1262!------------------------------------------------------------------------------ 1263 P = 0 1264 P0(1) = 1 1265 P1(1) = 1 1266 P1(2) = 0 1267 1268 DO i=1,n-1 1269 P(1:i+1) = (2*i+1) * P1(1:i+1) / (i+1) 1270 P(3:i+2) = P(3:i+2) - i*P0(1:i) / (i+1) 1271 P0(1:i+1) = P1(1:i+1); P1(1:i+2) = P(1:i+2) 1272 END DO 1273!------------------------------------------------------------------------------ 1274! Odd n implicates zero as one of the roots... 1275!------------------------------------------------------------------------------ 1276 np = n - MOD(n,2) 1277!------------------------------------------------------------------------------ 1278! Variable substitution: y=x^2 1279!------------------------------------------------------------------------------ 1280 np = np / 2 1281 DO i=1,np+1 1282 P(i) = P(2*i-1) 1283 END DO 1284!------------------------------------------------------------------------------ 1285! Solve the roots of the polynomial by forming a matrix whose characteristic 1286! polynomial is the n:th Legendre polynomial and solving for the eigenvalues. 1287! Dunno if this is a very good method.... 1288!------------------------------------------------------------------------------ 1289 A=0 1290 DO i=1,np-1 1291 A(i,i+1) = 1 1292 END DO 1293 1294 DO i=1,np 1295 A(np,i) = -P(np+2-i) / P(1) 1296 END DO 1297 CALL DGEEV( 'N','N',np,A,n/2,Points,P0,Work,1,Work,1,Work,8*n,info ) 1298! CALL EigenValues( A, np, Eigs ) 1299! Points(1:np) = REAL( Eigs(1:np) ) 1300 1301!------------------------------------------------------------------------------ 1302! Backsubstitute from y=x^2 1303!------------------------------------------------------------------------------ 1304 Q(1:np+1) = P(1:np+1) 1305 P = 0 1306 DO i=1,np+1 1307 P(2*i-1) = Q(i) 1308 END DO 1309 1310 Q(1:np) = Points(1:np) 1311 DO i=1,np 1312 Points(2*i-1) = +SQRT( Q(i) ) 1313 Points(2*i) = -SQRT( Q(i) ) 1314 END DO 1315 IF ( MOD(n,2) == 1 ) Points(n) = 0.0d0 1316 1317 CALL DerivPoly( n,Q,P ) 1318 CALL RefineRoots( n,P,Q,Points ) 1319!------------------------------------------------------------------------------ 1320! Check for roots 1321!------------------------------------------------------------------------------ 1322 DO i=1,n 1323 s = EvalPoly( n,P,Points(i) ) 1324 IF ( ABS(s) > 1.0d-12 ) THEN 1325 CALL Warn( 'ComputeGaussPoints1D', & 1326 '-------------------------------------------------------' ) 1327 CALL Warn( 'ComputeGaussPoints1D', 'Computed integration point' ) 1328 CALL Warn( 'ComputeGaussPoints1D', 'seems to be inaccurate: ' ) 1329 WRITE( Message, * ) 'Points req.: ',n 1330 CALL Warn( 'ComputeGaussPoints1D', Message ) 1331 WRITE( Message, * ) 'Residual: ',s 1332 CALL Warn( 'ComputeGaussPoints1D', Message ) 1333 WRITE( Message, * ) 'Point: +-', SQRT(Points(i)) 1334 CALL Warn( 'ComputeGaussPoints1D', Message ) 1335 CALL Warn( 'ComputeGaussPoints1D', & 1336 '-------------------------------------------------------' ) 1337 END IF 1338 END DO 1339!------------------------------------------------------------------------------ 1340! Finally, the integration weights equal to 1341! 1342! W_i = 2/( (1-x_i^2)*Q(x_i)^2 ), x_i is the i:th root, and Q(x) = dP(x) / dx 1343!------------------------------------------------------------------------------ 1344 CALL DerivPoly( n,Q,P ) 1345 1346 DO i=1,n 1347 s = EvalPoly( n-1,Q,Points(i) ) 1348 Weights(i) = 2 / ((1-Points(i)**2)*s**2); 1349 END DO 1350!------------------------------------------------------------------------------ 1351! ... make really sure the weights add up: 1352!------------------------------------------------------------------------------ 1353 Weights(1:n) = 2 * Weights(1:n) / SUM(Weights(1:n)) 1354 1355CONTAINS 1356 1357!-------------------------------------------------------------------------- 1358 1359 FUNCTION EvalPoly( n,P,x ) RESULT(s) 1360 INTEGER :: i,n 1361 REAL(KIND=dp) :: P(n+1),x,s 1362 1363 s = 0.0d0 1364 DO i=1,n+1 1365 s = s * x + P(i) 1366 END DO 1367 END FUNCTION EvalPoly 1368 1369!-------------------------------------------------------------------------- 1370 1371 SUBROUTINE DerivPoly( n,Q,P ) 1372 INTEGER :: i,n 1373 REAL(KIND=dp) :: Q(n),P(n+1) 1374 1375 DO i=1,n 1376 Q(i) = P(i)*(n-i+1) 1377 END DO 1378 END SUBROUTINE DerivPoly 1379 1380!-------------------------------------------------------------------------- 1381 1382 SUBROUTINE RefineRoots( n,P,Q,Points ) 1383 INTEGER :: i,j,n 1384 REAL(KIND=dp) :: P(n+1),Q(n),Points(n) 1385 1386 REAL(KIND=dp) :: x,s 1387 INTEGER, PARAMETER :: MaxIter = 100 1388 1389 DO i=1,n 1390 x = Points(i) 1391 DO j=1,MaxIter 1392 s = EvalPoly(n,P,x) / EvalPoly(n-1,Q,x) 1393 x = x - s 1394 IF ( ABS(s) <= ABS(x)*EPSILON(s) ) EXIT 1395 END DO 1396 IF ( ABS(EvalPoly(n,P,x))<ABS(EvalPoly(n,P,Points(i))) ) THEN 1397 IF ( ABS(x-Points(i))<1.0d-8 ) Points(i) = x 1398 END IF 1399 END DO 1400 END SUBROUTINE RefineRoots 1401 1402!-------------------------------------------------------------------------- 1403 END SUBROUTINE ComputeGaussPoints1D 1404!-------------------------------------------------------------------------- 1405 1406!-------------------------------------------------------------------------- 1407 FUNCTION GaussPointsInitialized() RESULT(g) 1408!-------------------------------------------------------------------------- 1409 LOGICAL :: g 1410!-------------------------------------------------------------------------- 1411 g=Ginit 1412!-------------------------------------------------------------------------- 1413 END FUNCTION GaussPointsInitialized 1414!-------------------------------------------------------------------------- 1415 1416! !------------------------------------------------------------------------------ 1417! SUBROUTINE GaussPointsInit 1418! !------------------------------------------------------------------------------ 1419! INTEGER :: i,n,istat 1420! INTEGER :: nthreads, thread, omp_get_max_threads, omp_get_thread_num 1421! 1422! IF ( .NOT. ALLOCATED(IntegStuff) ) THEN 1423! !$omp barrier 1424! !$omp critical 1425! IF ( .NOT. GInit ) THEN 1426! GInit = .TRUE. 1427! nthreads = 1 1428! !$ nthreads = omp_get_max_threads() 1429! ALLOCATE( IntegStuff(nthreads) ) 1430! DO i=1,nthreads 1431! IntegStuff(i) % u => NULL() 1432! IntegStuff(i) % v => NULL() 1433! IntegStuff(i) % w => NULL() 1434! IntegStuff(i) % s => NULL() 1435! END DO 1436! DO n=1,MAXN 1437! CALL ComputeGaussPoints1D( Points(1:n,n),Weights(1:n,n),n ) 1438! END DO 1439! END IF 1440! !$omp end critical 1441! END IF 1442! 1443! thread = 1 1444! !$ thread = omp_get_thread_num()+1 1445! ALLOCATE( IntegStuff(thread) % u(MAX_INTEGRATION_POINTS), & 1446! IntegStuff(thread) % v(MAX_INTEGRATION_POINTS), & 1447! IntegStuff(thread) % w(MAX_INTEGRATION_POINTS), & 1448! IntegStuff(thread) % s(MAX_INTEGRATION_POINTS), STAT=istat ) 1449! 1450! IF ( istat /= 0 ) THEN 1451! CALL Fatal( 'GaussPointsInit', 'Memory allocation error.' ) 1452! STOP 1453! END IF 1454! !------------------------------------------------------------------------------ 1455! ! END SUBROUTINE GaussPointsInit 1456! !------------------------------------------------------------------------------ 1457 1458!------------------------------------------------------------------------------ 1459 SUBROUTINE GaussPointsInit 1460!------------------------------------------------------------------------------ 1461 INTEGER :: i,n,istat 1462 1463 IF ( .NOT. GInit ) THEN 1464 DO n=1,MAXN 1465 CALL ComputeGaussPoints1D( Points(1:n,n),Weights(1:n,n),n ) 1466 END DO 1467 GInit = .TRUE. 1468 END IF 1469 1470 ALLOCATE( IntegStuff % u(MAX_INTEGRATION_POINTS), & 1471 IntegStuff % v(MAX_INTEGRATION_POINTS), & 1472 IntegStuff % w(MAX_INTEGRATION_POINTS), & 1473 IntegStuff % s(MAX_INTEGRATION_POINTS), STAT=istat ) 1474 IntegStuff % u = 0._dp 1475 IntegStuff % v = 0._dp 1476 IntegStuff % w = 0._dp 1477 IntegStuff % s = 0._dp 1478 IF ( istat /= 0 ) THEN 1479 CALL Fatal( 'GaussPointsInit', 'Memory allocation error.' ) 1480 END IF 1481!------------------------------------------------------------------------------ 1482 END SUBROUTINE GaussPointsInit 1483!------------------------------------------------------------------------------ 1484 1485 1486!------------------------------------------------------------------------------ 1487 FUNCTION GaussPoints0D( n ) RESULT(p) 1488!------------------------------------------------------------------------------ 1489 INTEGER :: n 1490 TYPE(GaussIntegrationPoints_t), POINTER :: p 1491! INTEGER :: thread, omp_get_thread_num 1492 1493 IF ( .NOT. GInit ) CALL GaussPointsInit 1494! thread = 1 1495! !$ thread = omp_get_thread_num()+1 1496! p => IntegStuff(thread) 1497 p => IntegStuff 1498 p % n = 1 1499 p % u(1) = 0 1500 p % v(1) = 0 1501 p % w(1) = 0 1502 p % s(1) = 1 1503!------------------------------------------------------------------------------ 1504 END FUNCTION GaussPoints0D 1505!------------------------------------------------------------------------------ 1506 1507 1508!------------------------------------------------------------------------------ 1509!> Return Gaussian integration points for 1D line element 1510!------------------------------------------------------------------------------ 1511 FUNCTION GaussPoints1D( n ) RESULT(p) 1512!------------------------------------------------------------------------------ 1513 INTEGER :: n !< number of points in the requested rule 1514 TYPE(GaussIntegrationPoints_t), POINTER :: p 1515!------------------------------------------------------------------------------ 1516! INTEGER :: thread, omp_get_thread_num 1517 1518 IF ( .NOT. GInit ) CALL GaussPointsInit 1519! thread = 1 1520! !$ thread = omp_get_thread_num()+1 1521! p => IntegStuff(thread) 1522 p => IntegStuff 1523 IF ( n < 1 .OR. n > MAXN ) THEN 1524 p % n = 0 1525 WRITE( Message, * ) 'Invalid number of points: ',n 1526 CALL Error( 'GaussPoints1D', Message ) 1527 RETURN 1528 END IF 1529 1530 p % n = n 1531 p % u(1:n) = Points(1:n,n) 1532 p % v(1:n) = 0.0d0 1533 p % w(1:n) = 0.0d0 1534 p % s(1:n) = Weights(1:n,n) 1535!------------------------------------------------------------------------------ 1536 END FUNCTION GaussPoints1D 1537!------------------------------------------------------------------------------ 1538 1539!------------------------------------------------------------------------------ 1540 FUNCTION GaussPointsPTriangle(n) RESULT(p) 1541!------------------------------------------------------------------------------ 1542 INTEGER :: i,n 1543 TYPE(GaussIntegrationPoints_t), POINTER :: p 1544 REAL (KIND=dp) :: uq, vq, sq 1545! INTEGER :: thread, omp_get_thread_num 1546!------------------------------------------------------------------------------ 1547 IF ( .NOT. GInit ) CALL GaussPointsInit 1548! thread = 1 1549! !$ thread = omp_get_thread_num()+1 1550! p => IntegStuff(thread) 1551 p => IntegStuff 1552 1553 ! Construct Gauss points for p (barycentric) triangle from 1554 ! Gauss points for quadrilateral 1555 p = GaussPointsQuad( n ) 1556 1557 ! For each point apply mapping from quad to triangle and 1558 ! multiply weight by detJ of mapping 1559!DIR$ IVDEP 1560 DO i=1,p % n 1561 uq = p % u(i) 1562 vq = p % v(i) 1563 sq = p % s(i) 1564 p % u(i) = 1d0/2*(uq-uq*vq) 1565 p % v(i) = SQRT(3d0)/2*(1d0+vq) 1566 p % s(i) = -SQRT(3d0)/4*(-1+vq)*sq 1567 END DO 1568 1569 p % w(1:n) = 0.0d0 1570!------------------------------------------------------------------------------ 1571 END FUNCTION GaussPointsPTriangle 1572!------------------------------------------------------------------------------ 1573 1574!------------------------------------------------------------------------------ 1575!> Return Gaussian integration points for 2D triangle element. Here the 1576!> reference element over which the integration is done can also be the 1577!> equilateral triangle used in the description of p-elements. In that case, 1578!> this routine may return a more economical set of integration points. 1579!------------------------------------------------------------------------------ 1580 FUNCTION GaussPointsTriangle( n, PReferenceElement ) RESULT(p) 1581!------------------------------------------------------------------------------ 1582 INTEGER :: n !< number of points in the requested rule 1583 LOGICAL, OPTIONAL :: PReferenceElement !< used for switching the reference element 1584 TYPE(GaussIntegrationPoints_t), POINTER :: p 1585!------------------------------------------------------------------------------ 1586 INTEGER :: i 1587 LOGICAL :: ConvertToPTriangle 1588 REAL (KIND=dp) :: uq, vq, sq 1589! INTEGER :: thread, omp_get_thread_num 1590!------------------------------------------------------------------------------- 1591 ConvertToPtriangle = .FALSE. 1592 IF ( PRESENT(PReferenceElement) ) THEN 1593 ConvertToPTriangle = PReferenceElement 1594 END IF 1595 1596 IF ( .NOT. GInit ) CALL GaussPointsInit 1597! thread = 1 1598! !$ thread = omp_get_thread_num()+1 1599! p => IntegStuff(thread) 1600 p => IntegStuff 1601 1602 SELECT CASE (n) 1603 CASE (1) 1604 p % u(1:n) = UTriangle1P 1605 p % v(1:n) = VTriangle1P 1606 p % s(1:n) = STriangle1P / 2.0D0 1607 p % n = 1 1608 CASE (3) 1609 p % u(1:n) = UTriangle3P 1610 p % v(1:n) = VTriangle3P 1611 p % s(1:n) = STriangle3P / 2.0D0 1612 p % n = 3 1613 CASE (4) 1614 p % u(1:n) = UTriangle4P 1615 p % v(1:n) = VTriangle4P 1616 p % s(1:n) = STriangle4P / 2.0D0 1617 p % n = 4 1618 CASE (6) 1619 p % u(1:n) = UTriangle6P 1620 p % v(1:n) = VTriangle6P 1621 p % s(1:n) = STriangle6P / 2.0D0 1622 p % n = 6 1623 CASE (7) 1624 p % u(1:n) = UTriangle7P 1625 p % v(1:n) = VTriangle7P 1626 p % s(1:n) = STriangle7P / 2.0D0 1627 p % n = 7 1628 CASE (11) 1629 p % u(1:n) = UTriangle11P 1630 p % v(1:n) = VTriangle11P 1631 p % s(1:n) = STriangle11P 1632 p % n = 11 1633 CASE (12) 1634 p % u(1:n) = UTriangle12P 1635 p % v(1:n) = VTriangle12P 1636 p % s(1:n) = STriangle12P 1637 p % n = 12 1638 CASE (17) 1639 p % u(1:n) = UTriangle17P 1640 p % v(1:n) = VTriangle17P 1641 p % s(1:n) = STriangle17P 1642 p % n = 17 1643 CASE (20) 1644 p % u(1:n) = UTriangle20P 1645 p % v(1:n) = VTriangle20P 1646 p % s(1:n) = STriangle20P 1647 p % n = 20 1648 CASE DEFAULT 1649! CALL Error( 'GaussPointsTriangle', 'Invalid number of points requested' ) 1650! p % n = 0 1651 1652 p = GaussPointsQuad( n ) 1653!DIR$ IVDEP 1654 DO i=1,p % n 1655 p % v(i) = (p % v(i) + 1) / 2 1656 p % u(i) = (p % u(i) + 1) / 2 * (1 - p % v(i)) 1657 p % s(i) = p % s(i) * ( 1 - p % v(i) ) 1658 END DO 1659 p % s(1:p % n) = 0.5d0 * p % s(1:p % n) / SUM( p % s(1:p % n) ) 1660 END SELECT 1661 1662 IF (ConvertToPTriangle) THEN 1663 !------------------------------------------------------------------- 1664 ! Apply an additional transformation if the actual reference element 1665 ! is the equilateral triangle used in the description of p-elements. 1666 ! We map the original integration points into their counterparts on the 1667 ! p-reference element and scale the weights by the determinant of the 1668 ! deformation gradient associated with the change of reference element. 1669 !------------------------------------------------------------------- 1670!DIR$ IVDEP 1671 DO i=1,P % n 1672 uq = P % u(i) 1673 vq = P % v(i) 1674 sq = P % s(i) 1675 P % u(i) = -1.0d0 + 2.0d0*uq + vq 1676 P % v(i) = SQRT(3.0d0)*vq 1677 P % s(i) = SQRT(3.0d0)*2.0d0*sq 1678 END DO 1679 END IF 1680 1681 p % w(1:n) = 0.0d0 1682!------------------------------------------------------------------------------ 1683 END FUNCTION GaussPointsTriangle 1684!------------------------------------------------------------------------------ 1685 1686 1687!------------------------------------------------------------------------------ 1688!> Return Gaussian integration points for 2D quad element 1689!------------------------------------------------------------------------------ 1690 FUNCTION GaussPointsQuad( np, PMethod) RESULT(p) 1691!------------------------------------------------------------------------------ 1692 INTEGER :: np !< number of points in the requested rule 1693 LOGICAL, OPTIONAL :: PMethod 1694 TYPE(GaussIntegrationPoints_t), POINTER :: p 1695!------------------------------------------------------------------------------ 1696 LOGICAL :: Economic 1697 INTEGER i,j,n,t 1698! INTEGER :: thread, omp_get_thread_num 1699!------------------------------------------------------------------------------ 1700 Economic = .FALSE. 1701 IF (PRESENT(PMethod)) Economic = PMethod 1702 1703 IF ( .NOT. GInit ) CALL GaussPointsInit 1704! thread = 1 1705! !$ thread = omp_get_thread_num()+1 1706! p => IntegStuff(thread) 1707 p => IntegStuff 1708 1709 IF (Economic .AND. (np > 4) .AND. (np <= 60)) THEN 1710 !PRINT *, 'SELECTING A SPECIAL QUADRATURE FOR p-ELEMENTS' 1711 SELECT CASE(np) 1712 CASE (8) 1713 p % u(1:np) = UPQuad8(1:np) 1714 p % v(1:np) = VPQuad8(1:np) 1715 p % s(1:np) = SPQuad8(1:np) 1716 p % n = 8 1717 CASE(12) 1718 p % u(1:np) = UPQuad12(1:np) 1719 p % v(1:np) = VPQuad12(1:np) 1720 p % s(1:np) = SPQuad12(1:np) 1721 p % n = 12 1722 CASE(20) 1723 p % u(1:np) = UPQuad20(1:np) 1724 p % v(1:np) = VPQuad20(1:np) 1725 p % s(1:np) = SPQuad20(1:np) 1726 p % n = 20 1727 CASE(25) 1728 p % u(1:np) = UPQuad25(1:np) 1729 p % v(1:np) = VPQuad25(1:np) 1730 p % s(1:np) = SPQuad25(1:np) 1731 p % n = 25 1732 CASE(36) 1733 p % u(1:np) = UPQuad36(1:np) 1734 p % v(1:np) = VPQuad36(1:np) 1735 p % s(1:np) = SPQuad36(1:np) 1736 p % n = 36 1737 CASE(45) 1738 p % u(1:np) = UPQuad45(1:np) 1739 p % v(1:np) = VPQuad45(1:np) 1740 p % s(1:np) = SPQuad45(1:np) 1741 p % n = 45 1742 CASE(60) 1743 p % u(1:np) = UPQuad60(1:np) 1744 p % v(1:np) = VPQuad60(1:np) 1745 p % s(1:np) = SPQuad60(1:np) 1746 p % n = 60 1747 CASE DEFAULT 1748 WRITE( Message, * ) 'Invalid number of points for p-quadrature: ', np 1749 CALL Fatal('GaussPointsQuad', Message ) 1750 END SELECT 1751 1752 p % w(1:np) = 0.0_dp 1753 !PRINT *, 'NUMBER OF G-POINTS=',p % n 1754 RETURN 1755 END IF 1756 1757 n = SQRT( REAL(np) ) + 0.5 1758 1759 ! WRITE (*,*) 'Integration:', n, np 1760 1761 IF ( n < 1 .OR. n > MAXN ) THEN 1762 p % n = 0 1763 WRITE( Message, * ) 'Invalid number of points: ', n 1764 CALL Error( 'GaussPointsQuad', Message ) 1765 RETURN 1766 END IF 1767 1768 t = 0 1769 DO i=1,n 1770!DIR$ IVDEP 1771 DO j=1,n 1772 t = t + 1 1773 p % u(t) = Points(j,n) 1774 p % v(t) = Points(i,n) 1775 p % s(t) = Weights(i,n)*Weights(j,n) 1776 END DO 1777 END DO 1778 p % n = t 1779!------------------------------------------------------------------------------ 1780 END FUNCTION GaussPointsQuad 1781!------------------------------------------------------------------------------ 1782 1783 1784!------------------------------------------------------------------------------ 1785 FUNCTION GaussPointsPTetra(np) RESULT(p) 1786!------------------------------------------------------------------------------ 1787 INTEGER :: i,np,n 1788 TYPE(GaussIntegrationPoints_t), POINTER :: p 1789 REAL(KIND=dp) :: uh, vh, wh, sh 1790! INTEGER :: thread, omp_get_thread_num 1791!------------------------------------------------------------------------------ 1792 IF ( .NOT. GInit ) CALL GaussPointsInit 1793! thread = 1 1794! !$ thread = omp_get_thread_num()+1 1795! p => IntegStuff(thread) 1796 p => IntegStuff 1797 n = DBLE(np)**(1.0D0/3.0D0) + 0.5D0 1798 1799 ! Get Gauss points of p brick 1800 ! (take into account term z^2) from jacobian determinant 1801 p = GaussPointsPBrick(n,n,n+1) 1802 ! p = GaussPointsBrick( np ) 1803 ! WRITE (*,*) 'Getting Gauss points for: ', n, p % n 1804 1805 ! For each point apply mapping from brick to 1806 ! tetrahedron and multiply each weight by detJ 1807 ! of mapping 1808!DIR$ IVDEP 1809 DO i=1,p % n 1810 uh = p % u(i) 1811 vh = p % v(i) 1812 wh = p % w(i) 1813 sh = p % s(i) 1814 1815 p % u(i)= 1d0/4*(uh - uh*vh - uh*wh + uh*vh*wh) 1816 p % v(i)= SQRT(3d0)/4*(5d0/3 + vh - wh/3 - vh*wh) 1817 p % w(i)= SQRT(6d0)/3*(1d0 + wh) 1818 p % s(i)= -sh * SQRT(2d0)/16 * (1d0 - vh - wh + vh*wh) * (-1d0 + wh) 1819 END DO 1820!------------------------------------------------------------------------------ 1821 END FUNCTION GaussPointsPTetra 1822!------------------------------------------------------------------------------ 1823 1824!------------------------------------------------------------------------------ 1825!> Return Gaussian integration points for 3D tetrahedral element. Here the 1826!> reference element over which the integration is done can also be the 1827!> regular tetrahedron used in the description of p-elements. In that case, 1828!> this routine may return a more economical set of integration points. 1829!------------------------------------------------------------------------------ 1830 FUNCTION GaussPointsTetra( n, PReferenceElement ) RESULT(p) 1831!------------------------------------------------------------------------------ 1832 INTEGER :: n !< number of points in the requested rule 1833 LOGICAL, OPTIONAL :: PReferenceElement !< used for switching the reference element 1834 TYPE(GaussIntegrationPoints_t), POINTER :: p 1835!------------------------------------------------------------------------------ 1836 REAL( KIND=dp ) :: ScaleFactor 1837 INTEGER :: i 1838 LOGICAL :: ConvertToPTetrahedron 1839 REAL (KIND=dp) :: uq, vq, wq, sq 1840! INTEGER :: thread, omp_get_thread_num 1841!---------------------------------------------------------------------------------- 1842 ConvertToPTetrahedron = .FALSE. 1843 IF ( PRESENT(PReferenceElement) ) THEN 1844 ConvertToPTetrahedron = PReferenceElement 1845 END IF 1846 1847 IF ( .NOT. GInit ) CALL GaussPointsInit 1848! thread = 1 1849! !$ thread = omp_get_thread_num()+1 1850! p => IntegStuff(thread) 1851 p => IntegStuff 1852 1853 SELECT CASE (n) 1854 CASE (1) 1855 p % u(1:n) = UTetra1P 1856 p % v(1:n) = VTetra1P 1857 p % w(1:n) = WTetra1P 1858 p % s(1:n) = STetra1P / 6.0D0 1859 p % n = 1 1860 CASE (4) 1861 p % u(1:n) = UTetra4P 1862 p % v(1:n) = VTetra4P 1863 p % w(1:n) = WTetra4P 1864 p % s(1:n) = STetra4P / 6.0D0 1865 p % n = 4 1866 CASE (5) 1867 p % u(1:n) = UTetra5P 1868 p % v(1:n) = VTetra5P 1869 p % w(1:n) = WTetra5P 1870 p % s(1:n) = STetra5P / 6.0D0 1871 p % n = 5 1872 CASE (11) 1873 p % u(1:n) = UTetra11P 1874 p % v(1:n) = VTetra11P 1875 p % w(1:n) = WTetra11P 1876 p % s(1:n) = STetra11P / 6.0D0 1877 p % n = 11 1878 CASE DEFAULT 1879! CALL Error( 'GaussPointsTetra', 'Invalid number of points requested.' ) 1880! p % n = 0 1881 p = GaussPointsBrick( n ) 1882 1883!DIR$ IVDEP 1884 DO i=1,p % n 1885 ScaleFactor = 0.5d0 1886 p % u(i) = ( p % u(i) + 1 ) * Scalefactor 1887 p % v(i) = ( p % v(i) + 1 ) * ScaleFactor 1888 p % w(i) = ( p % w(i) + 1 ) * ScaleFactor 1889 p % s(i) = p % s(i) * ScaleFactor**3 1890 1891 ScaleFactor = 1.0d0 - p % w(i) 1892 p % u(i) = p % u(i) * ScaleFactor 1893 p % v(i) = p % v(i) * ScaleFactor 1894 p % s(i) = p % s(i) * ScaleFactor**2 1895 1896 ScaleFactor = 1.0d0 - p % v(i) / ScaleFactor 1897 p % u(i) = p % u(i) * ScaleFactor 1898 p % s(i) = p % s(i) * ScaleFactor 1899 END DO 1900! p % s(1:p % n) = p % s(1:p % n) / SUM( p % s(1:p % n) ) / 6.0d0 1901 END SELECT 1902 1903 IF (ConvertToPTetrahedron) THEN 1904 !------------------------------------------------------------------- 1905 ! Apply an additional transformation if the actual reference element 1906 ! is the regular tetrahedron used in the description of p-elements 1907 ! We map the original integration points into their counterparts on the 1908 ! p-reference element and scale the weights by the determinant of the 1909 ! deformation gradient associated with the change of reference element. 1910 !------------------------------------------------------------------- 1911!DIR$ IVDEP 1912 DO i=1,P % n 1913 uq = P % u(i) 1914 vq = P % v(i) 1915 wq = P % w(i) 1916 sq = P % s(i) 1917 P % u(i) = -1.0d0 + 2.0d0*uq + vq + wq 1918 P % v(i) = SQRT(3.0d0)*vq + 1.0d0/SQRT(3.0d0)*wq 1919 P % w(i) = SQRT(8.0d0)/SQRT(3.0d0)*wq 1920 P % s(i) = SQRT(8.0d0)*2.0d0*sq 1921 END DO 1922 END IF 1923!------------------------------------------------------------------------------ 1924 END FUNCTION GaussPointsTetra 1925!------------------------------------------------------------------------------ 1926 1927!------------------------------------------------------------------------------ 1928 FUNCTION GaussPointsPPyramid( np ) RESULT(p) 1929!------------------------------------------------------------------------------ 1930 INTEGER :: np,n,i 1931 REAL(KIND=dp) :: uh,vh,wh,sh 1932 TYPE(GaussIntegrationPoints_t), POINTER :: p 1933! INTEGER :: thread, omp_get_thread_num 1934!------------------------------------------------------------------------------ 1935 IF ( .NOT. GInit ) CALL GaussPointsInit 1936! thread = 1 1937! !$ thread = omp_get_thread_num()+1 1938! p => IntegStuff(thread) 1939 p => IntegStuff 1940 1941 n = DBLE(np)**(1.0D0/3.0D0) + 0.5D0 1942 1943 ! Get Gauss points of p brick 1944 ! (take into account term (-1+z)^2) from jacobian determinant 1945 p = GaussPointsPBrick(n,n,n+1) 1946 1947 ! For each point apply mapping from brick to 1948 ! pyramid and multiply each weight by detJ 1949 ! of mapping 1950!DIR$ IVDEP 1951 DO i=1,p % n 1952 uh = p % u(i) 1953 vh = p % v(i) 1954 wh = p % w(i) 1955 sh = p % s(i) 1956 1957 p % u(i)= 1d0/2*uh*(1d0-wh) 1958 p % v(i)= 1d0/2*vh*(1d0-wh) 1959 p % w(i)= SQRT(2d0)/2*(1d0+wh) 1960 p % s(i)= sh * SQRT(2d0)/8 * (-1d0+wh)**2 1961 END DO 1962!------------------------------------------------------------------------------ 1963 END FUNCTION GaussPointsPPyramid 1964!------------------------------------------------------------------------------ 1965 1966 1967!------------------------------------------------------------------------------ 1968!> Return Gaussian integration points for 3D pyramid element. 1969!------------------------------------------------------------------------------ 1970 FUNCTION GaussPointsPyramid( np ) RESULT(p) 1971!------------------------------------------------------------------------------ 1972 INTEGER :: np !< number of points in the requested rule 1973 TYPE(GaussIntegrationPoints_t), POINTER :: p 1974!------------------------------------------------------------------------------ 1975 INTEGER :: i,j,k,n,t 1976! INTEGER :: thread, omp_get_thread_num 1977 1978 IF ( .NOT. GInit ) CALL GaussPointsInit 1979! thread = 1 1980! !$ thread = omp_get_thread_num()+1 1981! p => IntegStuff(thread) 1982 p => IntegStuff 1983 1984 n = REAL(np)**(1.0D0/3.0D0) + 0.5D0 1985 1986 IF ( n < 1 .OR. n > MAXN ) THEN 1987 p % n = 0 1988 WRITE( Message, * ) 'Invalid number of points: ', n 1989 CALL Error( 'GaussPointsPyramid', Message ) 1990 RETURN 1991 END IF 1992 1993 t = 0 1994 DO i=1,n 1995 DO j=1,n 1996!DIR$ IVDEP 1997 DO k=1,n 1998 t = t + 1 1999 p % u(t) = Points(k,n) 2000 p % v(t) = Points(j,n) 2001 p % w(t) = Points(i,n) 2002 p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n) 2003 END DO 2004 END DO 2005 END DO 2006 p % n = t 2007 2008!DIR$ IVDEP 2009 DO t=1,p % n 2010 p % w(t) = (p % w(t) + 1.0d0) / 2.0d0 2011 p % u(t) = p % u(t) * (1.0d0-p % w(t)) 2012 p % v(t) = p % v(t) * (1.0d0-p % w(t)) 2013 p % s(t) = p % s(t) * (1.0d0-p % w(t))**2/2 2014 END DO 2015!------------------------------------------------------------------------------ 2016 END FUNCTION GaussPointsPyramid 2017!------------------------------------------------------------------------------ 2018 2019!------------------------------------------------------------------------------ 2020 FUNCTION GaussPointsPWedge(n) RESULT(p) 2021!------------------------------------------------------------------------------ 2022 INTEGER :: n, i 2023 REAL(KIND=dp) :: uh,vh,wh,sh 2024 TYPE(GaussIntegrationPoints_t), POINTER :: p 2025! INTEGER :: thread, omp_get_thread_num 2026!------------------------------------------------------------------------------ 2027 IF ( .NOT. GInit ) CALL GaussPointsInit 2028! thread = 1 2029! !$ thread = omp_get_thread_num()+1 2030! p => IntegStuff(thread) 2031 p => IntegStuff 2032 2033 ! Get Gauss points of brick 2034 p = GaussPointsBrick(n) 2035 2036 ! For each point apply mapping from brick to 2037 ! wedge and multiply each weight by detJ 2038 ! of mapping 2039!DIR$ IVDEP 2040 DO i=1,p % n 2041 uh = p % u(i) 2042 vh = p % v(i) 2043 wh = p % w(i) 2044 sh = p % s(i) 2045 2046 p % u(i)= 1d0/2*(uh-uh*vh) 2047 p % v(i)= SQRT(3d0)/2*(1d0+vh) 2048 p % w(i)= wh 2049 p % s(i)= sh * SQRT(3d0)*(1-vh)/4 2050 END DO 2051!------------------------------------------------------------------------------ 2052 END FUNCTION GaussPointsPWedge 2053!------------------------------------------------------------------------------ 2054 2055 2056 2057!------------------------------------------------------------------------------ 2058!> Return Gaussian integration points for 3D wedge element 2059!------------------------------------------------------------------------------ 2060 FUNCTION GaussPointsWedge( np ) RESULT(p) 2061!------------------------------------------------------------------------------ 2062 INTEGER :: np !< number of points in the requested rule 2063 TYPE(GaussIntegrationPoints_t), POINTER :: p 2064!------------------------------------------------------------------------------ 2065 INTEGER :: i,j,k,n,t 2066! INTEGER :: thread, omp_get_thread_num 2067!------------------------------------------------------------------------------ 2068 IF ( .NOT. GInit ) CALL GaussPointsInit 2069! thread = 1 2070! !$ thread = omp_get_thread_num()+1 2071! p => IntegStuff(thread) 2072 p => IntegStuff 2073 2074 n = REAL(np)**(1.0d0/3.0d0) + 0.5d0 2075 2076 IF ( n < 1 .OR. n > MAXN ) THEN 2077 p % n = 0 2078 WRITE( Message, * ) 'Invalid number of points: ', n 2079 CALL Error( 'GaussPointsWedge', Message ) 2080 RETURN 2081 END IF 2082 2083 t = 0 2084 DO i=1,n 2085 DO j=1,n 2086!DIR$ IVDEP 2087 DO k=1,n 2088 t = t + 1 2089 p % u(t) = Points(k,n) 2090 p % v(t) = Points(j,n) 2091 p % w(t) = Points(i,n) 2092 p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n) 2093 END DO 2094 END DO 2095 END DO 2096 p % n = t 2097 2098!DIR$ IVDEP 2099 DO i=1,p % n 2100 p % v(i) = (p % v(i) + 1)/2 2101 p % u(i) = (p % u(i) + 1)/2 * (1 - p % v(i)) 2102 p % s(i) = p % s(i) * (1-p % v(i))/4 2103 END DO 2104!------------------------------------------------------------------------------ 2105 END FUNCTION GaussPointsWedge 2106!------------------------------------------------------------------------------ 2107 2108!------------------------------------------------------------------------------ 2109!> Return an optimized number of Gaussian points for integrating over prisms. 2110!> Here the reference element can also be that of the p-approximation. 2111!> A rule with m x n points is returned. 2112!------------------------------------------------------------------------------ 2113 FUNCTION GaussPointsWedge2(m,n,PReferenceElement) RESULT(p) 2114!------------------------------------------------------------------------------ 2115 TYPE(GaussIntegrationPoints_t), POINTER :: p 2116 INTEGER :: m !< The number of points over a triangular face 2117 INTEGER :: n !< The number of points in the orthogonal direction to the triangular faces 2118 LOGICAL, OPTIONAL :: PReferenceElement !< Used for switching the reference element 2119!------------------------------------------------------------------------------ 2120 INTEGER :: i,j,k,t 2121 LOGICAL :: ConvertToPPrism 2122 REAL(KIND=dp) :: uq(20), vq(20), sq(20) 2123 REAL(KIND=dp) :: uh,vh,wh,sh 2124!----------------------------------------------------------------------------- 2125 ConvertToPPrism = .FALSE. 2126 IF ( PRESENT(PReferenceElement) ) THEN 2127 ConvertToPPrism = PReferenceElement 2128 END IF 2129 IF ( .NOT. GInit ) CALL GaussPointsInit 2130 p => IntegStuff 2131 2132 SELECT CASE (m) 2133 CASE (1) 2134 IF (ConvertToPPrism) THEN 2135 uq(1) = -1.0d0 + 2.0d0*UTriangle1P(1) + VTriangle1P(1) 2136 vq(1) = SQRT(3.0d0) * VTriangle1P(1) 2137 sq(1) = SQRT(3.0d0) * STriangle1P(1) 2138 ELSE 2139 uq(1:m) = UTriangle1P 2140 vq(1:m) = VTriangle1P 2141 sq(1:m) = STriangle1P / 2.0D0 2142 END IF 2143 CASE (3) 2144 IF (ConvertToPPrism) THEN 2145!DIR$ IVDEP 2146 DO i=1,m 2147 uq(i) = -1.0d0 + 2.0d0*UTriangle3P(i) + VTriangle3P(i) 2148 vq(i) = SQRT(3.0d0) * VTriangle3P(i) 2149 sq(i) = SQRT(3.0d0) * STriangle3P(i) 2150 END DO 2151 ELSE 2152 uq(1:m) = UTriangle3P 2153 vq(1:m) = VTriangle3P 2154 sq(1:m) = STriangle3P / 2.0D0 2155 END IF 2156 CASE (4) 2157 IF (ConvertToPPrism) THEN 2158!DIR$ IVDEP 2159 DO i=1,m 2160 uq(i) = -1.0d0 + 2.0d0*UTriangle4P(i) + VTriangle4P(i) 2161 vq(i) = SQRT(3.0d0) * VTriangle4P(i) 2162 sq(i) = SQRT(3.0d0) * STriangle4P(i) 2163 END DO 2164 ELSE 2165 uq(1:m) = UTriangle4P 2166 vq(1:m) = VTriangle4P 2167 sq(1:m) = STriangle4P / 2.0D0 2168 END IF 2169 CASE (6) 2170 IF (ConvertToPPrism) THEN 2171!DIR$ IVDEP 2172 DO i=1,m 2173 uq(i) = -1.0d0 + 2.0d0*UTriangle6P(i) + VTriangle6P(i) 2174 vq(i) = SQRT(3.0d0) * VTriangle6P(i) 2175 sq(i) = SQRT(3.0d0) * STriangle6P(i) 2176 END DO 2177 ELSE 2178 uq(1:m) = UTriangle6P 2179 vq(1:m) = VTriangle6P 2180 sq(1:m) = STriangle6P / 2.0D0 2181 END IF 2182 CASE (7) 2183 IF (ConvertToPPrism) THEN 2184!DIR$ IVDEP 2185 DO i=1,m 2186 uq(i) = -1.0d0 + 2.0d0*UTriangle7P(i) + VTriangle7P(i) 2187 vq(i) = SQRT(3.0d0) * VTriangle7P(i) 2188 sq(i) = SQRT(3.0d0) * STriangle7P(i) 2189 END DO 2190 ELSE 2191 uq(1:m) = UTriangle7P 2192 vq(1:m) = VTriangle7P 2193 sq(1:m) = STriangle7P / 2.0D0 2194 END IF 2195 CASE (11) 2196 IF (ConvertToPPrism) THEN 2197!DIR$ IVDEP 2198 DO i=1,m 2199 uq(i) = -1.0d0 + 2.0d0*UTriangle11P(i) + VTriangle11P(i) 2200 vq(i) = SQRT(3.0d0) * VTriangle11P(i) 2201 sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle11P(i) 2202 END DO 2203 ELSE 2204 uq(1:m) = UTriangle11P 2205 vq(1:m) = VTriangle11P 2206 sq(1:m) = STriangle11P 2207 END IF 2208 CASE (12) 2209 IF (ConvertToPPrism) THEN 2210!DIR$ IVDEP 2211 DO i=1,m 2212 uq(i) = -1.0d0 + 2.0d0*UTriangle12P(i) + VTriangle12P(i) 2213 vq(i) = SQRT(3.0d0) * VTriangle12P(i) 2214 sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle12P(i) 2215 END DO 2216 ELSE 2217 uq(1:m) = UTriangle12P 2218 vq(1:m) = VTriangle12P 2219 sq(1:m) = STriangle12P 2220 END IF 2221 CASE (17) 2222 IF (ConvertToPPrism) THEN 2223!DIR$ IVDEP 2224 DO i=1,m 2225 uq(i) = -1.0d0 + 2.0d0*UTriangle17P(i) + VTriangle17P(i) 2226 vq(i) = SQRT(3.0d0) * VTriangle17P(i) 2227 sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle17P(i) 2228 END DO 2229 ELSE 2230 uq(1:m) = UTriangle17P 2231 vq(1:m) = VTriangle17P 2232 sq(1:m) = STriangle17P 2233 END IF 2234 CASE (20) 2235 IF (ConvertToPPrism) THEN 2236!DIR$ IVDEP 2237 DO i=1,m 2238 uq(i) = -1.0d0 + 2.0d0*UTriangle20P(i) + VTriangle20P(i) 2239 vq(i) = SQRT(3.0d0) * VTriangle20P(i) 2240 sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle20P(i) 2241 END DO 2242 ELSE 2243 uq(1:m) = UTriangle20P 2244 vq(1:m) = VTriangle20P 2245 sq(1:m) = STriangle20P 2246 END IF 2247 CASE DEFAULT 2248 !------------------------------------------------------------------------------------- 2249 ! The requested number of integration points associated with the triangular face 2250 ! cannot be created. Therefore use the given number of points in the third direction 2251 ! and generate silently n x n x n rule. 2252 !----------------------------------------------------------------------------------- 2253 IF (ConvertToPPrism) THEN 2254 p = GaussPointsBrick(n*n*n) 2255!DIR$ IVDEP 2256 DO i=1,p % n 2257 uh = p % u(i) 2258 vh = p % v(i) 2259 wh = p % w(i) 2260 sh = p % s(i) 2261 2262 p % u(i)= 1d0/2*(uh-uh*vh) 2263 p % v(i)= SQRT(3d0)/2*(1d0+vh) 2264 p % w(i)= wh 2265 p % s(i)= sh * SQRT(3d0)*(1-vh)/4 2266 END DO 2267 ELSE 2268 t = 0 2269 DO i=1,n 2270 DO j=1,n 2271!DIR$ IVDEP 2272 DO k=1,n 2273 t = t + 1 2274 p % u(t) = Points(k,n) 2275 p % v(t) = Points(j,n) 2276 p % w(t) = Points(i,n) 2277 p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n) 2278 END DO 2279 END DO 2280 END DO 2281 p % n = t 2282 2283!DIR$ IVDEP 2284 DO i=1,p % n 2285 p % v(i) = (p % v(i) + 1)/2 2286 p % u(i) = (p % u(i) + 1)/2 * (1 - p % v(i)) 2287 p % s(i) = p % s(i) * (1-p % v(i))/4 2288 END DO 2289 END IF 2290 RETURN 2291 END SELECT 2292 2293 t = 0 2294 DO i=1,m 2295!DIR$ IVDEP 2296 DO j=1,n 2297 t = t + 1 2298 p % u(t) = uq(i) 2299 p % v(t) = vq(i) 2300 p % w(t) = Points(j,n) 2301 p % s(t) = sq(i)*Weights(j,n) 2302 END DO 2303 END DO 2304 p % n = t 2305!------------------------------------------------------------------------------ 2306 END FUNCTION GaussPointsWedge2 2307!------------------------------------------------------------------------------ 2308 2309 2310!------------------------------------------------------------------------------ 2311!> Return Gaussian integration points for 3D wedge elements using 2312!> economical quadratures that are not product of segment and triangle rules. 2313!------------------------------------------------------------------------------ 2314 FUNCTION GaussPointsWedgeEconomic( n, PReferenceElement ) RESULT(p) 2315!------------------------------------------------------------------------------ 2316 INTEGER :: n !< number of points in the requested rule 2317 LOGICAL, OPTIONAL :: PReferenceElement !< used for switching the reference element 2318 TYPE(GaussIntegrationPoints_t), POINTER :: p 2319!------------------------------------------------------------------------------ 2320 REAL( KIND=dp ) :: ScaleFactor 2321 INTEGER :: i 2322 LOGICAL :: ConvertToPWedge 2323 REAL (KIND=dp) :: uq, vq, wq, sq 2324! INTEGER :: thread, omp_get_thread_num 2325!---------------------------------------------------------------------------------- 2326 ConvertToPWedge = .FALSE. 2327 IF ( PRESENT(PReferenceElement) ) THEN 2328 ConvertToPWedge = PReferenceElement 2329 END IF 2330 2331 IF ( .NOT. GInit ) CALL GaussPointsInit 2332! thread = 1 2333! !$ thread = omp_get_thread_num()+1 2334! p => IntegStuff(thread) 2335 p => IntegStuff 2336 2337 SELECT CASE (n) 2338 CASE (4) 2339 p % u(1:n) = UWedge4P 2340 p % v(1:n) = VWedge4P 2341 p % w(1:n) = WWedge4P 2342 p % s(1:n) = SWedge4P 2343 CASE (5) 2344 p % u(1:n) = UWedge5P 2345 p % v(1:n) = VWedge5P 2346 p % w(1:n) = WWedge5P 2347 p % s(1:n) = SWedge5P 2348 CASE (6) 2349 p % u(1:n) = UWedge6P 2350 p % v(1:n) = VWedge6P 2351 p % w(1:n) = WWedge6P 2352 p % s(1:n) = SWedge6P 2353 CASE (7) 2354 p % u(1:n) = UWedge7P 2355 p % v(1:n) = VWedge7P 2356 p % w(1:n) = WWedge7P 2357 p % s(1:n) = SWedge7P 2358 CASE (10) 2359 p % u(1:n) = UWedge10P 2360 p % v(1:n) = VWedge10P 2361 p % w(1:n) = WWedge10P 2362 p % s(1:n) = SWedge10P 2363 CASE (11) 2364 p % u(1:n) = UWedge11P 2365 p % v(1:n) = VWedge11P 2366 p % w(1:n) = WWedge11P 2367 p % s(1:n) = SWedge11P 2368 CASE (14) 2369 p % u(1:n) = UWedge14P 2370 p % v(1:n) = VWedge14P 2371 p % w(1:n) = WWedge14P 2372 p % s(1:n) = SWedge14P 2373 CASE (15) 2374 p % u(1:n) = UWedge15P 2375 p % v(1:n) = VWedge15P 2376 p % w(1:n) = WWedge15P 2377 p % s(1:n) = SWedge15P 2378 CASE (16) 2379 p % u(1:n) = UWedge16P 2380 p % v(1:n) = VWedge16P 2381 p % w(1:n) = WWedge16P 2382 p % s(1:n) = SWedge16P 2383 CASE (24) 2384 p % u(1:n) = UWedge24P 2385 p % v(1:n) = VWedge24P 2386 p % w(1:n) = WWedge24P 2387 p % s(1:n) = SWedge24P 2388 2389 CASE DEFAULT 2390 CALL Fatal( 'GaussPointsWedgeEconomic',& 2391 'Invalid number of points requested') 2392 END SELECT 2393 2394 p % n = n 2395 2396 IF (ConvertToPWedge ) THEN 2397 DO i=1,P % n 2398 uq = P % u(i) 2399 vq = P % v(i) 2400 sq = P % s(i) 2401 2402 ! first to classical 2403 uq = (uq+1.d0)/2.0d0 2404 vq = (vq+1.d0)/2.0d0 2405 2406 ! then to p-convention 2407 P % u(i) = -1.0d0 + 2.0d0*uq + vq 2408 P % v(i) = SQRT(3.0d0) * vq 2409 P % s(i) = SQRT(3.0d0) * sq / 2.0d0 2410 END DO 2411 ELSE 2412 ! Map to classical Elmer local coordinates in [0,1] 2413 p % u(1:n) = ( p % u(1:n)+1.0d0 ) / 2.0d0 2414 p % v(1:n) = ( p % v(1:n)+1.0d0 ) / 2.0d0 2415 p % s(1:n) = p % s(1:n) / 4.0d0 2416 END IF 2417!------------------------------------------------------------------------------ 2418 END FUNCTION GaussPointsWedgeEconomic 2419!------------------------------------------------------------------------------ 2420 2421 2422 2423!------------------------------------------------------------------------------ 2424!> Return Gaussian integration points for 3D brick element for 2425!> composite rule 2426!> sum_i=1^nx(sum_j=1^ny(sum_k=1^nz w_ijk f(x_{i,j,k},y_{i,j,k},z_{i,j,k}))). 2427!------------------------------------------------------------------------------ 2428 FUNCTION GaussPointsPBrick( nx, ny, nz ) RESULT(p) 2429!------------------------------------------------------------------------------ 2430 INTEGER :: nx !< number of points in the requested rule in x direction 2431 INTEGER :: ny !< number of points in the requested rule in y direction 2432 INTEGER :: nz !< number of points in the requested rule in z direction 2433 TYPE(GaussIntegrationPoints_t), POINTER :: p 2434!------------------------------------------------------------------------------ 2435 INTEGER i,j,k,t 2436! INTEGER :: thread, omp_get_thread_num 2437!------------------------------------------------------------------------------ 2438 IF ( .NOT. GInit ) CALL GaussPointsInit 2439! thread = 1 2440! !$ thread = omp_get_thread_num()+1 2441! p => IntegStuff(thread) 2442 p => IntegStuff 2443 2444 ! Check validity of number of integration points 2445 IF ( nx < 1 .OR. nx > MAXN .OR. & 2446 ny < 1 .OR. ny > MAXN .OR. & 2447 nz < 1 .OR. nz > MAXN) THEN 2448 p % n = 0 2449 WRITE( Message, * ) 'Invalid number of points: ', nx, ny, nz 2450 CALL Error( 'GaussPointsBrick', Message ) 2451 RETURN 2452 END IF 2453 2454 t = 0 2455 DO i=1,nx 2456 DO j=1,ny 2457!DIR$ IVDEP 2458 DO k=1,nz 2459 t = t + 1 2460 p % u(t) = Points(i,nx) 2461 p % v(t) = Points(j,ny) 2462 p % w(t) = Points(k,nz) 2463 p % s(t) = Weights(i,nx)*Weights(j,ny)*Weights(k,nz) 2464 END DO 2465 END DO 2466 END DO 2467 p % n = t 2468!------------------------------------------------------------------------------ 2469 END FUNCTION GaussPointsPBrick 2470!------------------------------------------------------------------------------ 2471 2472 2473!------------------------------------------------------------------------------ 2474!> Return Gaussian integration points for 3D brick element 2475!------------------------------------------------------------------------------ 2476 FUNCTION GaussPointsBrick( np ) RESULT(p) 2477!------------------------------------------------------------------------------ 2478 INTEGER :: np !< number of points in the requested rule 2479 TYPE(GaussIntegrationPoints_t), POINTER :: p 2480!------------------------------------------------------------------------------ 2481 INTEGER i,j,k,n,t 2482! INTEGER :: thread, omp_get_thread_num 2483 2484 IF ( .NOT. GInit ) CALL GaussPointsInit 2485! thread = 1 2486! !$ thread = omp_get_thread_num()+1 2487! p => IntegStuff(thread) 2488 p => IntegStuff 2489 2490 n = REAL(np)**(1.0D0/3.0D0) + 0.5D0 2491 2492 IF ( n < 1 .OR. n > MAXN ) THEN 2493 p % n = 0 2494 WRITE( Message, * ) 'Invalid number of points: ', n 2495 CALL Error( 'GaussPointsBrick', Message ) 2496 RETURN 2497 END IF 2498 2499 t = 0 2500 DO i=1,n 2501 DO j=1,n 2502!DIR$ IVDEP 2503 DO k=1,n 2504 t = t + 1 2505 p % u(t) = Points(k,n) 2506 p % v(t) = Points(j,n) 2507 p % w(t) = Points(i,n) 2508 p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n) 2509 END DO 2510 END DO 2511 END DO 2512 p % n = t 2513!------------------------------------------------------------------------------ 2514 END FUNCTION GaussPointsBrick 2515!------------------------------------------------------------------------------ 2516 2517 2518!------------------------------------------------------------------------------ 2519!> Given element structure return Gauss integration points for the element. 2520!---------------------------------------------------------------------------------------------- 2521 FUNCTION GaussPoints( elm, np, RelOrder, EdgeBasis, PReferenceElement, & 2522 EdgeBasisDegree) RESULT(IntegStuff) 2523!--------------------------------------------------------------------------------------------- 2524 USE PElementMaps, ONLY : isActivePElement 2525 TYPE( Element_t ) :: elm !< Structure holding element description 2526 INTEGER, OPTIONAL :: np !< Number of requested points 2527 INTEGER, OPTIONAL :: RelOrder !< Relative order of Gauss points (-1,0,+1) 2528 LOGICAL, OPTIONAL :: EdgeBasis !< For choosing quadrature suitable for edge elements 2529 LOGICAL, OPTIONAL :: PReferenceElement !< For switching to the p-version reference element 2530 INTEGER, OPTIONAL :: EdgeBasisDegree ! The degree of edge elements 2531 TYPE( GaussIntegrationPoints_t ) :: IntegStuff !< Structure holding the integration points 2532!------------------------------------------------------------------------------ 2533 LOGICAL :: pElement, UsePRefElement, Economic 2534 INTEGER :: n, eldim, p1d, ntri, nseg, necon 2535 TYPE(ElementType_t), POINTER :: elmt 2536!------------------------------------------------------------------------------ 2537 elmt => elm % TYPE 2538 2539 IF (PRESENT(EdgeBasis)) THEN 2540 IF (EdgeBasis) THEN 2541 UsePRefElement = .TRUE. 2542 IF (PRESENT(PReferenceElement)) UsePRefElement = PReferenceElement 2543 IF (PRESENT(EdgeBasisDegree)) THEN 2544 IntegStuff = EdgeElementGaussPoints(elmt % ElementCode/100, UsePRefElement, EdgeBasisDegree) 2545 ELSE 2546 IntegStuff = EdgeElementGaussPoints(elmt % ElementCode/100, UsePRefElement) 2547 END IF 2548 RETURN 2549 END IF 2550 END IF 2551 2552 IF( PRESENT(PReferenceElement)) THEN 2553 pElement = PReferenceElement 2554 ELSE 2555 pElement = isActivePElement(elm) 2556 END IF 2557 2558 IF ( PRESENT(np) ) THEN 2559 n = np 2560 ELSE IF( PRESENT( RelOrder ) ) THEN 2561 IF (pElement) THEN 2562 n = elm % PDefs % GaussPoints 2563 IF( RelOrder == 0 ) THEN 2564 CONTINUE 2565 ELSE 2566 eldim = elm % type % dimension 2567 p1d = NINT( n**(1.0_dp/eldim) ) + RelOrder 2568 IF( p1d < 1 ) THEN 2569 CALL Fatal('GaussPoints','Number of integration points must remain positive!') 2570 END IF 2571 n = p1d**eldim 2572 2573 ! Take into account the case of economic integration. Map a tensor 2574 ! product rule to a corresponding economic rule. 2575 IF (elmt % ElementCode / 100 == 4) THEN 2576 SELECT CASE(n) 2577 CASE(9) 2578 n = 8 2579 CASE(16) 2580 n = 12 2581 CASE(25) 2582 n = 20 2583 CASE(36) 2584 n = 25 2585 CASE(49) 2586 n = 36 2587 CASE(64) 2588 n = 45 2589 CASE(81) 2590 n = 60 2591 CASE DEFAULT 2592 CONTINUE 2593 END SELECT 2594 END IF 2595 END IF 2596 ELSE 2597 IF( RelOrder == 0 ) THEN 2598 n = elmt % GaussPoints 2599 ELSE IF( RelOrder == 1 ) THEN 2600 n = elmt % GaussPoints2 2601 ELSE IF( RelOrder == -1 ) THEN 2602 n = elmt % GaussPoints0 2603 ELSE 2604 PRINT *,'RelOrder can only be {-1, 0, 1} !' 2605 END IF 2606 END IF 2607 ELSE 2608 IF (pElement) THEN 2609 IF( ASSOCIATED( elm % PDefs ) ) THEN 2610 n = elm % PDefs % GaussPoints 2611 ELSE 2612 n = elmt % GaussPoints 2613 END IF 2614 IF( n == 0 ) n = elmt % GaussPoints 2615 ELSE 2616 n = elmt % GaussPoints 2617 END IF 2618 END IF 2619 2620 SELECT CASE( elmt % ElementCode / 100 ) 2621 CASE (1) 2622 IntegStuff = GaussPoints0D(n) 2623 2624 CASE (2) 2625 IntegStuff = GaussPoints1D(n) 2626 2627 CASE (3) 2628 IF (pElement) THEN 2629 IntegStuff = GaussPointsPTriangle(n) 2630 ELSE 2631 IntegStuff = GaussPointsTriangle(n) 2632 END IF 2633 2634 CASE (4) 2635 IF (pElement .AND. ASSOCIATED( elm % pdefs ) ) THEN 2636 Economic = .FALSE. 2637 ! For certain polynomial orders, economic quadratures may be used: 2638 IF (elm % PDefs % P > 1 .AND. elm % PDefs % P <= 8) Economic = .TRUE. 2639 ! An explicit bubble augmentation with lower-order methods switches to 2640 ! the standard rule: 2641 IF (elm % BDOFs > 0 .AND. elm % PDefs % P < 4) Economic = .FALSE. 2642 IntegStuff = GaussPointsQuad(n, Economic) 2643 ELSE 2644 IntegStuff = GaussPointsQuad(n) 2645 END IF 2646 2647 CASE (5) 2648 IF (pElement) THEN 2649 IntegStuff = GaussPointsPTetra(n) 2650 ELSE 2651 IntegStuff = GaussPointsTetra(n) 2652 END IF 2653 2654 CASE (6) 2655 IF (pElement) THEN 2656 IntegStuff = GaussPointsPPyramid(n) 2657 ELSE 2658 IntegStuff = GaussPointsPyramid(n) 2659 END IF 2660 2661 CASE (7) 2662 IF( PRESENT( np ) ) THEN 2663 ! possible values: 2664 ! triangle = 1, 3, 4, 6, 7, 11, 12, 17, 20 2665 ! segment = 1, 2, 3, 4, 5, 6, 7, 8, 9, 2666 2667 ntri = 0; nseg = 0; necon = 0 2668 2669 SELECT CASE( n ) 2670 2671 ! Cases ordered so that we have the segment x triangle rules first 2672 CASE( 1, 2, 3) 2673 nseg = n 2674 CASE( 6, 8 ) 2675 nseg = 2 2676 CASE( 12, 18, 21 ) 2677 nseg = 3 2678 CASE( 28, 44, 48 ) 2679 nseg = 4 2680 CASE( 85, 100 ) 2681 nseg = 5 2682 2683 ! The economical rules 2684 CASE( 4, 5, 7, 10, 11, 14, 15, 16, 24 ) 2685 ! Note: we would have 6 and 8 point rules from the economic family as well 2686 necon = n 2687 END SELECT 2688 2689 IF( nseg > 0 ) THEN 2690 ntri = n / nseg 2691 IntegStuff = GaussPointsWedge2(ntri,nseg,PReferenceElement=pElement) 2692 RETURN 2693 ELSE IF( necon > 0 ) THEN 2694 IntegStuff = GaussPointsWedgeEconomic(necon,PReferenceElement=pElement) 2695 RETURN 2696 END IF 2697 END IF 2698 2699 IF (pElement) THEN 2700 IntegStuff = GaussPointsPWedge(n) 2701 ELSE 2702 IntegStuff = GaussPointsWedge(n) 2703 END IF 2704 2705 CASE (8) 2706 IntegStuff = GaussPointsBrick(n) 2707 END SELECT 2708 2709 END FUNCTION GaussPoints 2710 2711 2712!---------------------------------------------------------------------------------- 2713!> Return a suitable version of the Gaussian numerical integration method for 2714!> H(curl)-conforming finite elements. The default here is that the edge basis 2715!> functions are obtained via calling the function EdgeElementInfo, so that 2716!> the reference element is chosen to be that used for p-approximation. 2717!> By giving the optional argument PiolaVersion = .FALSE. this function returns 2718!> an alternate (usually smaller) set of integration points on the classic 2719!> reference element. This option may be useful when the edge basis functions are 2720!> obtained via calling the alternate subroutine GetEdgeBasis. 2721!---------------------------------------------------------------------------------------- 2722 FUNCTION EdgeElementGaussPoints(ElementFamily, PiolaVersion, BasisDegree) RESULT(IP) 2723!--------------------------------------------------------------------------------------- 2724 INTEGER :: ElementFamily 2725 LOGICAL, OPTIONAL :: PiolaVersion 2726 INTEGER, OPTIONAL :: BasisDegree 2727 TYPE(GaussIntegrationPoints_t) :: IP 2728!------------------------------------------------------------------------------ 2729 LOGICAL :: PRefElement, SecondOrder 2730!------------------------------------------------------------------------------ 2731 PRefElement = .TRUE. 2732 SecondOrder = .FALSE. 2733 IF ( PRESENT(PiolaVersion) ) PRefElement = PiolaVersion 2734 IF ( PRESENT(BasisDegree) ) SecondOrder = BasisDegree > 1 2735 2736 SELECT CASE(ElementFamily) 2737 CASE (1) 2738 IP = GaussPoints0D(1) 2739 CASE (2) 2740 IP = GaussPoints1D(2) 2741 CASE(3) 2742 IF (SecondOrder) THEN 2743 IP = GaussPointsTriangle(6, PReferenceElement=PRefElement) 2744 ELSE 2745 IP = GaussPointsTriangle(3, PReferenceElement=PRefElement) 2746 END IF 2747 CASE(4) 2748 IF (PRefElement) THEN 2749 IP = GaussPointsQuad(9) 2750 ELSE 2751 IP = GaussPointsQuad(4) 2752 END IF 2753 CASE(5) 2754 IF (SecondOrder) THEN 2755 IP = GaussPointsTetra(11, PReferenceElement=PRefElement) 2756 ELSE 2757 IP = GaussPointsTetra(4, PReferenceElement=PRefElement) 2758 END IF 2759 CASE(6) 2760 IF (PRefElement) THEN 2761 IP = GaussPointsPPyramid(27) 2762 ELSE 2763 IP = GaussPointsPyramid(27) 2764 END IF 2765 CASE(7) 2766 IF (PRefElement) THEN 2767 IP = GaussPointsWedge2(6,3,PReferenceElement=PRefElement) 2768 ELSE 2769 IP = GaussPointsWedge2(3,2,PReferenceElement=PRefElement) 2770 END IF 2771 CASE(8) 2772 IF (PRefElement) THEN 2773 IP = GaussPointsPBrick(3,3,3) 2774 ELSE 2775 IP = GaussPointsBrick(8) 2776 END IF 2777 CASE DEFAULT 2778 CALL Fatal('Integration::EdgeElementGaussPoints','Unsupported element type') 2779 END SELECT 2780!------------------------------------------------------------------------------ 2781 END FUNCTION EdgeElementGaussPoints 2782!------------------------------------------------------------------------------ 2783 2784 2785!--------------------------------------------------------------------------- 2786END MODULE Integration 2787!--------------------------------------------------------------------------- 2788 2789!> \} ElmerLib 2790