CoCalc -- Integration.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fem/src/Integration.F90
³²⁰³ views
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:   [email protected]
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
!-----------------------------------------------------------------------------
44
MODULE Integration
45
   USE LoadMod
46

47
   IMPLICIT NONE
48

49
   INTEGER, PARAMETER, PRIVATE :: MAXN = 13, MAXNPAD = 16 ! Padded to 64-byte alignment
50
   INTEGER, PARAMETER, PRIVATE :: MAX_INTEGRATION_POINTS = MAXN**3
51

52
   LOGICAL, PRIVATE :: GInit = .FALSE.
53
   !$OMP THREADPRIVATE(GInit)
54

55
!------------------------------------------------------------------------------
56
   TYPE GaussIntegrationPoints_t
57
      INTEGER :: N
58
      REAL(KIND=dp), POINTER CONTIG :: u(:),v(:),w(:),s(:)
59
!DIR$ ATTRIBUTES ALIGN:64 :: u, v, w, s
60
   END TYPE GaussIntegrationPoints_t
61

62
   TYPE(GaussIntegrationPoints_t), TARGET, PRIVATE, SAVE :: IntegStuff
63
   !$OMP THREADPRIVATE(IntegStuff)
64
!------------------------------------------------------------------------------
65

66
!------------------------------------------------------------------------------
67
! Storage for 1d Gauss points, and weights. The values are computed on the
68
! fly (see ComputeGaussPoints1D below). These values are used for quads and
69
! bricks as well. To avoid NUMA issues, Points and Weights are private for each
70
! thread.
71
!------------------------------------------------------------------------------
72
   REAL(KIND=dp), PRIVATE, SAVE :: Points(MAXNPAD,MAXN),Weights(MAXNPAD,MAXN)
73
   !$OMP THREADPRIVATE(Points, Weights)
74
!DIR$ ATTRIBUTES ALIGN:64::Points, Weights
75
!------------------------------------------------------------------------------
76

77
!------------------------------------------------------------------------------
78
! Triangle - 1 point rule; exact integration of x^py^q, p+q<=1
79
!------------------------------------------------------------------------------
80
   REAL(KIND=dp),DIMENSION(1),PRIVATE :: UTriangle1P=(/ 0.3333333333333333D0 /)
81
   REAL(KIND=dp),DIMENSION(1),PRIVATE :: VTriangle1P=(/ 0.3333333333333333D0 /)
82
   REAL(KIND=dp),DIMENSION(1),PRIVATE :: STriangle1P=(/ 1.0000000000000000D0 /)
83

84
!------------------------------------------------------------------------------
85
! Triangle - 3 point rule; exact integration of x^py^q, p+q<=2
86
!------------------------------------------------------------------------------
87
   REAL(KIND=dp), DIMENSION(3), PRIVATE :: UTriangle3P = &
88
    (/ 0.16666666666666667D0, 0.66666666666666667D0, 0.16666666666666667D0 /)
89

90
   REAL(KIND=dp), DIMENSION(3), PRIVATE :: VTriangle3P = &
91
    (/ 0.16666666666666667D0, 0.16666666666666667D0, 0.66666666666666667D0 /)
92

93
   REAL(KIND=dp), DIMENSION(3), PRIVATE :: STriangle3P = &
94
    (/ 0.33333333333333333D0, 0.33333333333333333D0, 0.33333333333333333D0 /)
95

96
!------------------------------------------------------------------------------
97
! Triangle - 4 point rule; exact integration of x^py^q, p+q<=3
98
!------------------------------------------------------------------------------
99
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: UTriangle4P = &
100
         (/  0.33333333333333333D0, 0.2000000000000000D0, &
101
             0.60000000000000000D0, 0.20000000000000000D0 /)
102

103
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: VTriangle4P = &
104
         (/  0.33333333333333333D0, 0.2000000000000000D0, &
105
             0.20000000000000000D0, 0.60000000000000000D0 /)
106

107
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: STriangle4P = &
108
         (/ -0.56250000000000000D0, 0.52083333333333333D0, &
109
             0.52083333333333333D0, 0.52083333333333333D0 /)
110

111
!------------------------------------------------------------------------------
112
! Triangle - 6 point rule; exact integration of x^py^q, p+q<=4
113
!------------------------------------------------------------------------------
114
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: UTriangle6P = &
115
       (/ 0.091576213509771D0, 0.816847572980459D0, 0.091576213509771D0, &
116
          0.445948490915965D0, 0.108103018168070D0, 0.445948490915965D0 /)
117

118
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: VTriangle6P = &
119
        (/ 0.091576213509771D0, 0.091576213509771D0, 0.816847572980459D0, &
120
           0.445948490915965D0, 0.445948490915965D0, 0.108103018168070D0 /)
121

122
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: STriangle6P = &
123
        (/ 0.109951743655322D0, 0.109951743655322D0, 0.109951743655322D0, &
124
           0.223381589678011D0, 0.223381589678011D0, 0.223381589678011D0 /)
125

126
!------------------------------------------------------------------------------
127
! Triangle - 7 point rule; exact integration of x^py^q, p+q<=5
128
!------------------------------------------------------------------------------
129
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: UTriangle7P = &
130
        (/ 0.333333333333333D0, 0.101286507323456D0, 0.797426985353087D0, &
131
           0.101286507323456D0, 0.470142064105115D0, 0.059715871789770D0, &
132
           0.470142064105115D0 /)
133

134
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: VTriangle7P = &
135
        (/ 0.333333333333333D0, 0.101286507323456D0, 0.101286507323456D0, &
136
           0.797426985353087D0, 0.470142064105115D0, 0.470142064105115D0, &
137
           0.059715871789770D0 /)
138

139
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: STriangle7P = &
140
        (/ 0.225000000000000D0, 0.125939180544827D0, 0.125939180544827D0, &
141
           0.125939180544827D0, 0.132394152788506D0, 0.132394152788506D0, &
142
           0.132394152788506D0 /)
143

144
!------------------------------------------------------------------------------
145
! Triangle - 11 point rule; exact integration of x^py^q, p+q<=6
146
!------------------------------------------------------------------------------
147
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: UTriangle11P = &
148
    (/ 0.3019427231413448D-01, 0.5298143569082113D-01, &
149
       0.4972454892773975D-01, 0.7697772693248785D-01, &
150
       0.7008117469890058D+00, 0.5597774797709894D+00, &
151
       0.5428972301980696D+00, 0.3437947421925572D+00, &
152
       0.2356669356664465D+00, 0.8672623210691472D+00, &
153
       0.2151020995173866D+00 /)
154

155
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: VTriangle11P = &
156
    (/ 0.2559891985673773D+00, 0.1748087863744473D-01, &
157
       0.6330812033358987D+00, 0.8588528075577063D+00, &
158
       0.2708520519075563D+00, 0.1870602768957014D-01, &
159
       0.2027008533579804D+00, 0.5718576583152437D+00, &
160
       0.1000777578531811D+00, 0.4654861310422605D-01, &
161
       0.3929681357810497D+00 /)
162

163
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: STriangle11P = &
164
    (/ 0.3375321205342688D-01, 0.1148426034648707D-01, &
165
       0.4197958777582435D-01, 0.3098130358202468D-01, &
166
       0.2925899761167147D-01, 0.2778515729102349D-01, &
167
       0.8323049608963519D-01, 0.6825761580824108D-01, &
168
       0.6357334991651026D-01, 0.2649352562792455D-01, &
169
       0.8320249389723097D-01 /)
170

171
!------------------------------------------------------------------------------
172
! Triangle - 12 point rule; exact integration of x^py^q, p+q<=7
173
!------------------------------------------------------------------------------
174
   REAL(KIND=dp), DIMENSION(12), PRIVATE :: UTriangle12P = &
175
    (/ 0.6232720494911090D+00, 0.3215024938520235D+00, &
176
       0.5522545665686063D-01, 0.2777161669760336D+00, &
177
       0.5158423343535236D+00, 0.2064414986704435D+00, &
178
       0.3432430294535058D-01, 0.3047265008682535D+00, &
179
       0.6609491961864082D+00, 0.6238226509441210D-01, &
180
       0.8700998678316921D+00, 0.6751786707389329D-01 /)
181

182
   REAL(KIND=dp), DIMENSION(12), PRIVATE :: VTriangle12P = &
183
    (/ 0.3215024938520235D+00, 0.5522545665686063D-01, &
184
       0.6232720494911090D+00, 0.5158423343535236D+00, &
185
       0.2064414986704435D+00, 0.2777161669760336D+00, &
186
       0.3047265008682535D+00, 0.6609491961864082D+00, &
187
       0.3432430294535058D-01, 0.8700998678316921D+00, &
188
       0.6751786707389329D-01, 0.6238226509441210D-01 /)
189

190
   REAL(KIND=dp), DIMENSION(12), PRIVATE :: STriangle12P = &
191
    (/ 0.4388140871440586D-01, 0.4388140871440586D-01, &
192
       0.4388140871440587D-01, 0.6749318700971417D-01, &
193
       0.6749318700971417D-01, 0.6749318700971417D-01, &
194
       0.2877504278510970D-01, 0.2877504278510970D-01, &
195
       0.2877504278510969D-01, 0.2651702815743698D-01, &
196
       0.2651702815743698D-01, 0.2651702815743698D-01 /)
197

198
!------------------------------------------------------------------------------
199
! Triangle - 17 point rule; exact integration of x^py^q, p+q<=8
200
!------------------------------------------------------------------------------
201
   REAL(KIND=dp), DIMENSION(17), PRIVATE :: UTriangle17P = &
202
    (/ 0.2292423642627924D+00, 0.4951220175479885D-01, &
203
       0.3655948407066446D+00, 0.4364350639589269D+00, &
204
       0.1596405673569602D+00, 0.9336507149305228D+00, &
205
       0.5219569066777245D+00, 0.7110782758797098D+00, &
206
       0.5288509041694864D+00, 0.1396967677642513D-01, &
207
       0.4205421906708996D-01, 0.4651359156686354D-01, &
208
       0.1975981349257204D+00, 0.7836841874017514D+00, &
209
       0.4232808751402256D-01, 0.4557097415216423D+00, &
210
       0.2358934246935281D+00 /)
211

212
   REAL(KIND=dp), DIMENSION(17), PRIVATE :: VTriangle17P = &
213
    (/ 0.5117407211006358D+00, 0.7589103637479163D+00, &
214
       0.1529647481767193D+00, 0.3151398735074337D-01, &
215
       0.5117868393288316D-01, 0.1964516824106966D-01, &
216
       0.2347490459725670D+00, 0.4908682577187765D-01, &
217
       0.4382237537321878D+00, 0.3300210677033395D-01, &
218
       0.2088758614636060D+00, 0.9208929246654702D+00, &
219
       0.2742740954674795D+00, 0.1654585179097472D+00, &
220
       0.4930011699833554D+00, 0.4080804967846944D+00, &
221
       0.7127872162741824D+00 /)
222

223
   REAL(KIND=dp), DIMENSION(17), PRIVATE :: STriangle17P = &
224
    (/ 0.5956595662857148D-01, 0.2813390230006461D-01, &
225
       0.3500735477096827D-01, 0.2438077450393263D-01, &
226
       0.2843374448051010D-01, 0.7822856634218779D-02, &
227
       0.5179111341004783D-01, 0.3134229539096806D-01, &
228
       0.2454951584925144D-01, 0.5371382557647114D-02, &
229
       0.2571565514768072D-01, 0.1045933340802507D-01, &
230
       0.4937780841212319D-01, 0.2824772362317584D-01, &
231
       0.3218881684015661D-01, 0.2522089247693226D-01, &
232
       0.3239087356572598D-01 /)
233

234
!------------------------------------------------------------------------------
235
! Triangle - 20 point rule; exact integration of x^py^q, p+q<=9
236
!------------------------------------------------------------------------------
237
   REAL(KIND=dp), DIMENSION(20), PRIVATE :: UTriangle20P = &
238
   (/ 0.2469118866487856D-01, 0.3348782965514246D-00, &
239
      0.4162560937597861D-00, 0.1832492889417431D-00, &
240
      0.2183952668281443D-00, 0.4523362527443628D-01, &
241
      0.4872975112073226D-00, 0.7470127381316580D-00, &
242
      0.7390287107520658D-00, 0.3452260444515281D-01, &
243
      0.4946745467572288D-00, 0.3747439678780460D-01, &
244
      0.2257524791528391D-00, 0.9107964437563798D-00, &
245
      0.4254445629399445D-00, 0.1332215072275240D-00, &
246
      0.5002480151788234D-00, 0.4411517722238269D-01, &
247
      0.1858526744057914D-00, 0.6300024376672695D-00 /)
248

249
   REAL(KIND=dp), DIMENSION(20), PRIVATE :: VTriangle20P = &
250
   (/ 0.4783451248176442D-00, 0.3373844236168988D-00, &
251
      0.1244378463254732D-00, 0.6365569723648120D-00, &
252
      0.3899759363237886D-01, 0.9093437140096456D-00, &
253
      0.1968266037596590D-01, 0.2191311129347709D-00, &
254
      0.3833588560240875D-01, 0.7389795063475102D-00, &
255
      0.4800989285800525D-00, 0.2175137165318176D-00, &
256
      0.7404716820879975D-00, 0.4413531509926682D-01, &
257
      0.4431292142978816D-00, 0.4440953593652837D-00, &
258
      0.1430831401051367D-00, 0.4392970158411878D-01, &
259
      0.1973209364545017D-00, 0.1979381059170009D-00 /)
260

261
   REAL(KIND=dp), DIMENSION(20), PRIVATE :: STriangle20P = &
262
   (/ 0.1776913091122958D-01, 0.4667544936904065D-01, &
263
      0.2965283331432967D-01, 0.3880447634997608D-01, &
264
      0.2251511457011248D-01, 0.1314162394636178D-01, &
265
      0.1560341736610505D-01, 0.1967065434689744D-01, &
266
      0.2247962849501080D-01, 0.2087108394969067D-01, &
267
      0.1787661200700672D-01, 0.2147695865607915D-01, &
268
      0.2040998247303970D-01, 0.1270342300533680D-01, &
269
      0.3688713099356314D-01, 0.3813199811535777D-01, &
270
      0.1508642325812160D-01, 0.1238422287692121D-01, &
271
      0.3995072336992735D-01, 0.3790911262589247D-01 /)
272

273
!------------------------------------------------------------------------------
274
! Tetrahedron - 1 point rule; exact integration of x^py^qz^r, p+q+r<=1
275
!------------------------------------------------------------------------------
276
   REAL(KIND=dp), DIMENSION(1), PRIVATE :: UTetra1P = (/ 0.25D0 /)
277
   REAL(KIND=dp), DIMENSION(1), PRIVATE :: VTetra1P = (/ 0.25D0 /)
278
   REAL(KIND=dp), DIMENSION(1), PRIVATE :: WTetra1P = (/ 0.25D0 /)
279
   REAL(KIND=dp), DIMENSION(1), PRIVATE :: STetra1P = (/ 1.00D0 /)
280

281
!------------------------------------------------------------------------------
282
! Tetrahedron - 4 point rule; exact integration of x^py^qz^r, p+q+r<=2
283
!------------------------------------------------------------------------------
284
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: UTetra4P = &
285
    (/ 0.1757281246520584D0, 0.2445310270213291D0, &
286
       0.5556470949048655D0, 0.0240937534217468D0 /)
287

288
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: VTetra4P = &
289
    (/ 0.5656137776620919D0, 0.0501800797762026D0, &
290
       0.1487681308666864D0, 0.2354380116950194D0 /)
291

292
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: WTetra4P = &
293
    (/ 0.2180665126782654D0, 0.5635595064952189D0, &
294
       0.0350112499848832D0, 0.1833627308416330D0 /)
295

296
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: STetra4P = &
297
    (/ 0.2500000000000000D0, 0.2500000000000000D0, &
298
       0.2500000000000000D0, 0.2500000000000000D0 /)
299
!------------------------------------------------------------------------------
300
! Tetrahedron - 5 point rule; exact integration of x^py^qz^r, p+q+r<=3
301
!------------------------------------------------------------------------------
302
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: UTetra5P =  &
303
    (/ 0.25000000000000000D0, 0.50000000000000000D0, &
304
       0.16666666666666667D0, 0.16666666666666667D0, &
305
       0.16666666666666667D0 /)
306

307
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: VTetra5P =  &
308
    (/ 0.25000000000000000D0, 0.16666666666666667D0, &
309
       0.50000000000000000D0, 0.16666666666666667D0, &
310
       0.16666666666666667D0 /)
311

312
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: WTetra5P =  &
313
    (/ 0.25000000000000000D0, 0.16666666666666667D0, &
314
       0.16666666666666667D0, 0.50000000000000000D0, &
315
       0.16666666666666667D0 /)
316

317
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: STetra5P =  &
318
    (/-0.80000000000000000D0, 0.45000000000000000D0, &
319
       0.45000000000000000D0, 0.45000000000000000D0, &
320
       0.45000000000000000D0 /)
321
!------------------------------------------------------------------------------
322
! Tetrahedron - 11 point rule; exact integration of x^py^qz^r, p+q+r<=4
323
!------------------------------------------------------------------------------
324
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: UTetra11P =  &
325
    (/ 0.3247902050850455D+00, 0.4381969657060433D+00, &
326
       0.8992592373310454D-01, 0.1092714936292849D+00, &
327
       0.3389119319942253D-01, 0.5332363613904868D-01, &
328
       0.1935618747806815D+00, 0.4016250624424964D-01, &
329
       0.3878132182319405D+00, 0.7321489692875428D+00, &
330
       0.8066342495294049D-01 /)
331

332
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: VTetra11P =  &
333
    (/ 0.4573830181783998D+00, 0.9635325047480842D-01, &
334
       0.3499588148445295D+00, 0.1228957438582778D+00, &
335
       0.4736224692062527D-01, 0.4450376952468180D+00, &
336
       0.2165626476982170D+00, 0.8033385922433729D+00, &
337
       0.7030897281814283D-01, 0.1097836536360084D+00, &
338
       0.1018859284267242D-01 /)
339

340
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: WTetra11P =  &
341
    (/ 0.1116787541193331D+00, 0.6966288385119494D-01, &
342
       0.5810783971325720D-01, 0.3424607753785182D+00, &
343
       0.7831772466208499D+00, 0.3688112094344830D+00, &
344
       0.5872345323698884D+00, 0.6178518963560731D-01, &
345
       0.4077342860913465D+00, 0.9607290317342082D-01, &
346
       0.8343823045787845D-01 /)
347

348
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: STetra11P =  &
349
    (/ 0.1677896627448221D+00, 0.1128697325878004D+00, &
350
       0.1026246621329828D+00, 0.1583002576888426D+00, &
351
       0.3847841737508437D-01, 0.1061709382037234D+00, &
352
       0.5458124994014422D-01, 0.3684475128738168D-01, &
353
       0.1239234851349682D+00, 0.6832098141300300D-01, &
354
       0.3009586149124714D-01 /)
355
!------------------------------------------------------------------------------
356

357
!------------------------------------------------------------------------------
358
! A SELECTION OF QUADRATURE RULES UP TO A HIGH ORDER DESIGNED FOR ECONOMIC
359
! INTEGRATION OF COMPLETE POLYNOMIALS
360
!
361
! These rules have been reproduced from the source "P. Solin, K. Segeth
362
! and I. Dolezel: Higher-Order Finite Element Methods",
363
! Chapman & Hall/CRC Press, 2003.
364
!------------------------------------------------------------------------------
365

366
!------------------------------------------------------------------------------
367
! Quadrilateral - 8-point rule for complete polynomials of order p<=5
368
!------------------------------------------------------------------------------
369
   REAL(KIND=dp), DIMENSION(8), PRIVATE :: UPQuad8 = &
370
       (/ 0.683130051063973d0,-0.683130051063973d0, 0.000000000000000d0, &
371
          0.000000000000000d0, 0.881917103688197d0, 0.881917103688197d0, &
372
          -0.881917103688197d0,-0.881917103688197d0 /)
373

374
   REAL(KIND=dp), DIMENSION(8), PRIVATE :: VPQuad8 = &
375
       (/ 0.000000000000000d0, 0.000000000000000d0, 0.683130051063973d0, &
376
         -0.683130051063973d0, 0.881917103688197d0,-0.881917103688197d0, &
377
          0.881917103688197d0,-0.881917103688197d0 /)
378

379
   REAL(KIND=dp), DIMENSION(8), PRIVATE :: SPQuad8 = &
380
       (/ 0.816326530612245d0, 0.816326530612245d0, 0.816326530612245d0, &
381
          0.816326530612245d0, 0.183673469387755d0, 0.183673469387755d0, &
382
          0.183673469387755d0, 0.183673469387755d0 /)
383
!------------------------------------------------------------------------------
384
! Quadrilateral - 12-point rule for complete polynomials of order p<=7
385
! NOTE: It seems that this rule may give somehow faulty results at least
386
!       for some serendipity polynomials of degree 3. Therefore
387
!       we never use it. 
388
!------------------------------------------------------------------------------
389
   REAL(KIND=dp), DIMENSION(12), PRIVATE :: UPQuad12 = &
390
       (/ 0.925820099772551d0, -0.925820099772551d0,  0.000000000000000d0, &
391
          0.000000000000000d0,  0.805979782918599d0,  0.805979782918599d0, &
392
         -0.805979782918599d0, -0.805979782918599d0,  0.380554433208316d0, &
393
          0.380554433208316d0, -0.380554433208316d0, -0.380554433208316d0 /)
394

395
   REAL(KIND=dp), DIMENSION(12), PRIVATE :: VPQuad12 = &
396
       (/ 0.000000000000000d0,  0.000000000000000d0,  0.925820099772551d0, &
397
         -0.925820099772551d0,  0.805979782918599d0, -0.805979782918599d0, &
398
          0.805979782918599d0, -0.805979782918599d0,  0.380554433208316d0, &
399
         -0.380554433208316d0,  0.380554433208316d0, -0.380554433208316d0 /)
400

401
   REAL(KIND=dp), DIMENSION(12), PRIVATE :: SPQuad12 = &
402
       (/ 0.241975308641975d0, 0.241975308641975d0, 0.241975308641975d0, &
403
          0.241975308641975d0, 0.237431774690630d0, 0.237431774690630d0, &
404
          0.237431774690630d0, 0.237431774690630d0, 0.520592916667394d0, &
405
          0.520592916667394d0, 0.520592916667394d0, 0.520592916667394d0 /)
406
!------------------------------------------------------------------------------
407
! Quadrilateral - 20-point rule for complete polynomials of order p<=9
408
! Note: Some points lie outside the reference element
409
!------------------------------------------------------------------------------
410
   REAL(KIND=dp), DIMENSION(20), PRIVATE :: UPQuad20 = &
411
       (/ 0.112122576386656d1,  -0.112122576386656d1,   0.000000000000000d0,  &
412
          0.000000000000000d0,   0.451773049920657d0,  -0.451773049920657d0,  &
413
          0.000000000000000d0,   0.000000000000000d0,   0.891849420851512d0,  &
414
          0.891849420851512d0,  -0.891849420851512d0,  -0.891849420851512d0,  &
415
          0.824396370749276d0,   0.824396370749276d0,  -0.824396370749276d0,  &
416
         -0.824396370749276d0,   0.411623426336542d0,   0.411623426336542d0,  &
417
         -0.411623426336542d0,  -0.411623426336542d0 /)
418

419
   REAL(KIND=dp), DIMENSION(20), PRIVATE :: VPQuad20 = &
420
       (/ 0.000000000000000d0,  0.000000000000000d0,  0.112122576386656d1, &
421
         -0.112122576386656d1,  0.000000000000000d0,  0.000000000000000d0, &
422
          0.451773049920657d0, -0.451773049920657d0,  0.891849420851512d0, &
423
         -0.891849420851512d0,  0.891849420851512d0, -0.891849420851512d0, &
424
          0.411623426336542d0, -0.411623426336542d0,  0.411623426336542d0, &
425
         -0.411623426336542d0,  0.824396370749276d0, -0.824396370749276d0, &
426
         0.824396370749276d0, -0.824396370749276d0 /)
427

428
   REAL(KIND=dp), DIMENSION(20), PRIVATE :: SPQuad20 = &
429
       (/ 0.184758425074910d-1,  0.184758425074910d-1,  0.184758425074910d-1,  &
430
          0.184758425074910d-1,  0.390052939160735d0,   0.390052939160735d0,   &
431
          0.390052939160735d0,   0.390052939160735d0,   0.830951780264820d-1,  &
432
          0.830951780264820d-1,  0.830951780264820d-1,  0.830951780264820d-1,  &
433
          0.254188020152646d0,   0.254188020152646d0,   0.254188020152646d0,   &
434
          0.254188020152646d0,   0.254188020152646d0,   0.254188020152646d0,   &
435
          0.254188020152646d0,   0.254188020152646d0 /)
436
!------------------------------------------------------------------------------
437
! Quadrilateral - 25-point rule for complete polynomials of order p<=11
438
! Note: Some points lie outside the reference element
439
!------------------------------------------------------------------------------
440
   REAL(KIND=dp), DIMENSION(25), PRIVATE :: UPQuad25 = &
441
       (/0.000000000000000d0,   0.104440291540981d1,  -0.104440291540981d1,  &
442
         0.000000000000000d0,   0.000000000000000d0,   0.769799068396649d0,  &
443
        -0.769799068396649d0,   0.000000000000000d0,   0.000000000000000d0,  &
444
         0.935787012440540d0,   0.935787012440540d0,  -0.935787012440540d0,  &
445
        -0.935787012440540d0,   0.413491953449114d0,   0.413491953449114d0,  &
446
        -0.413491953449114d0,  -0.413491953449114d0,   0.883025508525690d0,  &
447
         0.883025508525690d0,  -0.883025508525690d0,  -0.883025508525690d0,  &
448
         0.575653595840465d0,   0.575653595840465d0,  -0.575653595840465d0,  &
449
        -0.575653595840465d0 /)
450

451
   REAL(KIND=dp), DIMENSION(25), PRIVATE :: VPQuad25 = &
452
       (/ 0.000000000000000d0,   0.000000000000000d0,   0.000000000000000d0, &
453
          0.104440291540981d1,  -0.104440291540981d1,   0.000000000000000d0, &
454
          0.000000000000000d0,   0.769799068396649d0,  -0.769799068396649d0, &
455
          0.935787012440540d0,  -0.935787012440540d0,   0.935787012440540d0, &
456
         -0.935787012440540d0,   0.413491953449114d0,  -0.413491953449114d0, &
457
          0.413491953449114d0,  -0.413491953449114d0,   0.575653595840465d0, &
458
         -0.575653595840465d0,   0.575653595840465d0,  -0.575653595840465d0, &
459
          0.883025508525690d0,  -0.883025508525690d0,   0.883025508525690d0, &
460
          -0.883025508525690d0 /)
461

462
   REAL(KIND=dp), DIMENSION(25), PRIVATE :: SPQuad25 = &
463
       (/ 0.365379525585903d0,   0.277561655642040d-1,   0.277561655642040d-1, &
464
          0.277561655642040d-1,  0.277561655642040d-1,   0.244272057751754d0,  &
465
          0.244272057751754d0,   0.244272057751754d0,    0.244272057751754d0,  &
466
          0.342651038512290d-1,  0.342651038512290d-1,   0.342651038512290d-1, &
467
          0.342651038512290d-1,  0.308993036133713d0,    0.308993036133713d0,  &
468
          0.308993036133713d0,   0.308993036133713d0,    0.146684377651312d0,  &
469
          0.146684377651312d0,   0.146684377651312d0,    0.146684377651312d0,  &
470
          0.146684377651312d0,   0.146684377651312d0,    0.146684377651312d0,  &
471
          0.146684377651312d0  /)
472
!------------------------------------------------------------------------------
473
! Quadrilateral - 36-point rule for complete polynomials of order p<=13
474
! Note 1: Some points lie outside the reference element
475
! Note 2: Weights somewhat inaccurate?
476
!------------------------------------------------------------------------------
477
   REAL(KIND=dp), DIMENSION(36), PRIVATE :: UPQuad36 = &
478
       (/  0.108605615857397d1,  -0.108605615857397d1,   0.000000000000000d0, &
479
        0.000000000000000d0,   0.658208197042585d0,  -0.658208197042585d0,    &
480
        0.000000000000000d0,   0.000000000000000d0,   0.100130060299173d1,    &
481
        0.100130060299173d1,  -0.100130060299173d1,  -0.100130060299173d1,    &
482
        0.584636168775946d0,   0.584636168775946d0,  -0.584636168775946d0,    &
483
       -0.584636168775946d0,   0.246795612720261d0,   0.246795612720261d0,    &
484
       -0.246795612720261d0,  -0.246795612720261d0,   0.900258815287201d0,    &
485
        0.900258815287201d0,  -0.900258815287201d0,  -0.900258815287201d0,    &
486
        0.304720678579870d0,   0.304720678579870d0,  -0.304720678579870d0,    &
487
       -0.304720678579870d0,   0.929866705560780d0,   0.929866705560780d0,    &
488
       -0.929866705560780d0,  -0.929866705560780d0,   0.745052720131169d0,    &
489
        0.745052720131169d0,  -0.745052720131169d0,  -0.745052720131169d0  /)
490

491
   REAL(KIND=dp), DIMENSION(36), PRIVATE :: VPQuad36 = &
492
       (/ 0.000000000000000d0,   0.000000000000000d0,   0.108605615857397d1,  &
493
       -0.108605615857397d1,   0.000000000000000d0,   0.000000000000000d0, &
494
        0.658208197042585d0,  -0.658208197042585d0,   0.100130060299173d1, &
495
       -0.100130060299173d1,   0.100130060299173d1,  -0.100130060299173d1, &
496
        0.584636168775946d0,  -0.584636168775946d0,   0.584636168775946d0, &
497
       -0.584636168775946d0,   0.246795612720261d0,  -0.246795612720261d0, &
498
        0.246795612720261d0,  -0.246795612720261d0,   0.304720678579870d0, &
499
       -0.304720678579870d0,   0.304720678579870d0,  -0.304720678579870d0, &
500
        0.900258815287201d0,  -0.900258815287201d0,   0.900258815287201d0, &
501
       -0.900258815287201d0,   0.745052720131169d0,  -0.745052720131169d0, &
502
        0.745052720131169d0,  -0.745052720131169d0,   0.929866705560780d0, &
503
       -0.929866705560780d0,   0.929866705560780d0,  -0.929866705560780d0 /)
504

505
   REAL(KIND=dp), DIMENSION(36), PRIVATE :: SPQuad36 = &
506
       (/ 0.565616969376400d-2,   0.565616969376400d-2,   0.565616969376400d-2,   &
507
       0.565616969376400d-2,  0.192443967470396d0,   0.192443967470396d0,  &
508
       0.192443967470396d0,   0.192443967470396d0,   0.516683297977300d-2, &
509
       0.516683297977300d-2,  0.516683297977300d-2,  0.516683297977300d-2, &
510
       0.200302559622138d0,   0.200302559622138d0,   0.200302559622138d0,  &
511
       0.200302559622138d0,   0.228125175912536d0,   0.228125175912536d0,  &
512
       0.228125175912536d0,   0.228125175912536d0,   0.117496926974491d0,  &
513
       0.117496926974491d0,   0.117496926974491d0,   0.117496926974491d0,  &
514
       0.117496926974491d0,   0.117496926974491d0,   0.117496926974491d0,  &
515
       0.117496926974491d0,   0.666557701862050d-1,  0.666557701862050d-1, &
516
       0.666557701862050d-1,  0.666557701862050d-1,  0.666557701862050d-1, &
517
       0.666557701862050d-1,  0.666557701862050d-1,  0.666557701862050d-1  /)
518
!------------------------------------------------------------------------------
519
! Quadrilateral - 45-point rule for complete polynomials of order p<=15
520
! Note: Some points lie outside the reference element
521
!------------------------------------------------------------------------------
522
   REAL(KIND=dp), DIMENSION(45), PRIVATE :: UPQuad45 = &
523
       (/ 0.000000000000000d0,   0.102731435771937d1,  -0.102731435771937d1,   &
524
        0.000000000000000d0,   0.000000000000000d0,   0.856766776147643d0,  &
525
       -0.856766776147643d0,   0.000000000000000d0,   0.000000000000000d0,  &
526
        0.327332998189723d0,  -0.327332998189723d0,   0.000000000000000d0,  &
527
        0.000000000000000d0,   0.967223740028505d0,   0.967223740028505d0,  &
528
       -0.967223740028505d0,  -0.967223740028505d0,   0.732168901749711d0,  &
529
        0.732168901749711d0,  -0.732168901749711d0,  -0.732168901749711d0,  &
530
        0.621974427996805d0,   0.621974427996805d0,  -0.621974427996805d0,  &
531
       -0.621974427996805d0,   0.321696694921009d0,   0.321696694921009d0,  &
532
       -0.321696694921009d0,  -0.321696694921009d0,   0.928618480068352d0,  &
533
        0.928618480068352d0,  -0.928618480068352d0,  -0.928618480068352d0,  &
534
        0.455124178121179d0,   0.455124178121179d0,  -0.455124178121179d0,  &
535
       -0.455124178121179d0,   0.960457474887516d0,   0.960457474887516d0,  &
536
       -0.960457474887516d0,  -0.960457474887516d0,   0.809863684081217d0,  &
537
        0.809863684081217d0,  -0.809863684081217d0,  -0.809863684081217d0   /)
538

539
   REAL(KIND=dp), DIMENSION(45), PRIVATE :: VPQuad45 = &
540
       (/ 0.000000000000000d0,   0.000000000000000d0,   0.000000000000000d0,   &
541
        0.102731435771937d1,  -0.102731435771937d1,   0.000000000000000d0, &
542
        0.000000000000000d0,   0.856766776147643d0,  -0.856766776147643d0, &
543
        0.000000000000000d0,   0.000000000000000d0,   0.327332998189723d0, &
544
       -0.327332998189723d0,   0.967223740028505d0,  -0.967223740028505d0, &
545
        0.967223740028505d0,  -0.967223740028505d0,   0.732168901749711d0, &
546
       -0.732168901749711d0,   0.732168901749711d0,  -0.732168901749711d0, &
547
        0.321696694921009d0,  -0.321696694921009d0,   0.321696694921009d0, &
548
       -0.321696694921009d0,   0.621974427996805d0,  -0.621974427996805d0, &
549
        0.621974427996805d0,  -0.621974427996805d0,   0.455124178121179d0, &
550
       -0.455124178121179d0,   0.455124178121179d0,  -0.455124178121179d0, &
551
        0.928618480068352d0,  -0.928618480068352d0,   0.928618480068352d0, &
552
       -0.928618480068352d0,   0.809863684081217d0,  -0.809863684081217d0, &
553
        0.809863684081217d0,  -0.809863684081217d0,   0.960457474887516d0, &
554
       -0.960457474887516d0,   0.960457474887516d0,  -0.960457474887516d0 /)
555

556
   REAL(KIND=dp), DIMENSION(45), PRIVATE :: SPQuad45 = &
557
       (/ -0.176897982720700d-2,   0.128167266175120d-1,   0.128167266175120d-1,   &
558
       0.128167266175120d-1,  0.128167266175120d-1,  0.119897873101347d0,   &
559
       0.119897873101347d0,   0.119897873101347d0,   0.119897873101347d0,   &
560
       0.210885452208801d0,   0.210885452208801d0,   0.210885452208801d0,   &
561
       0.210885452208801d0,   0.639272012821500d-2,  0.639272012821500d-2,  &
562
       0.639272012821500d-2,  0.639272012821500d-2,  0.104415680788580d0,   &
563
       0.104415680788580d0,   0.104415680788580d0,   0.104415680788580d0,   &
564
       0.168053047203816d0,   0.168053047203816d0,   0.168053047203816d0,   &
565
       0.168053047203816d0,   0.168053047203816d0,   0.168053047203816d0,   &
566
       0.168053047203816d0,   0.168053047203816d0,   0.761696944522940d-1,  &
567
       0.761696944522940d-1,  0.761696944522940d-1,  0.761696944522940d-1,  &
568
       0.761696944522940d-1,  0.761696944522940d-1,  0.761696944522940d-1,  &
569
       0.761696944522940d-1,  0.287941544000640d-1,  0.287941544000640d-1,  &
570
       0.287941544000640d-1,  0.287941544000640d-1,  0.287941544000640d-1,  &
571
       0.287941544000640d-1,  0.287941544000640d-1,  0.287941544000640d-1 /)
572
!------------------------------------------------------------------------------
573
! Quadrilateral - 60-point rule for complete polynomials of order p<=17
574
!------------------------------------------------------------------------------
575
   REAL(KIND=dp), DIMENSION(60), PRIVATE :: UPQuad60 = &
576
       (/ 0.989353074512600d0,  -0.989353074512600d0,   0.000000000000000d0,   &
577
        0.000000000000000d0,   0.376285207157973d0,  -0.376285207157973d0,  &
578
        0.000000000000000d0,   0.000000000000000d0,   0.978848279262233d0,  &
579
        0.978848279262233d0,  -0.978848279262233d0,  -0.978848279262233d0,  &
580
        0.885794729164116d0,   0.885794729164116d0,  -0.885794729164116d0,  &
581
       -0.885794729164116d0,   0.171756123838348d0,   0.171756123838348d0,  &
582
       -0.171756123838348d0,  -0.171756123838348d0,   0.590499273806002d0,  &
583
        0.590499273806002d0,  -0.590499273806002d0,  -0.590499273806002d0,  &
584
        0.319505036634574d0,   0.319505036634574d0,  -0.319505036634574d0,  &
585
       -0.319505036634574d0,   0.799079131916863d0,   0.799079131916863d0,  &
586
       -0.799079131916863d0,  -0.799079131916863d0,   0.597972451929457d0,  &
587
        0.597972451929457d0,  -0.597972451929457d0,  -0.597972451929457d0,  &
588
        0.803743962958745d0,   0.803743962958745d0,  -0.803743962958745d0,  &
589
       -0.803743962958745d0,   0.583444817765510d-1,  0.583444817765510d-1, &
590
       -0.583444817765510d-1, -0.583444817765510d-1,  0.936506276127495d0,  &
591
        0.936506276127495d0,  -0.936506276127495d0,  -0.936506276127495d0,  &
592
        0.347386316166203d0,   0.347386316166203d0,  -0.347386316166203d0,  &
593
       -0.347386316166203d0,   0.981321179805452d0,   0.981321179805452d0,  &
594
       -0.981321179805452d0,  -0.981321179805452d0,   0.706000287798646d0,  &
595
        0.706000287798646d0,  -0.706000287798646d0,  -0.706000287798646d0 /)
596

597
   REAL(KIND=dp), DIMENSION(60), PRIVATE :: VPQuad60 = &
598
       (/0.000000000000000d0,   0.000000000000000d0,   0.989353074512600d0,  &
599
       -0.989353074512600d0,   0.000000000000000d0,   0.000000000000000d0,  &
600
        0.376285207157973d0,  -0.376285207157973d0,   0.978848279262233d0,  &
601
       -0.978848279262233d0,   0.978848279262233d0,  -0.978848279262233d0,  &
602
        0.885794729164116d0,  -0.885794729164116d0,   0.885794729164116d0,  &
603
       -0.885794729164116d0,   0.171756123838348d0,  -0.171756123838348d0,  &
604
        0.171756123838348d0,  -0.171756123838348d0,   0.319505036634574d0,  &
605
       -0.319505036634574d0,   0.319505036634574d0,  -0.319505036634574d0,  &
606
        0.590499273806002d0,  -0.590499273806002d0,   0.590499273806002d0,  &
607
       -0.590499273806002d0,   0.597972451929457d0,  -0.597972451929457d0,  &
608
        0.597972451929457d0,  -0.597972451929457d0,   0.799079131916863d0,  &
609
       -0.799079131916863d0,   0.799079131916863d0,  -0.799079131916863d0,  &
610
        0.583444817765510d-1, -0.583444817765510d-1,  0.583444817765510d-1, &
611
       -0.583444817765510d-1,  0.803743962958745d0,  -0.803743962958745d0,  &
612
        0.803743962958745d0,  -0.803743962958745d0,   0.347386316166203d0,  &
613
       -0.347386316166203d0,   0.347386316166203d0,  -0.347386316166203d0,  &
614
        0.936506276127495d0,  -0.936506276127495d0,   0.936506276127495d0,  &
615
       -0.936506276127495d0,   0.706000287798646d0,  -0.706000287798646d0,  &
616
        0.706000287798646d0,  -0.706000287798646d0,   0.981321179805452d0,  &
617
       -0.981321179805452d0,   0.981321179805452d0,  -0.981321179805452d0  /)
618

619
   REAL(KIND=dp), DIMENSION(60), PRIVATE :: SPQuad60 = &
620
       (/ 0.206149159199910d-1,   0.206149159199910d-1,   0.206149159199910d-1,   &
621
       0.206149159199910d-1,  0.128025716179910d0,    0.128025716179910d0,  &
622
       0.128025716179910d0,   0.128025716179910d0,    0.551173953403200d-2, &
623
       0.551173953403200d-2,  0.551173953403200d-2,   0.551173953403200d-2, &
624
       0.392077124571420d-1,  0.392077124571420d-1,   0.392077124571420d-1, &
625
       0.392077124571420d-1,  0.763969450798630d-1,   0.763969450798630d-1, &
626
       0.763969450798630d-1,  0.763969450798630d-1,   0.141513729949972d0,  &
627
       0.141513729949972d0,   0.141513729949972d0,    0.141513729949972d0,  &
628
       0.141513729949972d0,   0.141513729949972d0,    0.141513729949972d0,  &
629
       0.141513729949972d0,   0.839032793637980d-1,   0.839032793637980d-1, &
630
       0.839032793637980d-1,  0.839032793637980d-1,   0.839032793637980d-1, &
631
       0.839032793637980d-1,  0.839032793637980d-1,   0.839032793637980d-1, &
632
       0.603941636496850d-1,  0.603941636496850d-1,   0.603941636496850d-1, &
633
       0.603941636496850d-1,  0.603941636496850d-1,   0.603941636496850d-1, &
634
       0.603941636496850d-1,  0.603941636496850d-1,   0.573877529692130d-1, &
635
       0.573877529692130d-1,  0.573877529692130d-1,   0.573877529692130d-1, &
636
       0.573877529692130d-1,  0.573877529692130d-1,   0.573877529692130d-1, &
637
       0.573877529692130d-1,  0.219225594818640d-1,   0.219225594818640d-1, &
638
       0.219225594818640d-1,  0.219225594818640d-1,   0.219225594818640d-1, &
639
       0.219225594818640d-1,  0.219225594818640d-1,   0.219225594818640d-1  /)
640

641

642
!------------------------------------------------------------------------------
643
! These rules have been reproduced from the paper
644
! Ethan J. Kubatko, Benjamin A. Yeager, Ashley L. Magg and I. Dolezel:
645
! "New computationally efficient quadrature formulas for triangular
646
! prism elements.", Computers & Fluids 73 (2013) 187–201.
647
!------------------------------------------------------------------------------
648

649
!------------------------------------------------------------------------------
650
! Wedge - 4 point rule, 2nd order polynoms
651
!------------------------------------------------------------------------------
652
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: SWedge4P = [ &
653
       1.000000000000000d0, &
654
       1.000000000000000d0, &
655
       1.000000000000000d0, &
656
       1.000000000000000d0]
657
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: UWedge4P = [ &
658
       0.483163247594393d0, &
659
       -0.605498860309242d0, &
660
       -0.605498860309242d0, &
661
       -0.605498860309242d0]
662
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: VWedge4P = [ &
663
       -0.741581623797196d0, &
664
       0.469416096821288d0, &
665
       -0.530583903178712d0, &
666
       -0.530583903178712d0]
667
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: WWedge4P = [ &
668
       0.000000000000000d0, &
669
       0.000000000000000d0, &
670
       0.816496580927726d0, &
671
       -0.816496580927726d0 ]
672

673
!------------------------------------------------------------------------------
674
! Wedge - n=5, p=2
675
!------------------------------------------------------------------------------
676
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: SWedge5P = [ &
677
       0.333333333333333d0, &
678
       0.333333333333333d0, &
679
       0.333333333333333d0, &
680
       1.500000000000000d0, &
681
       1.500000000000000d0 ]
682
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: UWedge5P = [ &
683
       -1.000000000000000d0, &
684
       -1.000000000000000d0, &
685
       1.000000000000000d0, &
686
       -0.333333333333333d0, &
687
       -0.333333333333333d0 ]
688
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: VWedge5P = [ &
689
       -1.000000000000000d0, &
690
       1.000000000000000d0, &
691
       -1.000000000000000d0, &
692
       -0.333333333333333d0, &
693
       -0.333333333333333d0 ]
694
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: WWedge5P = [ &
695
       0.000000000000000d0, &
696
       0.000000000000000d0, &
697
       0.000000000000000d0, &
698
       0.666666666666667d0, &
699
       -0.666666666666667d0 ]
700

701
!------------------------------------------------------------------------------
702
! Wedge - n=6, p=3
703
!------------------------------------------------------------------------------
704
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: SWedge6P = [ &
705
       0.742534890852309d0, &
706
       0.375143463443327d0, &
707
       0.495419047908462d0, &
708
       0.523999970843238d0, &
709
       0.980905839025611d0, &
710
       0.881996787927053d0 ]
711
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: UWedge6P = [ &
712
       0.240692796349625d0, &
713
       -0.968326281451138d0, &
714
       0.467917833640195d0, &
715
       -0.786144119530819d0, &
716
       -0.484844897886675d0, &
717
       -0.559053711932125d0 ]
718
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: VWedge6P = [ &
719
       -0.771991660873412d0, &
720
       -0.568046512457875d0, &
721
       -0.549342790168347d0, &
722
       0.362655041695561d0, &
723
       -0.707931130717342d0, &
724
       0.260243325246813d0 ]
725
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: WWedge6P = [ &
726
       0.614747128207527d0, &
727
       0.676689529541421d0, &
728
       -0.599905857322635d0, &
729
       -0.677609795694786d0, &
730
       -0.502482717716373d0, &
731
       0.493010512161538d0 ]
732

733
!------------------------------------------------------------------------------
734
! Wedge - n=7, p=3
735
!------------------------------------------------------------------------------
736
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: SWedge7P = [ &
737
       -2.250000000000000d0, &
738
       1.041666666666667d0, &
739
       1.041666666666667d0, &
740
       1.041666666666667d0, &
741
       1.041666666666667d0, &
742
       1.041666666666667d0, &
743
       1.041666666666667d0 ]
744
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: UWedge7P = [ &
745
       -0.333333333333333d0, &
746
       -0.600000000000000d0, &
747
       -0.600000000000000d0, &
748
       0.200000000000000d0, &
749
       -0.600000000000000d0, &
750
       -0.600000000000000d0, &
751
       0.200000000000000d0 ]
752
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: VWedge7P = [ &
753
       -0.333333333333333d0, &
754
       -0.600000000000000d0, &
755
       0.200000000000000d0, &
756
       -0.600000000000000d0, &
757
       -0.600000000000000d0, &
758
       0.200000000000000d0, &
759
       -0.600000000000000d0 ]
760
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: WWedge7P = [ &
761
       0.000000000000000d0, &
762
       0.461880215351701d0, &
763
       0.461880215351701d0, &
764
       0.461880215351701d0, &
765
       -0.461880215351701d0, &
766
       -0.461880215351701d0, &
767
       -0.461880215351701d0 ]
768

769
!------------------------------------------------------------------------------
770
! Wedge - n=10, p=4
771
!------------------------------------------------------------------------------
772
   REAL(KIND=dp), DIMENSION(10), PRIVATE :: SWedge10P = [ &
773
       0.111155943811228d0, &
774
       0.309060899887509d0, &
775
       0.516646862442958d0, &
776
       0.567975205132714d0, &
777
       0.382742555939017d0, &
778
       0.355960928492268d0, &
779
       0.108183228294342d0, &
780
       0.126355242780924d0, &
781
       0.587370828592853d0, &
782
       0.934548304626188d0 ]
783
   REAL(KIND=dp), DIMENSION(10), PRIVATE :: UWedge10P = [ &
784
       0.812075900047562d0, &
785
       -0.792166223585545d0, &
786
       -0.756726179789306d0, &
787
       -0.552495167978340d0, &
788
       -0.357230019521233d0, &
789
       -0.987225392999058d0, &
790
       -0.816603728785918d0, &
791
       0.423489172633859d0, &
792
       0.363041084609230d0, &
793
       -0.175780343149613d0 ]
794
   REAL(KIND=dp), DIMENSION(10), PRIVATE :: VWedge10P = [ &
795
       -0.986242751499303d0, &
796
       0.687201105597868d0, &
797
       -0.731311840596107d0, &
798
       0.015073398439985d0, &
799
       0.126888850505978d0, &
800
       0.082647545710800d0, &
801
       -0.915066171481315d0, &
802
       -1.112801167237130d0, &
803
       -0.499011410082669d0, &
804
       -0.654971142379686d0 ]
805
   REAL(KIND=dp), DIMENSION(10), PRIVATE :: WWedge10P = [ &
806
       0.850716248413834d0, &
807
       -0.115214772515700d0, &
808
       -0.451491675441927d0, &
809
       -0.824457000064439d0, &
810
       0.855349689995606d0, &
811
       0.452976444667786d0, &
812
       0.997939285245240d0, &
813
       -0.963298774205756d0, &
814
       -0.299892769705443d0, &
815
       0.367947041936472d0 ]
816

817
!------------------------------------------------------------------------------
818
! Wedge - n=11, p=4
819
!------------------------------------------------------------------------------
820
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: SWedge11P = [ &
821
       0.545658450421913d0, &
822
       0.545658450421913d0, &
823
       0.545658450421913d0, &
824
       0.431647899262139d0, &
825
       0.249954808368331d0, &
826
       0.249954808368331d0, &
827
       0.249954808368331d0, &
828
       0.431647899262139d0, &
829
       0.249954808368331d0, &
830
       0.249954808368331d0, &
831
       0.249954808368331d0 ]
832
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: UWedge11P = [ &
833
       -0.062688380276010d0, &
834
       -0.062688380276010d0, &
835
       -0.874623239447980d0, &
836
       -0.333333333333333d0, &
837
       -0.798519188402179d0, &
838
       -0.798519188402179d0, &
839
       0.597038376804357d0, &
840
       -0.333333333333333d0, &
841
       -0.798519188402179d0, &
842
       -0.798519188402179d0, &
843
       0.597038376804357d0 ]
844
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: VWedge11P = [ &
845
       -0.062688380276010d0, &
846
       -0.874623239447980d0, &
847
       -0.062688380276010d0, &
848
       -0.333333333333333d0, &
849
       -0.798519188402179d0, &
850
       0.597038376804357d0, &
851
       -0.798519188402179d0, &
852
       -0.333333333333333d0, &
853
       -0.798519188402179d0, &
854
       0.597038376804357d0, &
855
       -0.798519188402179d0 ]
856
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: WWedge11P = [ &
857
       0.000000000000000d0, &
858
       0.000000000000000d0, &
859
       0.000000000000000d0, &
860
       0.866861974009030d0, &
861
       0.675639823682265d0, &
862
       0.675639823682265d0, &
863
       0.675639823682265d0, &
864
       -0.866861974009030d0, &
865
       -0.675639823682265d0, &
866
       -0.675639823682265d0, &
867
       -0.675639823682265d0 ]
868

869
!------------------------------------------------------------------------------
870
! Wedge - n=14, p=5
871
!------------------------------------------------------------------------------
872
   REAL(KIND=dp), DIMENSION(14), PRIVATE :: SWedge14P = [ &
873
       0.087576186678730d0, &
874
       0.229447629454892d0, &
875
       0.229447629454891d0, &
876
       0.833056798542985d0, &
877
       0.166443428304729d0, &
878
       0.166443428304729d0, &
879
       0.376993270712316d0, &
880
       0.170410864470884d0, &
881
       0.298194157223163d0, &
882
       0.298194157223162d0, &
883
       0.376993270712316d0, &
884
       0.170410864470884d0, &
885
       0.298194157223163d0, &
886
       0.298194157223162d0 ]
887
   REAL(KIND=dp), DIMENSION(14), PRIVATE :: UWedge14P = [ &
888
       -0.955901336867574d0, &
889
       -0.051621305926029d0, &
890
       -1.017063354038640d0, &
891
       -0.297388746460523d0, &
892
       0.774849157169622d0, &
893
       -0.849021640062097d0, &
894
       -0.665292008657551d0, &
895
       -0.012171148087346d0, &
896
       -0.734122164680096d0, &
897
       0.334122164680066d0, &
898
       -0.665292008657551d0, &
899
       -0.012171148087346d0, &
900
       -0.734122164680096d0, &
901
       0.334122164680066d0 ]
902
   REAL(KIND=dp), DIMENSION(14), PRIVATE :: VWedge14P = [ &
903
       -0.955901336867577d0, &
904
       -1.017063354038640d0, &
905
       -0.051621305926032d0, &
906
       -0.297388746460521d0, &
907
       -0.849021640062096d0, &
908
       0.774849157169620d0, &
909
       -0.665292008657551d0, &
910
       -0.012171148087349d0, &
911
       0.334122164680066d0, &
912
       -0.734122164680096d0, &
913
       -0.665292008657551d0, &
914
       -0.012171148087349d0, &
915
       0.334122164680066d0, &
916
       -0.734122164680096d0 ]
917
   REAL(KIND=dp), DIMENSION(14), PRIVATE :: WWedge14P = [ &
918
       0.000000000000286d0, &
919
       0.000000000000038d0, &
920
       0.000000000000038d0, &
921
       0.000000000000008d0, &
922
       0.000000000000080d0, &
923
       0.000000000000080d0, &
924
       0.756615409654429d0, &
925
       0.600149379161583d0, &
926
       0.808115770496521d0, &
927
       0.808115770496521d0, &
928
       -0.756615409654429d0, &
929
       -0.600149379161583d0, &
930
       -0.808115770496521d0, &
931
       -0.808115770496521d0 ]
932

933
!------------------------------------------------------------------------------
934
! Wedge - n=15, p=5
935
!------------------------------------------------------------------------------
936
   REAL(KIND=dp), DIMENSION(15), PRIVATE :: SWedge15P = [ &
937
       0.213895020288765d0, &
938
       0.141917375616806d0, &
939
       0.295568859378071d0, &
940
       0.256991945593379d0, &
941
       0.122121979248381d0, &
942
       0.175194917962627d0, &
943
       0.284969106392719d0, &
944
       0.323336131783334d0, &
945
       0.159056110329063d0, &
946
       0.748067388709644d0, &
947
       0.280551223607231d0, &
948
       0.147734016552639d0, &
949
       0.259874920383688d0, &
950
       0.235144061421191d0, &
951
       0.355576942732463d0 ]
952
   REAL(KIND=dp), DIMENSION(15), PRIVATE :: UWedge15P = [ &
953
       -0.820754415297359d0, &
954
       0.611831616907812d0, &
955
       -0.951379065092975d0, &
956
       0.200535109198601d0, &
957
       -0.909622841605196d0, &
958
       0.411514133287729d0, &
959
       -0.127534496411879d0, &
960
       -0.555217727817199d0, &
961
       0.706942532529193d0, &
962
       -0.278092963133809d0, &
963
       -0.057824844208300d0, &
964
       -0.043308436222116d0, &
965
       -0.774478920734726d0, &
966
       -0.765638443571458d0, &
967
       -0.732830649614460d0 ]
968
   REAL(KIND=dp), DIMENSION(15), PRIVATE :: VWedge15P = [ &
969
       0.701020947925133d0, &
970
       -0.869995576950389d0, &
971
       -0.087091980055873d0, &
972
       -0.783721735474016d0, &
973
       -0.890218158063352d0, &
974
       -0.725374126531787d0, &
975
       -0.953467283619037d0, &
976
       -0.530472194607007d0, &
977
       -0.782248553944540d0, &
978
       -0.291936530517119d0, &
979
       -0.056757587543798d0, &
980
       0.012890722780611d0, &
981
       0.476188541042454d0, &
982
       0.177195164202219d0, &
983
       -0.737447982744191d0 ]
984
   REAL(KIND=dp), DIMENSION(15), PRIVATE :: WWedge15P = [ &
985
       -0.300763696502910d0, &
986
       0.348546607420888d0, &
987
       0.150656042323906d0, &
988
       -0.844285153176719d0, &
989
       0.477120081549168d0, &
990
       0.864653509536562d0, &
991
       0.216019762875977d0, &
992
       0.873409672725819d0, &
993
       -0.390653804976705d0, &
994
       -0.126030507204870d0, &
995
       0.539907869785112d0, &
996
       -0.776314479909204d0, &
997
       0.745875967497062d0, &
998
       -0.888355356215127d0, &
999
       -0.651653242952189d0 ]
1000

1001
!------------------------------------------------------------------------------
1002
! Wedge - n=16, p=5
1003
!------------------------------------------------------------------------------
1004
   REAL(KIND=dp), DIMENSION(16), PRIVATE :: SWedge16P = [ &
1005
       0.711455555931488d0, &
1006
       0.224710067228267d0, &
1007
       0.224710067228267d0, &
1008
       0.224710067228267d0, &
1009
       0.185661421316158d0, &
1010
       0.185661421316158d0, &
1011
       0.185661421316158d0, &
1012
       0.250074285747794d0, &
1013
       0.250074285747794d0, &
1014
       0.250074285747794d0, &
1015
       0.185661421316158d0, &
1016
       0.185661421316158d0, &
1017
       0.185661421316158d0, &
1018
       0.250074285747794d0, &
1019
       0.250074285747794d0, &
1020
       0.250074285747794d0 ]
1021
   REAL(KIND=dp), DIMENSION(16), PRIVATE :: UWedge16P = [ &
1022
       -0.333333333333333d0, &
1023
       -0.025400070899509d0, &
1024
       -0.025400070899509d0, &
1025
       -0.949199858200983d0, &
1026
       -0.108803790659256d0, &
1027
       -0.108803790659256d0, &
1028
       -0.782392418681488d0, &
1029
       -0.798282108034583d0, &
1030
       -0.798282108034583d0, &
1031
       0.596564216069166d0, &
1032
       -0.108803790659256d0, &
1033
       -0.108803790659256d0, &
1034
       -0.782392418681488d0, &
1035
       -0.798282108034583d0, &
1036
       -0.798282108034583d0, &
1037
       0.596564216069166d0 ]
1038
   REAL(KIND=dp), DIMENSION(16), PRIVATE :: VWedge16P = [ &
1039
       -0.333333333333333d0, &
1040
       -0.025400070899509d0, &
1041
       -0.949199858200983d0, &
1042
       -0.025400070899509d0, &
1043
       -0.108803790659256d0, &
1044
       -0.782392418681488d0, &
1045
       -0.108803790659256d0, &
1046
       -0.798282108034583d0, &
1047
       0.596564216069166d0, &
1048
       -0.798282108034583d0, &
1049
       -0.108803790659256d0, &
1050
       -0.782392418681488d0, &
1051
       -0.108803790659256d0, &
1052
       -0.798282108034583d0, &
1053
       0.596564216069166d0, &
1054
       -0.798282108034583d0 ]
1055
   REAL(KIND=dp), DIMENSION(16), PRIVATE :: WWedge16P = [ &
1056
       0.000000000000000d0, &
1057
       0.000000000000000d0, &
1058
       0.000000000000000d0, &
1059
       0.000000000000000d0, &
1060
       0.871002934865444d0, &
1061
       0.871002934865444d0, &
1062
       0.871002934865444d0, &
1063
       0.570426980705159d0, &
1064
       0.570426980705159d0, &
1065
       0.570426980705159d0, &
1066
       -0.871002934865444d0, &
1067
       -0.871002934865444d0, &
1068
       -0.871002934865444d0, &
1069
       -0.570426980705159d0, &
1070
       -0.570426980705159d0, &
1071
       -0.570426980705159d0 ]
1072

1073
!------------------------------------------------------------------------------
1074
! Wedge - n=24, p=6
1075
!------------------------------------------------------------------------------
1076
   REAL(KIND=dp), DIMENSION(24), PRIVATE :: SWedge24P = [ &
1077
       0.168480079940386d0, &
1078
       0.168480079940386d0, &
1079
       0.168480079940386d0, &
1080
       0.168480079940386d0, &
1081
       0.168480079940386d0, &
1082
       0.168480079940386d0, &
1083
       0.000079282874851d0, &
1084
       0.000079282874851d0, &
1085
       0.000079282874851d0, &
1086
       0.544286440652304d0, &
1087
       0.026293733850586d0, &
1088
       0.283472344926041d0, &
1089
       0.283472344926041d0, &
1090
       0.283472344926041d0, &
1091
       0.115195615637235d0, &
1092
       0.115195615637235d0, &
1093
       0.115195615637235d0, &
1094
       0.026293733850586d0, &
1095
       0.283472344926041d0, &
1096
       0.283472344926041d0, &
1097
       0.283472344926041d0, &
1098
       0.115195615637235d0, &
1099
       0.115195615637235d0, &
1100
       0.115195615637235d0 ]
1101
   REAL(KIND=dp), DIMENSION(24), PRIVATE :: UWedge24P = [ &
1102
       0.513019949700545d0, &
1103
       -0.582925557762337d0, &
1104
       -0.582925557762337d0, &
1105
       0.513019949700545d0, &
1106
       -0.930094391938207d0, &
1107
       -0.930094391938207d0, &
1108
       -1.830988812620400d0, &
1109
       -1.830988812620400d0, &
1110
       2.661977625240810d0, &
1111
       -0.333333333333333d0, &
1112
       -0.333333333333333d0, &
1113
       -0.098283514203544d0, &
1114
       -0.098283514203544d0, &
1115
       -0.803432971592912d0, &
1116
       -0.812603471654584d0, &
1117
       -0.812603471654584d0, &
1118
       0.625206943309168d0, &
1119
       -0.333333333333333d0, &
1120
       -0.098283514203544d0, &
1121
       -0.098283514203544d0, &
1122
       -0.803432971592912d0, &
1123
       -0.812603471654584d0, &
1124
       -0.812603471654584d0, &
1125
       0.625206943309168d0 ]
1126
   REAL(KIND=dp), DIMENSION(24), PRIVATE :: VWedge24P = [ &
1127
       -0.930094391938207d0, &
1128
       -0.930094391938207d0, &
1129
       0.513019949700545d0, &
1130
       -0.582925557762337d0, &
1131
       -0.582925557762337d0, &
1132
       0.513019949700545d0, &
1133
       -1.830988812620400d0, &
1134
       2.661977625240810d0, &
1135
       -1.830988812620400d0, &
1136
       -0.333333333333333d0, &
1137
       -0.333333333333333d0, &
1138
       -0.098283514203544d0, &
1139
       -0.803432971592912d0, &
1140
       -0.098283514203544d0, &
1141
       -0.812603471654584d0, &
1142
       0.625206943309168d0, &
1143
       -0.812603471654584d0, &
1144
       -0.333333333333333d0, &
1145
       -0.098283514203544d0, &
1146
       -0.803432971592912d0, &
1147
       -0.098283514203544d0, &
1148
       -0.812603471654584d0, &
1149
       0.625206943309168d0, &
1150
       -0.812603471654584d0 ]
1151
   REAL(KIND=dp), DIMENSION(24), PRIVATE :: WWedge24P = [ &
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
       0.000000000000000d0, &
1162
       1.250521622121900d0, &
1163
       0.685008566774710d0, &
1164
       0.685008566774710d0, &
1165
       0.685008566774710d0, &
1166
       0.809574716992997d0, &
1167
       0.809574716992997d0, &
1168
       0.809574716992997d0, &
1169
       -1.250521622121900d0, &
1170
       -0.685008566774710d0, &
1171
       -0.685008566774710d0, &
1172
       -0.685008566774710d0, &
1173
       -0.809574716992997d0, &
1174
       -0.809574716992997d0, &
1175
       -0.809574716992997d0 ]
1176

1177
! continue wedge rules here
1178
!------------------------------------------------------------------------------
1179
! Wedge - n=27, p=6
1180
!------------------------------------------------------------------------------
1181
!   REAL(KIND=dp), DIMENSION(27), PRIVATE :: SWedge27P = [ &
1182
!   REAL(KIND=dp), DIMENSION(27), PRIVATE :: UWedge27P = [ &
1183
!   REAL(KIND=dp), DIMENSION(27), PRIVATE :: VWedge27P = [ &
1184
!   REAL(KIND=dp), DIMENSION(27), PRIVATE :: WWedge27P = [ &
1185

1186

1187

1188
 CONTAINS
1189

1190

1191
!------------------------------------------------------------------------------
1192
!> Subroutine to compute Fejer integration points by FFT.
1193
!------------------------------------------------------------------------------
1194
   SUBROUTINE ComputeFejerPoints1D( Points,Weights,n )
1195
!------------------------------------------------------------------------------
1196
    INTEGER :: n
1197
    REAL(KIND=dp) :: Points(n),Weights(n)
1198
!------------------------------------------------------------------------------
1199
    INTEGER :: i,j,k,l,m
1200
    TYPE(C_FFTCMPLX) :: W(n+1)
1201
    REAL(KIND=dp) :: arr((n+1)/2+1), V(n+2)
1202

1203
    DO i=1,(n+1)/2
1204
      Points(i) = COS(i*PI/(n+1._dp))
1205
      Points(n-i+1) = -Points(i)
1206
    END DO
1207

1208
    k = 0
1209
    DO i=1,n+1,2
1210
      k = k + 1
1211
      arr(k) = i
1212
    END DO
1213

1214
    V = 0._dp
1215
    DO i=1,k
1216
      V(i) = 2._dp/(arr(i)*(arr(i)-2))
1217
    END DO
1218
    V(k+1) = 1._dp/arr(k)
1219

1220
    DO i=1,n+1
1221
      W(i) % rval = -(V(i)+V(n-i+3))
1222
    END DO
1223
    CALL frfftb(n+1,W,V)
1224

1225
    Weights(1:n) = V(2:n+1)/(n+1)
1226
    DO i=1,n/2
1227
      Weights(i) = (Weights(i) + Weights(n-i+1))/2
1228
      Weights(n-i+1) = Weights(i)
1229
    END DO
1230
    Weights(1:n) = 2 * Weights(1:n) / SUM(Weights(1:n))
1231
!------------------------------------------------------------------------------
1232
  END SUBROUTINE ComputeFejerPoints1D
1233
!------------------------------------------------------------------------------
1234

1235

1236
!------------------------------------------------------------------------------
1237
!> Subroutine to compute Gaussian integration points and weights in [-1,1]
1238
!> as roots of Legendre polynomials.
1239
!------------------------------------------------------------------------------
1240
   SUBROUTINE ComputeGaussPoints1D( Points,Weights,n )
1241
!------------------------------------------------------------------------------
1242
     INTEGER :: n
1243
     REAL(KIND=dp) :: Points(n),Weights(n)
1244
!------------------------------------------------------------------------------
1245
     REAL(KIND=dp)   :: A(n/2,n/2),s,x,Work(8*n)
1246
     COMPLEX(KIND=dp) :: Eigs(n/2)
1247
     REAL(KIND=dp)   :: P(n+1),Q(n),P0(n),P1(n+1)
1248
     INTEGER :: i,j,k,np,info
1249
!------------------------------------------------------------------------------
1250
! One point is trivial
1251
!------------------------------------------------------------------------------
1252
     IF ( n <= 1 ) THEN
1253
       Points(1)  = 0.0d0
1254
       Weights(1) = 2.0d0
1255
       RETURN
1256
     END IF
1257
!------------------------------------------------------------------------------
1258
! Compute coefficients of n:th Legendre polynomial from the recurrence:
1259
!
1260
! (i+1)P_{i+1}(x) = (2i+1)*x*P_i(x) - i*P_{i-1}(x), P_{0} = 1; P_{1} = x;
1261
!
1262
! CAVEAT: Computed coefficients inaccurate for n > ~15
1263
!------------------------------------------------------------------------------
1264
     P = 0
1265
     P0(1) = 1
1266
     P1(1) = 1
1267
     P1(2) = 0
1268

1269
     DO i=1,n-1
1270
       P(1:i+1) = (2*i+1) * P1(1:i+1)  / (i+1)
1271
       P(3:i+2) = P(3:i+2) - i*P0(1:i) / (i+1)
1272
       P0(1:i+1) = P1(1:i+1); P1(1:i+2) = P(1:i+2)
1273
     END DO
1274
!------------------------------------------------------------------------------
1275
! Odd n implicates zero as one of the roots...
1276
!------------------------------------------------------------------------------
1277
     np = n - MOD(n,2)
1278
!------------------------------------------------------------------------------
1279
!  Variable substitution: y=x^2
1280
!------------------------------------------------------------------------------
1281
     np = np / 2
1282
     DO i=1,np+1
1283
       P(i) = P(2*i-1)
1284
     END DO
1285
!------------------------------------------------------------------------------
1286
! Solve the roots of the polynomial by forming a matrix whose characteristic
1287
! polynomial is the n:th Legendre polynomial and solving for the eigenvalues.
1288
! Dunno if this is a very good method....
1289
!------------------------------------------------------------------------------
1290
     A=0
1291
     DO i=1,np-1
1292
       A(i,i+1) = 1
1293
     END DO
1294

1295
     DO i=1,np
1296
       A(np,i) = -P(np+2-i) / P(1)
1297
     END DO
1298
     CALL DGEEV( 'N','N',np,A,n/2,Points,P0,Work,1,Work,1,Work,8*n,info )
1299
!    CALL EigenValues( A, np, Eigs )
1300
!    Points(1:np) = REAL( Eigs(1:np) )
1301

1302
!------------------------------------------------------------------------------
1303
! Backsubstitute from y=x^2
1304
!------------------------------------------------------------------------------
1305
     Q(1:np+1) = P(1:np+1)
1306
     P = 0
1307
     DO i=1,np+1
1308
       P(2*i-1) = Q(i)
1309
     END DO
1310

1311
     Q(1:np) = Points(1:np)
1312
     DO i=1,np
1313
       Points(2*i-1) = +SQRT( Q(i) )
1314
       Points(2*i)   = -SQRT( Q(i) )
1315
     END DO
1316
     IF ( MOD(n,2) == 1 ) Points(n) = 0.0d0
1317

1318
     CALL DerivPoly( n,Q,P )
1319
     CALL RefineRoots( n,P,Q,Points )
1320
!------------------------------------------------------------------------------
1321
! Check for roots
1322
!------------------------------------------------------------------------------
1323
     DO i=1,n
1324
       s = EvalPoly( n,P,Points(i) )
1325
       IF ( ABS(s) > 1.0d-12 ) THEN
1326
         CALL Warn( 'ComputeGaussPoints1D', &
1327
                 '-------------------------------------------------------' )
1328
         CALL Warn( 'ComputeGaussPoints1D', 'Computed integration point' )
1329
         CALL Warn( 'ComputeGaussPoints1D', 'seems to be inaccurate: ' )
1330
         WRITE( Message, * ) 'Points req.: ',n
1331
         CALL Warn( 'ComputeGaussPoints1D', Message )
1332
         WRITE( Message, * ) 'Residual: ',s
1333
         CALL Warn( 'ComputeGaussPoints1D', Message )
1334
         WRITE( Message, * ) 'Point: +-', SQRT(Points(i))
1335
         CALL Warn( 'ComputeGaussPoints1D', Message )
1336
         CALL Warn( 'ComputeGaussPoints1D', &
1337
                 '-------------------------------------------------------' )
1338
       END IF
1339
     END DO
1340
!------------------------------------------------------------------------------
1341
! Finally, the integration weights equal to
1342
!
1343
! W_i = 2/( (1-x_i^2)*Q(x_i)^2 ), x_i is the i:th root, and Q(x) = dP(x) / dx
1344
!------------------------------------------------------------------------------
1345
     CALL DerivPoly( n,Q,P )
1346

1347
     DO i=1,n
1348
       s = EvalPoly( n-1,Q,Points(i) )
1349
       Weights(i) = 2 / ((1-Points(i)**2)*s**2);
1350
     END DO
1351
!------------------------------------------------------------------------------
1352
! ... make really sure the weights add up:
1353
!------------------------------------------------------------------------------
1354
     Weights(1:n) = 2 * Weights(1:n) / SUM(Weights(1:n))
1355

1356
CONTAINS
1357

1358
!--------------------------------------------------------------------------
1359

1360
   FUNCTION EvalPoly( n,P,x ) RESULT(s)
1361
     INTEGER :: i,n
1362
     REAL(KIND=dp) :: P(n+1),x,s
1363

1364
     s = 0.0d0
1365
     DO i=1,n+1
1366
       s = s * x + P(i)
1367
     END DO
1368
   END FUNCTION EvalPoly
1369

1370
!--------------------------------------------------------------------------
1371

1372
   SUBROUTINE DerivPoly( n,Q,P )
1373
     INTEGER :: i,n
1374
     REAL(KIND=dp) :: Q(n),P(n+1)
1375

1376
     DO i=1,n
1377
       Q(i) = P(i)*(n-i+1)
1378
     END DO
1379
   END SUBROUTINE DerivPoly
1380

1381
!--------------------------------------------------------------------------
1382

1383
   SUBROUTINE RefineRoots( n,P,Q,Points )
1384
     INTEGER :: i,j,n
1385
     REAL(KIND=dp) :: P(n+1),Q(n),Points(n)
1386

1387
     REAL(KIND=dp) :: x,s
1388
     INTEGER, PARAMETER :: MaxIter = 100
1389

1390
     DO i=1,n
1391
       x = Points(i)
1392
       DO j=1,MaxIter
1393
         s = EvalPoly(n,P,x) / EvalPoly(n-1,Q,x)
1394
         x = x - s
1395
         IF ( ABS(s) <= ABS(x)*EPSILON(s) ) EXIT
1396
       END DO
1397
       IF ( ABS(EvalPoly(n,P,x))<ABS(EvalPoly(n,P,Points(i))) ) THEN
1398
         IF ( ABS(x-Points(i))<1.0d-8 ) Points(i) = x
1399
       END IF
1400
     END DO
1401
   END SUBROUTINE RefineRoots
1402

1403
!--------------------------------------------------------------------------
1404
 END SUBROUTINE ComputeGaussPoints1D
1405
!--------------------------------------------------------------------------
1406

1407
!--------------------------------------------------------------------------
1408
   FUNCTION GaussPointsInitialized() RESULT(g)
1409
!--------------------------------------------------------------------------
1410
     LOGICAL :: g
1411
!--------------------------------------------------------------------------
1412
     g=Ginit
1413
!--------------------------------------------------------------------------
1414
   END FUNCTION GaussPointsInitialized
1415
!--------------------------------------------------------------------------
1416

1417
! !------------------------------------------------------------------------------
1418
!    SUBROUTINE GaussPointsInit
1419
! !------------------------------------------------------------------------------
1420
!      INTEGER :: i,n,istat
1421
!      INTEGER :: nthreads, thread, omp_get_max_threads, omp_get_thread_num
1422
!
1423
!      IF ( .NOT. ALLOCATED(IntegStuff) ) THEN
1424
! !$omp barrier
1425
! !$omp critical
1426
!        IF ( .NOT. GInit ) THEN
1427
!          GInit = .TRUE.
1428
!          nthreads = 1
1429
! !$       nthreads = omp_get_max_threads()
1430
!          ALLOCATE( IntegStuff(nthreads) )
1431
!          DO i=1,nthreads
1432
!            IntegStuff(i) % u => NULL()
1433
!            IntegStuff(i) % v => NULL()
1434
!            IntegStuff(i) % w => NULL()
1435
!            IntegStuff(i) % s => NULL()
1436
!          END DO
1437
!          DO n=1,MAXN
1438
!            CALL ComputeGaussPoints1D( Points(1:n,n),Weights(1:n,n),n )
1439
!          END DO
1440
!        END IF
1441
! !$omp end critical
1442
!      END IF
1443
!
1444
!      thread = 1
1445
! !$   thread = omp_get_thread_num()+1
1446
!     ALLOCATE( IntegStuff(thread) % u(MAX_INTEGRATION_POINTS), &
1447
!               IntegStuff(thread) % v(MAX_INTEGRATION_POINTS), &
1448
!               IntegStuff(thread) % w(MAX_INTEGRATION_POINTS), &
1449
!               IntegStuff(thread) % s(MAX_INTEGRATION_POINTS), STAT=istat )
1450
!
1451
!     IF ( istat /= 0 ) THEN
1452
!       CALL Fatal( 'GaussPointsInit', 'Memory allocation error.' )
1453
!       STOP
1454
!     END IF
1455
! !------------------------------------------------------------------------------
1456
! !  END SUBROUTINE GaussPointsInit
1457
! !------------------------------------------------------------------------------
1458

1459
!------------------------------------------------------------------------------
1460
   SUBROUTINE GaussPointsInit
1461
!------------------------------------------------------------------------------
1462
     INTEGER :: i,n,istat
1463

1464
     IF ( .NOT. GInit ) THEN
1465
        DO n=1,MAXN
1466
          CALL ComputeGaussPoints1D( Points(1:n,n),Weights(1:n,n),n )
1467
        END DO
1468
        GInit = .TRUE.
1469
     END IF
1470

1471
     ALLOCATE( IntegStuff % u(MAX_INTEGRATION_POINTS), &
1472
               IntegStuff % v(MAX_INTEGRATION_POINTS), &
1473
               IntegStuff % w(MAX_INTEGRATION_POINTS), &
1474
               IntegStuff % s(MAX_INTEGRATION_POINTS), STAT=istat )
1475
     IntegStuff % u = 0._dp
1476
     IntegStuff % v = 0._dp
1477
     IntegStuff % w = 0._dp
1478
     IntegStuff % s = 0._dp
1479
     IF ( istat /= 0 ) THEN
1480
       CALL Fatal( 'GaussPointsInit', 'Memory allocation error.' )
1481
     END IF
1482
!------------------------------------------------------------------------------
1483
  END SUBROUTINE GaussPointsInit
1484
!------------------------------------------------------------------------------
1485

1486

1487
!------------------------------------------------------------------------------
1488
   FUNCTION GaussPoints0D( n ) RESULT(p)
1489
!------------------------------------------------------------------------------
1490
      INTEGER :: n
1491
      TYPE(GaussIntegrationPoints_t), POINTER :: p
1492
!     INTEGER :: thread, omp_get_thread_num
1493

1494
      IF ( .NOT. GInit ) CALL GaussPointsInit
1495
!     thread = 1
1496
! !$    thread = omp_get_thread_num()+1
1497
!     p => IntegStuff(thread)
1498
      p => IntegStuff
1499
      p % n = 1
1500
      p % u(1) = 0
1501
      p % v(1) = 0
1502
      p % w(1) = 0
1503
      p % s(1) = 1
1504
!------------------------------------------------------------------------------
1505
   END FUNCTION GaussPoints0D
1506
!------------------------------------------------------------------------------
1507

1508

1509
!------------------------------------------------------------------------------
1510
!>    Return Gaussian integration points for 1D line element
1511
!------------------------------------------------------------------------------
1512
   FUNCTION GaussPoints1D( n ) RESULT(p)
1513
!------------------------------------------------------------------------------
1514
      INTEGER :: n   !< number of points in the requested rule
1515
      TYPE(GaussIntegrationPoints_t), POINTER :: p
1516
!------------------------------------------------------------------------------
1517
!     INTEGER :: thread, omp_get_thread_num
1518

1519
      IF ( .NOT. GInit ) CALL GaussPointsInit
1520
!     thread = 1
1521
! !$    thread = omp_get_thread_num()+1
1522
!      p => IntegStuff(thread)
1523
      p => IntegStuff
1524
      IF ( n < 1 .OR. n > MAXN ) THEN
1525
        p % n = 0
1526
        WRITE( Message, * ) 'Invalid number of points: ',n
1527
        CALL Fatal( 'GaussPoints1D', Message )
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,'(A,I0,A,I0)') 'Invalid number of points (max=',MAXN,'): ', n
1764
        CALL Fatal( 'GaussPointsQuad', Message )
1765
      END IF
1766

1767
      t = 0
1768
      DO i=1,n
1769
!DIR$ IVDEP
1770
        DO j=1,n
1771
          t = t + 1
1772
          p % u(t) = Points(j,n)
1773
          p % v(t) = Points(i,n)
1774
          p % s(t) = Weights(i,n)*Weights(j,n)
1775
        END DO
1776
      END DO
1777
      p % n = t
1778
!------------------------------------------------------------------------------
1779
   END FUNCTION GaussPointsQuad
1780
!------------------------------------------------------------------------------
1781

1782

1783
!------------------------------------------------------------------------------
1784
   FUNCTION GaussPointsPTetra(np) RESULT(p)
1785
!------------------------------------------------------------------------------
1786
   INTEGER :: i,np,n
1787
   TYPE(GaussIntegrationPoints_t), POINTER :: p
1788
   REAL(KIND=dp) :: uh, vh, wh, sh
1789
!  INTEGER :: thread, omp_get_thread_num
1790
!------------------------------------------------------------------------------
1791
   IF ( .NOT. GInit ) CALL GaussPointsInit
1792
!    thread = 1
1793
! !$ thread = omp_get_thread_num()+1
1794
!    p => IntegStuff(thread)
1795
   p => IntegStuff
1796
   n = DBLE(np)**(1.0D0/3.0D0) + 0.5D0
1797

1798
   ! Get Gauss points of p brick
1799
   ! (take into account term z^2) from jacobian determinant
1800
   p = GaussPointsPBrick(n,n,n+1)
1801
   ! p = GaussPointsBrick( np )
1802
   ! WRITE (*,*) 'Getting Gauss points for: ', n, p % n
1803

1804
   ! For each point apply mapping from brick to
1805
   ! tetrahedron and multiply each weight by detJ
1806
   ! of mapping
1807
!DIR$ IVDEP
1808
   DO i=1,p % n
1809
      uh = p % u(i)
1810
      vh = p % v(i)
1811
      wh = p % w(i)
1812
      sh = p % s(i)
1813

1814
      p % u(i)= 1d0/4*(uh - uh*vh - uh*wh + uh*vh*wh)
1815
      p % v(i)= SQRT(3d0)/4*(5d0/3 + vh - wh/3 - vh*wh)
1816
      p % w(i)= SQRT(6d0)/3*(1d0 + wh)
1817
      p % s(i)= -sh * SQRT(2d0)/16 * (1d0 - vh - wh + vh*wh) * (-1d0 + wh)
1818
   END DO
1819
!------------------------------------------------------------------------------
1820
 END FUNCTION GaussPointsPTetra
1821
!------------------------------------------------------------------------------
1822

1823
!------------------------------------------------------------------------------
1824
!>    Return Gaussian integration points for 3D tetrahedral element. Here the
1825
!>    reference element over which the integration is done can also be the
1826
!>    regular tetrahedron used in the description of p-elements. In that case,
1827
!>    this routine may return a more economical set of integration points.
1828
!------------------------------------------------------------------------------
1829
   FUNCTION GaussPointsTetra( n, PReferenceElement ) RESULT(p)
1830
!------------------------------------------------------------------------------
1831
      INTEGER :: n      !< number of points in the requested rule
1832
      LOGICAL, OPTIONAL ::  PReferenceElement !< used for switching the reference element
1833
      TYPE(GaussIntegrationPoints_t), POINTER :: p
1834
!------------------------------------------------------------------------------
1835
      REAL( KIND=dp ) :: ScaleFactor
1836
      INTEGER :: i
1837
      LOGICAL :: ConvertToPTetrahedron
1838
      REAL (KIND=dp) :: uq, vq, wq, sq
1839
!     INTEGER :: thread, omp_get_thread_num
1840
!----------------------------------------------------------------------------------
1841
      ConvertToPTetrahedron = .FALSE.
1842
      IF ( PRESENT(PReferenceElement) ) THEN
1843
         ConvertToPTetrahedron =  PReferenceElement
1844
      END IF
1845

1846
      IF ( .NOT. GInit ) CALL GaussPointsInit
1847
!       thread = 1
1848
! !$    thread = omp_get_thread_num()+1
1849
!       p => IntegStuff(thread)
1850
      p => IntegStuff
1851

1852
      SELECT CASE (n)
1853
      CASE (1)
1854
         p % u(1:n) = UTetra1P
1855
         p % v(1:n) = VTetra1P
1856
         p % w(1:n) = WTetra1P
1857
         p % s(1:n) = STetra1P / 6.0D0
1858
         p % n = 1
1859
      CASE (4)
1860
         p % u(1:n) = UTetra4P
1861
         p % v(1:n) = VTetra4P
1862
         p % w(1:n) = WTetra4P
1863
         p % s(1:n) = STetra4P / 6.0D0
1864
         p % n = 4
1865
      CASE (5)
1866
         p % u(1:n) = UTetra5P
1867
         p % v(1:n) = VTetra5P
1868
         p % w(1:n) = WTetra5P
1869
         p % s(1:n) = STetra5P / 6.0D0
1870
         p % n = 5
1871
      CASE (11)
1872
         p % u(1:n) = UTetra11P
1873
         p % v(1:n) = VTetra11P
1874
         p % w(1:n) = WTetra11P
1875
         p % s(1:n) = STetra11P / 6.0D0
1876
         p % n = 11
1877
      CASE DEFAULT
1878
!        CALL Error( 'GaussPointsTetra', 'Invalid number of points requested.' )
1879
!        p % n = 0
1880
         p = GaussPointsBrick( n )
1881

1882
!DIR$ IVDEP
1883
         DO i=1,p % n
1884
            ScaleFactor = 0.5d0
1885
            p % u(i) = ( p % u(i) + 1 ) * Scalefactor
1886
            p % v(i) = ( p % v(i) + 1 ) * ScaleFactor
1887
            p % w(i) = ( p % w(i) + 1 ) * ScaleFactor
1888
            p % s(i) = p % s(i) * ScaleFactor**3
1889

1890
            ScaleFactor = 1.0d0 - p % w(i)
1891
            p % u(i) = p % u(i) * ScaleFactor
1892
            p % v(i) = p % v(i) * ScaleFactor
1893
            p % s(i) = p % s(i) * ScaleFactor**2
1894

1895
            ScaleFactor = 1.0d0 - p % v(i) / ScaleFactor
1896
            p % u(i) = p % u(i) * ScaleFactor
1897
            p % s(i) = p % s(i) * ScaleFactor
1898
         END DO
1899
!         p % s(1:p % n) = p % s(1:p % n) / SUM( p % s(1:p % n) ) / 6.0d0
1900
      END SELECT
1901

1902
      IF (ConvertToPTetrahedron) THEN
1903
         !-------------------------------------------------------------------
1904
         ! Apply an additional transformation if the actual reference element
1905
         ! is the regular tetrahedron used in the description of p-elements
1906
	 ! We map the original integration points into their counterparts on the
1907
	 ! p-reference element and scale the weights by the determinant of the
1908
	 ! deformation gradient associated with the change of reference element.
1909
         !-------------------------------------------------------------------
1910
!DIR$ IVDEP
1911
        DO i=1,P % n
1912
            uq = P % u(i)
1913
            vq = P % v(i)
1914
            wq = P % w(i)
1915
            sq = P % s(i)
1916
            P % u(i) = -1.0d0 + 2.0d0*uq + vq + wq
1917
            P % v(i) = SQRT(3.0d0)*vq + 1.0d0/SQRT(3.0d0)*wq
1918
            P % w(i) = SQRT(8.0d0)/SQRT(3.0d0)*wq
1919
            P % s(i) = SQRT(8.0d0)*2.0d0*sq
1920
         END DO
1921
      END IF
1922
!------------------------------------------------------------------------------
1923
   END FUNCTION GaussPointsTetra
1924
!------------------------------------------------------------------------------
1925

1926
!------------------------------------------------------------------------------
1927
   FUNCTION GaussPointsPPyramid( np ) RESULT(p)
1928
!------------------------------------------------------------------------------
1929
   INTEGER :: np,n,i
1930
   REAL(KIND=dp) :: uh,vh,wh,sh
1931
   TYPE(GaussIntegrationPoints_t), POINTER :: p
1932
!  INTEGER :: thread, omp_get_thread_num
1933
!------------------------------------------------------------------------------
1934
   IF ( .NOT. GInit ) CALL GaussPointsInit
1935
!    thread = 1
1936
! !$ thread = omp_get_thread_num()+1
1937
!    p => IntegStuff(thread)
1938
   p => IntegStuff
1939

1940
   n = DBLE(np)**(1.0D0/3.0D0) + 0.5D0
1941

1942
   ! Get Gauss points of p brick
1943
   ! (take into account term (-1+z)^2) from jacobian determinant
1944
   p = GaussPointsPBrick(n,n,n+1)
1945

1946
   ! For each point apply mapping from brick to
1947
   ! pyramid and multiply each weight by detJ
1948
   ! of mapping
1949
!DIR$ IVDEP
1950
   DO i=1,p % n
1951
      uh = p % u(i)
1952
      vh = p % v(i)
1953
      wh = p % w(i)
1954
      sh = p % s(i)
1955

1956
      p % u(i)= 1d0/2*uh*(1d0-wh)
1957
      p % v(i)= 1d0/2*vh*(1d0-wh)
1958
      p % w(i)= SQRT(2d0)/2*(1d0+wh)
1959
      p % s(i)= sh * SQRT(2d0)/8 * (-1d0+wh)**2
1960
   END DO
1961
!------------------------------------------------------------------------------
1962
   END FUNCTION GaussPointsPPyramid
1963
!------------------------------------------------------------------------------
1964

1965

1966
!------------------------------------------------------------------------------
1967
!>    Return Gaussian integration points for 3D pyramid element.
1968
!------------------------------------------------------------------------------
1969
   FUNCTION GaussPointsPyramid( np ) RESULT(p)
1970
!------------------------------------------------------------------------------
1971
      INTEGER :: np     !< number of points in the requested rule
1972
      TYPE(GaussIntegrationPoints_t), POINTER :: p
1973
!------------------------------------------------------------------------------
1974
      INTEGER :: i,j,k,n,t
1975
!       INTEGER :: thread, omp_get_thread_num
1976

1977
      IF ( .NOT. GInit ) CALL GaussPointsInit
1978
!       thread = 1
1979
! !$    thread = omp_get_thread_num()+1
1980
!       p => IntegStuff(thread)
1981
      p => IntegStuff
1982

1983
      n = REAL(np)**(1.0D0/3.0D0) + 0.5D0
1984

1985
      IF ( n < 1 .OR. n > MAXN ) THEN
1986
         p % n = 0
1987
         WRITE( Message, * ) 'Invalid number of points: ', n
1988
         CALL Fatal( 'GaussPointsPyramid', Message )
1989
      END IF
1990

1991
      t = 0
1992
      DO i=1,n
1993
        DO j=1,n
1994
!DIR$ IVDEP
1995
          DO k=1,n
1996
             t = t + 1
1997
             p % u(t) = Points(k,n)
1998
             p % v(t) = Points(j,n)
1999
             p % w(t) = Points(i,n)
2000
             p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n)
2001
          END DO
2002
        END DO
2003
      END DO
2004
      p % n = t
2005

2006
!DIR$ IVDEP
2007
      DO t=1,p % n
2008
        p % w(t) = (p % w(t) + 1.0d0) / 2.0d0
2009
        p % u(t) = p % u(t) * (1.0d0-p % w(t))
2010
        p % v(t) = p % v(t) * (1.0d0-p % w(t))
2011
        p % s(t) = p % s(t) * (1.0d0-p % w(t))**2/2
2012
      END DO
2013
!------------------------------------------------------------------------------
2014
   END FUNCTION GaussPointsPyramid
2015
!------------------------------------------------------------------------------
2016

2017
!------------------------------------------------------------------------------
2018
   FUNCTION GaussPointsPWedge(n) RESULT(p)
2019
!------------------------------------------------------------------------------
2020
   INTEGER :: n, i
2021
   REAL(KIND=dp) :: uh,vh,wh,sh
2022
   TYPE(GaussIntegrationPoints_t), POINTER :: p
2023
!   INTEGER :: thread, omp_get_thread_num
2024
!------------------------------------------------------------------------------
2025
   IF ( .NOT. GInit ) CALL GaussPointsInit
2026
!    thread = 1
2027
! !$ thread = omp_get_thread_num()+1
2028
!    p => IntegStuff(thread)
2029
   p => IntegStuff
2030

2031
   ! Get Gauss points of brick
2032
   p = GaussPointsBrick(n)
2033

2034
   ! For each point apply mapping from brick to
2035
   ! wedge and multiply each weight by detJ
2036
   ! of mapping
2037
!DIR$ IVDEP
2038
   DO i=1,p % n
2039
      uh = p % u(i)
2040
      vh = p % v(i)
2041
      wh = p % w(i)
2042
      sh = p % s(i)
2043

2044
      p % u(i)= 1d0/2*(uh-uh*vh)
2045
      p % v(i)= SQRT(3d0)/2*(1d0+vh)
2046
      p % w(i)= wh
2047
      p % s(i)= sh * SQRT(3d0)*(1-vh)/4
2048
   END DO
2049
!------------------------------------------------------------------------------
2050
   END FUNCTION GaussPointsPWedge
2051
!------------------------------------------------------------------------------
2052

2053

2054
   
2055
!------------------------------------------------------------------------------
2056
!>    Return Gaussian integration points for 3D wedge element
2057
!------------------------------------------------------------------------------
2058
   FUNCTION GaussPointsWedge( np ) RESULT(p)
2059
!------------------------------------------------------------------------------
2060
      INTEGER :: np     !< number of points in the requested rule
2061
      TYPE(GaussIntegrationPoints_t), POINTER :: p
2062
!------------------------------------------------------------------------------
2063
      INTEGER :: i,j,k,n,t
2064
!       INTEGER :: thread, omp_get_thread_num
2065
!------------------------------------------------------------------------------
2066
      IF ( .NOT. GInit ) CALL GaussPointsInit
2067
!       thread = 1
2068
! !$    thread = omp_get_thread_num()+1
2069
!       p => IntegStuff(thread)
2070
      p => IntegStuff
2071

2072
      n = REAL(np)**(1.0d0/3.0d0) + 0.5d0
2073

2074
      IF ( n < 1 .OR. n > MAXN ) THEN
2075
         p % n = 0
2076
         WRITE( Message, * ) 'Invalid number of points: ', n
2077
         CALL Fatal( 'GaussPointsWedge', Message )
2078
      END IF
2079

2080
      t = 0
2081
      DO i=1,n
2082
        DO j=1,n
2083
!DIR$ IVDEP
2084
          DO k=1,n
2085
             t = t + 1
2086
             p % u(t) = Points(k,n)
2087
             p % v(t) = Points(j,n)
2088
             p % w(t) = Points(i,n)
2089
             p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n)
2090
          END DO
2091
        END DO
2092
      END DO
2093
      p % n = t
2094

2095
!DIR$ IVDEP
2096
      DO i=1,p % n
2097
        p % v(i) = (p % v(i) + 1)/2
2098
        p % u(i) = (p % u(i) + 1)/2 * (1 - p % v(i))
2099
        p % s(i) = p % s(i) * (1-p % v(i))/4
2100
      END DO
2101
!------------------------------------------------------------------------------
2102
   END FUNCTION GaussPointsWedge
2103
!------------------------------------------------------------------------------
2104

2105
!------------------------------------------------------------------------------
2106
!>  Return an optimized number of Gaussian points for integrating over prisms.
2107
!>  Here the reference element can also be that of the p-approximation.
2108
!>  A rule with m x n points is returned.
2109
!------------------------------------------------------------------------------
2110
   FUNCTION GaussPointsWedge2(m,n,PReferenceElement) RESULT(p)
2111
!------------------------------------------------------------------------------
2112
      TYPE(GaussIntegrationPoints_t), POINTER :: p
2113
      INTEGER :: m     !< The number of points over a triangular face
2114
      INTEGER :: n     !< The number of points in the orthogonal direction to the triangular faces
2115
      LOGICAL, OPTIONAL ::  PReferenceElement !< Used for switching the reference element
2116
!------------------------------------------------------------------------------
2117
      INTEGER :: i,j,k,t
2118
      LOGICAL :: ConvertToPPrism
2119
      REAL(KIND=dp) :: uq(20), vq(20), sq(20)
2120
      REAL(KIND=dp) :: uh,vh,wh,sh
2121
!-----------------------------------------------------------------------------
2122
      ConvertToPPrism = .FALSE.
2123
      IF ( PRESENT(PReferenceElement) ) THEN
2124
         ConvertToPPrism =  PReferenceElement
2125
      END IF
2126
      IF ( .NOT. GInit ) CALL GaussPointsInit
2127
      p => IntegStuff
2128

2129
      SELECT CASE (m)
2130
      CASE (1)
2131
         IF (ConvertToPPrism) THEN
2132
            uq(1) = -1.0d0 + 2.0d0*UTriangle1P(1) + VTriangle1P(1)
2133
            vq(1) = SQRT(3.0d0) * VTriangle1P(1)
2134
            sq(1) = SQRT(3.0d0) * STriangle1P(1)
2135
         ELSE
2136
            uq(1:m) = UTriangle1P
2137
            vq(1:m) = VTriangle1P
2138
            sq(1:m) = STriangle1P / 2.0D0
2139
         END IF
2140
      CASE (3)
2141
        IF (ConvertToPPrism) THEN
2142
!DIR$ IVDEP
2143
            DO i=1,m
2144
               uq(i) = -1.0d0 + 2.0d0*UTriangle3P(i) + VTriangle3P(i)
2145
               vq(i) = SQRT(3.0d0) * VTriangle3P(i)
2146
               sq(i) = SQRT(3.0d0) * STriangle3P(i)
2147
            END DO
2148
         ELSE
2149
            uq(1:m) = UTriangle3P
2150
            vq(1:m) = VTriangle3P
2151
            sq(1:m) = STriangle3P / 2.0D0
2152
         END IF
2153
      CASE (4)
2154
        IF (ConvertToPPrism) THEN
2155
!DIR$ IVDEP
2156
            DO i=1,m
2157
               uq(i) = -1.0d0 + 2.0d0*UTriangle4P(i) + VTriangle4P(i)
2158
               vq(i) = SQRT(3.0d0) * VTriangle4P(i)
2159
               sq(i) = SQRT(3.0d0) * STriangle4P(i)
2160
            END DO
2161
         ELSE
2162
            uq(1:m) = UTriangle4P
2163
            vq(1:m) = VTriangle4P
2164
            sq(1:m) = STriangle4P / 2.0D0
2165
         END IF
2166
      CASE (6)
2167
        IF (ConvertToPPrism) THEN
2168
!DIR$ IVDEP
2169
            DO i=1,m
2170
               uq(i) = -1.0d0 + 2.0d0*UTriangle6P(i) + VTriangle6P(i)
2171
               vq(i) = SQRT(3.0d0) * VTriangle6P(i)
2172
               sq(i) = SQRT(3.0d0) * STriangle6P(i)
2173
            END DO
2174
         ELSE
2175
            uq(1:m) = UTriangle6P
2176
            vq(1:m) = VTriangle6P
2177
            sq(1:m) = STriangle6P / 2.0D0
2178
         END IF
2179
      CASE (7)
2180
        IF (ConvertToPPrism) THEN
2181
!DIR$ IVDEP
2182
            DO i=1,m
2183
               uq(i) = -1.0d0 + 2.0d0*UTriangle7P(i) + VTriangle7P(i)
2184
               vq(i) = SQRT(3.0d0) * VTriangle7P(i)
2185
               sq(i) = SQRT(3.0d0) * STriangle7P(i)
2186
            END DO
2187
         ELSE
2188
            uq(1:m) = UTriangle7P
2189
            vq(1:m) = VTriangle7P
2190
            sq(1:m) = STriangle7P / 2.0D0
2191
         END IF
2192
      CASE (11)
2193
        IF (ConvertToPPrism) THEN
2194
!DIR$ IVDEP
2195
            DO i=1,m
2196
               uq(i) = -1.0d0 + 2.0d0*UTriangle11P(i) + VTriangle11P(i)
2197
               vq(i) = SQRT(3.0d0) * VTriangle11P(i)
2198
               sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle11P(i)
2199
            END DO
2200
         ELSE
2201
            uq(1:m) = UTriangle11P
2202
            vq(1:m) = VTriangle11P
2203
            sq(1:m) = STriangle11P
2204
         END IF
2205
      CASE (12)
2206
        IF (ConvertToPPrism) THEN
2207
!DIR$ IVDEP
2208
            DO i=1,m
2209
               uq(i) = -1.0d0 + 2.0d0*UTriangle12P(i) + VTriangle12P(i)
2210
               vq(i) = SQRT(3.0d0) * VTriangle12P(i)
2211
               sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle12P(i)
2212
            END DO
2213
         ELSE
2214
            uq(1:m) = UTriangle12P
2215
            vq(1:m) = VTriangle12P
2216
            sq(1:m) = STriangle12P
2217
         END IF
2218
      CASE (17)
2219
        IF (ConvertToPPrism) THEN
2220
!DIR$ IVDEP
2221
            DO i=1,m
2222
               uq(i) = -1.0d0 + 2.0d0*UTriangle17P(i) + VTriangle17P(i)
2223
               vq(i) = SQRT(3.0d0) * VTriangle17P(i)
2224
               sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle17P(i)
2225
            END DO
2226
         ELSE
2227
            uq(1:m) = UTriangle17P
2228
            vq(1:m) = VTriangle17P
2229
            sq(1:m) = STriangle17P
2230
         END IF
2231
      CASE (20)
2232
        IF (ConvertToPPrism) THEN
2233
!DIR$ IVDEP
2234
            DO i=1,m
2235
               uq(i) = -1.0d0 + 2.0d0*UTriangle20P(i) + VTriangle20P(i)
2236
               vq(i) = SQRT(3.0d0) * VTriangle20P(i)
2237
               sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle20P(i)
2238
            END DO
2239
         ELSE
2240
            uq(1:m) = UTriangle20P
2241
            vq(1:m) = VTriangle20P
2242
            sq(1:m) = STriangle20P
2243
         END IF
2244
      CASE DEFAULT
2245
         !-------------------------------------------------------------------------------------
2246
         ! The requested number of integration points associated with the triangular face
2247
         ! cannot be created. Therefore use the given number of points in the third direction
2248
         ! and generate silently n x n x n rule.
2249
         !-----------------------------------------------------------------------------------
2250
         IF (ConvertToPPrism) THEN
2251
           p = GaussPointsBrick(n*n*n)
2252
!DIR$ IVDEP
2253
            DO i=1,p % n
2254
               uh = p % u(i)
2255
               vh = p % v(i)
2256
               wh = p % w(i)
2257
               sh = p % s(i)
2258

2259
               p % u(i)= 1d0/2*(uh-uh*vh)
2260
               p % v(i)= SQRT(3d0)/2*(1d0+vh)
2261
               p % w(i)= wh
2262
               p % s(i)= sh * SQRT(3d0)*(1-vh)/4
2263
            END DO
2264
         ELSE
2265
            t = 0
2266
            DO i=1,n
2267
              DO j=1,n
2268
!DIR$ IVDEP
2269
                  DO k=1,n
2270
                     t = t + 1
2271
                     p % u(t) = Points(k,n)
2272
                     p % v(t) = Points(j,n)
2273
                     p % w(t) = Points(i,n)
2274
                     p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n)
2275
                  END DO
2276
               END DO
2277
            END DO
2278
            p % n = t
2279

2280
!DIR$ IVDEP
2281
            DO i=1,p % n
2282
               p % v(i) = (p % v(i) + 1)/2
2283
               p % u(i) = (p % u(i) + 1)/2 * (1 - p % v(i))
2284
               p % s(i) = p % s(i) * (1-p % v(i))/4
2285
            END DO
2286
         END IF
2287
         RETURN
2288
      END SELECT
2289

2290
      t = 0
2291
      DO i=1,m
2292
!DIR$ IVDEP
2293
         DO j=1,n
2294
            t = t + 1
2295
            p % u(t) = uq(i)
2296
            p % v(t) = vq(i)
2297
            p % w(t) = Points(j,n)
2298
            p % s(t) = sq(i)*Weights(j,n)
2299
         END DO
2300
      END DO
2301
      p % n = t
2302
!------------------------------------------------------------------------------
2303
    END FUNCTION GaussPointsWedge2
2304
!------------------------------------------------------------------------------
2305

2306

2307
!------------------------------------------------------------------------------
2308
!>    Return Gaussian integration points for 3D wedge elements using
2309
!> economical quadratures that are not product of segment and triangle rules.
2310
!------------------------------------------------------------------------------
2311
   FUNCTION GaussPointsWedgeEconomic( n, PReferenceElement ) RESULT(p)
2312
!------------------------------------------------------------------------------
2313
      INTEGER :: n      !< number of points in the requested rule
2314
      LOGICAL, OPTIONAL ::  PReferenceElement !< used for switching the reference element 
2315
      TYPE(GaussIntegrationPoints_t), POINTER :: p
2316
!------------------------------------------------------------------------------
2317
      REAL( KIND=dp ) :: ScaleFactor
2318
      INTEGER :: i
2319
      LOGICAL :: ConvertToPWedge
2320
      REAL (KIND=dp) :: uq, vq, wq, sq
2321
!     INTEGER :: thread, omp_get_thread_num
2322
!----------------------------------------------------------------------------------
2323
      ConvertToPWedge = .FALSE.
2324
      IF ( PRESENT(PReferenceElement) ) THEN
2325
        ConvertToPWedge = PReferenceElement
2326
      END IF
2327

2328
      IF ( .NOT. GInit ) CALL GaussPointsInit
2329
!       thread = 1
2330
! !$    thread = omp_get_thread_num()+1
2331
!       p => IntegStuff(thread)
2332
      p => IntegStuff
2333

2334
      SELECT CASE (n)
2335
      CASE (4)
2336
         p % u(1:n) = UWedge4P
2337
         p % v(1:n) = VWedge4P
2338
         p % w(1:n) = WWedge4P
2339
         p % s(1:n) = SWedge4P
2340
      CASE (5)
2341
         p % u(1:n) = UWedge5P
2342
         p % v(1:n) = VWedge5P
2343
         p % w(1:n) = WWedge5P
2344
         p % s(1:n) = SWedge5P 
2345
      CASE (6)
2346
         p % u(1:n) = UWedge6P
2347
         p % v(1:n) = VWedge6P
2348
         p % w(1:n) = WWedge6P
2349
         p % s(1:n) = SWedge6P 
2350
      CASE (7)
2351
         p % u(1:n) = UWedge7P
2352
         p % v(1:n) = VWedge7P
2353
         p % w(1:n) = WWedge7P
2354
         p % s(1:n) = SWedge7P
2355
      CASE (10)
2356
         p % u(1:n) = UWedge10P
2357
         p % v(1:n) = VWedge10P
2358
         p % w(1:n) = WWedge10P
2359
         p % s(1:n) = SWedge10P 
2360
      CASE (11)
2361
         p % u(1:n) = UWedge11P
2362
         p % v(1:n) = VWedge11P
2363
         p % w(1:n) = WWedge11P
2364
         p % s(1:n) = SWedge11P 
2365
      CASE (14)
2366
         p % u(1:n) = UWedge14P
2367
         p % v(1:n) = VWedge14P
2368
         p % w(1:n) = WWedge14P
2369
         p % s(1:n) = SWedge14P 
2370
      CASE (15)
2371
         p % u(1:n) = UWedge15P
2372
         p % v(1:n) = VWedge15P
2373
         p % w(1:n) = WWedge15P
2374
         p % s(1:n) = SWedge15P 
2375
      CASE (16)
2376
         p % u(1:n) = UWedge16P
2377
         p % v(1:n) = VWedge16P
2378
         p % w(1:n) = WWedge16P
2379
         p % s(1:n) = SWedge16P
2380
      CASE (24)
2381
         p % u(1:n) = UWedge24P
2382
         p % v(1:n) = VWedge24P
2383
         p % w(1:n) = WWedge24P
2384
         p % s(1:n) = SWedge24P
2385
         
2386
      CASE DEFAULT
2387
        CALL Fatal( 'GaussPointsWedgeEconomic',&
2388
            'Invalid number of points requested')
2389
      END SELECT
2390

2391
      p % n = n
2392

2393
      IF (ConvertToPWedge ) THEN
2394
        DO i=1,P % n  
2395
          uq = P % u(i) 
2396
          vq = P % v(i)
2397
          sq = P % s(i)
2398

2399
          ! first to classical
2400
          uq = (uq+1.d0)/2.0d0
2401
          vq = (vq+1.d0)/2.0d0
2402

2403
          ! then to p-convention
2404
          P % u(i) = -1.0d0 + 2.0d0*uq + vq
2405
          P % v(i) = SQRT(3.0d0) * vq            
2406
          P % s(i) = SQRT(3.0d0) * sq / 2.0d0
2407
        END DO
2408
      ELSE
2409
        ! Map to classical Elmer local coordinates in [0,1]
2410
        p % u(1:n) = ( p % u(1:n)+1.0d0 ) / 2.0d0
2411
        p % v(1:n) = ( p % v(1:n)+1.0d0 ) / 2.0d0 
2412
        p % s(1:n) = p % s(1:n) / 4.0d0
2413
      END IF
2414
!------------------------------------------------------------------------------
2415
    END FUNCTION GaussPointsWedgeEconomic
2416
!------------------------------------------------------------------------------
2417

2418
    
2419

2420
!------------------------------------------------------------------------------
2421
!>    Return Gaussian integration points for 3D brick element for
2422
!>    composite rule
2423
!>    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}))).
2424
!------------------------------------------------------------------------------
2425
   FUNCTION GaussPointsPBrick( nx, ny, nz ) RESULT(p)
2426
!------------------------------------------------------------------------------
2427
      INTEGER :: nx    !< number of points in the requested rule in x direction
2428
      INTEGER :: ny    !< number of points in the requested rule in y direction
2429
      INTEGER :: nz    !< number of points in the requested rule in z direction
2430
      TYPE(GaussIntegrationPoints_t), POINTER :: p
2431
!------------------------------------------------------------------------------
2432
      INTEGER i,j,k,t
2433
!       INTEGER :: thread, omp_get_thread_num
2434
!------------------------------------------------------------------------------
2435
      IF ( .NOT. GInit ) CALL GaussPointsInit
2436
!       thread = 1
2437
! !$    thread = omp_get_thread_num()+1
2438
!       p => IntegStuff(thread)
2439
      p => IntegStuff
2440

2441
      ! Check validity of number of integration points
2442
      IF ( nx < 1 .OR. nx > MAXN .OR. &
2443
           ny < 1 .OR. ny > MAXN .OR. &
2444
           nz < 1 .OR. nz > MAXN) THEN
2445
        p % n = 0
2446
        WRITE( Message, * ) 'Invalid number of points: ', nx, ny, nz
2447
        CALL Fatal( 'GaussPointsBrick', Message )
2448
      END IF
2449

2450
      t = 0
2451
      DO i=1,nx
2452
        DO j=1,ny
2453
!DIR$ IVDEP
2454
          DO k=1,nz
2455
            t = t + 1
2456
            p % u(t) = Points(i,nx)
2457
            p % v(t) = Points(j,ny)
2458
            p % w(t) = Points(k,nz)
2459
            p % s(t) = Weights(i,nx)*Weights(j,ny)*Weights(k,nz)
2460
          END DO
2461
        END DO
2462
      END DO
2463
      p % n = t
2464
!------------------------------------------------------------------------------
2465
    END FUNCTION GaussPointsPBrick
2466
!------------------------------------------------------------------------------
2467

2468

2469
!------------------------------------------------------------------------------
2470
!>    Return Gaussian integration points for 3D brick element
2471
!------------------------------------------------------------------------------
2472
   FUNCTION GaussPointsBrick( np ) RESULT(p)
2473
!------------------------------------------------------------------------------
2474
      INTEGER :: np     !< number of points in the requested rule
2475
      TYPE(GaussIntegrationPoints_t), POINTER :: p
2476
!------------------------------------------------------------------------------
2477
      INTEGER i,j,k,n,t
2478
!      INTEGER :: thread, omp_get_thread_num
2479

2480
      IF ( .NOT. GInit ) CALL GaussPointsInit
2481
!      thread = 1
2482
! !$    thread = omp_get_thread_num()+1
2483
!       p => IntegStuff(thread)
2484
      p => IntegStuff
2485

2486
      SELECT CASE( np )
2487
      CASE( 8 )
2488
        n = 2
2489
      CASE( 27 )
2490
        n = 3
2491
      CASE( 64 )
2492
        n = 4
2493
      CASE DEFAULT
2494
        n = REAL(np)**(1.0D0/3.0D0) + 0.5D0
2495
        
2496
        IF ( n < 1 .OR. n > MAXN ) THEN
2497
          p % n = 0
2498
          WRITE( Message, * ) 'Invalid number of points: ', n
2499
          CALL Fatal( 'GaussPointsBrick', Message )
2500
        END IF
2501
      END SELECT
2502
        
2503
      t = 0
2504
      DO i=1,n
2505
        DO j=1,n
2506
!DIR$ IVDEP
2507
          DO k=1,n
2508
            t = t + 1
2509
            p % u(t) = Points(k,n)
2510
            p % v(t) = Points(j,n)
2511
            p % w(t) = Points(i,n)
2512
            p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n)
2513
          END DO
2514
        END DO
2515
      END DO
2516
      p % n = t
2517
!------------------------------------------------------------------------------
2518
   END FUNCTION GaussPointsBrick
2519
!------------------------------------------------------------------------------
2520

2521

2522
!------------------------------------------------------------------------------
2523
!>    Given element structure return Gauss integration points for the element.
2524
!----------------------------------------------------------------------------------------------
2525
   FUNCTION GaussPoints( elm, np, RelOrder, EdgeBasis, PReferenceElement, &
2526
        EdgeBasisDegree) RESULT(IntegStuff)
2527
!---------------------------------------------------------------------------------------------
2528
     USE PElementMaps, ONLY : isActivePElement
2529
     TYPE( Element_t ) :: elm   !< Structure holding element description
2530
     INTEGER, OPTIONAL :: np    !< Number of requested points
2531
     INTEGER, OPTIONAL :: RelOrder  !< Relative order of Gauss points (-1,0,+1)
2532
     LOGICAL, OPTIONAL :: EdgeBasis !< For choosing quadrature suitable for edge elements
2533
     LOGICAL, OPTIONAL :: PReferenceElement !< For switching to the p-version reference element
2534
     INTEGER, OPTIONAL :: EdgeBasisDegree ! The degree of edge elements
2535
     TYPE( GaussIntegrationPoints_t ) :: IntegStuff   !< Structure holding the integration points
2536
!------------------------------------------------------------------------------
2537
     LOGICAL :: pElement, UsePRefElement, Economic
2538
     INTEGER :: n, eldim, p1d, ntri, nseg, necon
2539
     TYPE(ElementType_t), POINTER :: elmt
2540
!------------------------------------------------------------------------------
2541
     elmt => elm % TYPE
2542

2543
     IF (PRESENT(EdgeBasis)) THEN
2544
       IF (EdgeBasis) THEN
2545
         UsePRefElement = .TRUE.
2546
         IF (PRESENT(PReferenceElement)) UsePRefElement = PReferenceElement
2547
         IF (PRESENT(EdgeBasisDegree)) THEN
2548
           IntegStuff = EdgeElementGaussPoints(elmt % ElementCode/100, UsePRefElement, EdgeBasisDegree)
2549
         ELSE
2550
           IntegStuff = EdgeElementGaussPoints(elmt % ElementCode/100, UsePRefElement)
2551
         END IF
2552
         RETURN
2553
       END IF
2554
     END IF
2555

2556
     IF( PRESENT(PReferenceElement)) THEN
2557
       pElement = PReferenceElement
2558
     ELSE
2559
       pElement = isActivePElement(elm)
2560
     END IF
2561
     
2562
     IF ( PRESENT(np) ) THEN
2563
       n = np
2564
     ELSE IF( PRESENT( RelOrder ) ) THEN
2565
       IF (pElement) THEN
2566
         n = elm % PDefs % GaussPoints
2567
         IF( RelOrder == 0 ) THEN
2568
           CONTINUE
2569
         ELSE
2570
           eldim = elm % type % dimension
2571
           p1d = NINT( n**(1.0_dp/eldim) ) + RelOrder
2572
           IF( p1d < 1 ) THEN
2573
             CALL Fatal('GaussPoints','Number of integration points must remain positive!')
2574
           END IF
2575
           n = p1d**eldim
2576

2577
           ! Take into account the case of economic integration. Map a tensor
2578
           ! product rule to a corresponding economic rule.
2579
           IF (elm % PDefs % Serendipity .AND. elmt % ElementCode / 100 == 4) THEN
2580
             SELECT CASE(n)
2581
             CASE(9)
2582
               n = 8
2583
             CASE(16)
2584
               n = 12
2585
             CASE(25)
2586
               n = 20
2587
             CASE(36)
2588
               n = 25
2589
             CASE(49)
2590
               n = 36
2591
             CASE(64)
2592
               n = 45
2593
             CASE(81)
2594
               n = 60
2595
             CASE DEFAULT
2596
               CONTINUE
2597
             END SELECT
2598
           END IF
2599
         END IF
2600
       ELSE
2601
         IF( RelOrder == 0 ) THEN
2602
           n = elmt % GaussPoints
2603
         ELSE IF( RelOrder == 1 ) THEN
2604
           n = elmt % GaussPoints2
2605
         ELSE IF( RelOrder == -1 ) THEN
2606
           n = elmt % GaussPoints0
2607
         ELSE
2608
           PRINT *,'RelOrder can only be {-1, 0, 1} !'
2609
         END IF
2610
       END IF
2611
     ELSE
2612
       IF (pElement) THEN
2613
         IF( ASSOCIATED( elm % PDefs ) ) THEN
2614
           n = elm % PDefs % GaussPoints
2615
         ELSE
2616
           n = elmt % GaussPoints
2617
         END IF
2618
         IF( n == 0 ) n = elmt % GaussPoints
2619
       ELSE
2620
         n = elmt % GaussPoints
2621
       END IF
2622
     END IF
2623

2624
     SELECT CASE( elmt % ElementCode / 100 )
2625
     CASE (1)
2626
        IntegStuff = GaussPoints0D(n)
2627

2628
     CASE (2)
2629
        IntegStuff = GaussPoints1D(n)
2630

2631
     CASE (3)
2632
        IF (pElement) THEN
2633
          IntegStuff = GaussPointsPTriangle(n)
2634
        ELSE
2635
          IntegStuff = GaussPointsTriangle(n)
2636
        END IF
2637

2638
     CASE (4)
2639
       IF (pElement .AND. ASSOCIATED( elm % pdefs ) ) THEN
2640
         Economic = .FALSE.
2641
         IF(elm % pdefs % Serendipity) THEN
2642
           ! For certain polynomial orders, economic quadratures may be used:
2643
           IF (elm % PDefs % P > 1 .AND.  elm % PDefs % P <= 7) Economic = .TRUE.
2644
           ! An explicit bubble augmentation with lower-order methods switches to
2645
           ! the standard rule:
2646
           IF (elm % BDOFs > 0 .AND. elm % PDefs % P < 4) Economic = .FALSE.
2647
           ! The economic 12-point rule appears to be somehow faulty, so we 
2648
           ! shall never call it
2649
         END IF
2650

2651
         IF (Economic .AND. n==12) THEN
2652
           IntegStuff = GaussPointsQuad(16)
2653
         ELSE
2654
           IntegStuff = GaussPointsQuad(n, Economic)
2655
         END IF
2656
       ELSE
2657
         IntegStuff = GaussPointsQuad(n)
2658
       END IF
2659

2660
     CASE (5)
2661
        IF (pElement) THEN
2662
           IntegStuff = GaussPointsPTetra(n)
2663
        ELSE
2664
           IntegStuff = GaussPointsTetra(n)
2665
        END IF
2666

2667
     CASE (6)
2668
        IF (pElement) THEN
2669
           IntegStuff = GaussPointsPPyramid(n)
2670
        ELSE
2671
           IntegStuff = GaussPointsPyramid(n)
2672
        END IF
2673

2674
      CASE (7)
2675
        IF( PRESENT( np ) ) THEN
2676
          ! possible values:
2677
          ! triangle = 1, 3, 4, 6, 7, 11, 12, 17, 20
2678
          ! segment  = 1, 2, 3, 4, 5, 6,  7,  8,  9,
2679

2680
          ntri = 0; nseg = 0; necon = 0
2681

2682
          SELECT CASE( n )
2683

2684
            ! Cases ordered so that we have the segment x triangle rules first
2685
          CASE( 1, 2, 3)
2686
            nseg = n
2687
          CASE( 6, 8 )
2688
            nseg = 2
2689
          CASE( 12, 18, 21 )
2690
            nseg = 3
2691
          CASE( 28, 44, 48 )
2692
            nseg = 4
2693
          CASE( 85, 100 )
2694
            nseg = 5
2695

2696
            ! The economical rules
2697
          CASE( 4, 5, 7, 10, 11, 14, 15, 16, 24 )
2698
            ! Note: we would have 6 and 8 point rules from the economic family as well
2699
            necon = n 
2700
          END SELECT
2701

2702
          IF( nseg > 0 ) THEN
2703
            ntri =  n / nseg
2704
            IntegStuff = GaussPointsWedge2(ntri,nseg,PReferenceElement=pElement)
2705
            RETURN
2706
          ELSE IF( necon > 0 ) THEN
2707
            IntegStuff = GaussPointsWedgeEconomic(necon,PReferenceElement=pElement)
2708
            RETURN
2709
          END IF
2710
        END IF
2711
        
2712
        IF (pElement) THEN
2713
           IntegStuff = GaussPointsPWedge(n)
2714
        ELSE
2715
           IntegStuff = GaussPointsWedge(n)
2716
        END IF
2717

2718
     CASE (8)       
2719
       IntegStuff = GaussPointsBrick(n)
2720
       
2721
     END SELECT
2722
       
2723
   END FUNCTION GaussPoints
2724

2725

2726

2727
!------------------------------------------------------------------------------
2728
!>  Given element structure return Gauss integration points at the element corners.
2729
!----------------------------------------------------------------------------------------------
2730
   FUNCTION CornerGaussPoints( elm, EdgeBasis, PReferenceElement ) RESULT(CornerStuff)
2731
!---------------------------------------------------------------------------------------------
2732
     USE PElementMaps, ONLY : isActivePElement
2733
     TYPE( Element_t ) :: elm
2734
     LOGICAL, OPTIONAL :: EdgeBasis
2735
     LOGICAL, OPTIONAL :: PReferenceElement 
2736
     TYPE( GaussIntegrationPoints_t ) :: CornerStuff   
2737
!------------------------------------------------------------------------------
2738
     LOGICAL :: pElement
2739
     INTEGER :: n, i, j, k, t, ecode
2740
     TYPE( GaussIntegrationPoints_t), POINTER :: ip
2741
!------------------------------------------------------------------------------
2742
     ecode = elm % TYPE % ElementCode
2743

2744
     IF (PRESENT(PReferenceElement)) THEN
2745
       pElement = PReferenceElement
2746
     ELSE IF (PRESENT(EdgeBasis)) THEN
2747
       pElement = .TRUE.
2748
     ELSE
2749
       pElement = isActivePElement(elm)
2750
     END IF
2751

2752
     IF(.NOT. Ginit) CALL GaussPointsInit()
2753
     ip => IntegStuff
2754
     
2755
     ! Compute the number of corner nodes
2756
     n = ecode / 100
2757
     IF( n >= 5 .AND. n <= 7 ) n = n-1 
2758
     ip % n = n
2759
     ip % s(1:n) = 1.0_dp / n
2760
     
2761
     SELECT CASE( ecode  / 100 )
2762

2763
     CASE( 3 )
2764
       ip % u(1) = 0.0_dp; ip % v(1) = 0.0_dp
2765
       ip % u(2) = 1.0_dp; ip % v(2) = 0.0_dp
2766
       ip % u(3) = 0.0_dp; ip % v(3) = 1.0_dp
2767
       
2768
     CASE( 4 )
2769
       t = 0
2770
       DO i=1,2
2771
         DO j=1,2
2772
           t = t+1
2773
           ip % u(t) = (-1.0_dp)**i
2774
           ip % v(t) = (-1.0_dp)**j
2775
         END DO
2776
       END DO
2777

2778
     CASE( 5 )
2779
       ip % u(1) = 0.0_dp; ip % v(1) = 0.0_dp; ip % w(1) = 0.0_dp
2780
       ip % u(2) = 1.0_dp; ip % v(2) = 0.0_dp; ip % w(2) = 0.0_dp
2781
       ip % u(3) = 0.0_dp; ip % v(3) = 1.0_dp; ip % w(3) = 0.0_dp
2782
       ip % u(4) = 0.0_dp; ip % v(4) = 0.0_dp; ip % w(4) = 1.0_dp
2783
       
2784
     CASE( 8 )
2785
       t = 0
2786
       DO i=1,2
2787
         DO j=1,2
2788
           DO k=1,2
2789
             t = t+1
2790
             ip % u(t) = (-1.0_dp)**i
2791
             ip % v(t) = (-1.0_dp)**j
2792
             ip % w(t) = (-1.0_dp)**k
2793
           END DO
2794
         END DO
2795
       END DO
2796

2797
     CASE DEFAULT
2798
       WRITE(Message,'(A,I0)') 'Not implemented for element type: ', ecode
2799
       CALL Fatal('CornerGaussPoints',Message)
2800
     END SELECT
2801

2802
     IF( pElement ) THEN
2803
       DO i=1,n
2804
         CALL ConvertToPReference(ecode,ip % u(i),ip % v(i),ip % w(i))            
2805
       END DO
2806
     END IF
2807
     
2808
#if 0
2809
     PRINT *,'Corner Gauss Points:',n
2810
     PRINT *,'u:',IP % u(1:n)
2811
     PRINT *,'v:',IP % v(1:n)
2812
     IF(ecode > 500) PRINT *,'w:',IP % w(1:n)
2813
#endif
2814
     
2815
     CornerStuff = ip
2816
     
2817
   END FUNCTION CornerGaussPoints
2818

2819

2820
!------------------------------------------------------------------------------
2821
!>  Given element structure return integration points at the element center only.
2822
!----------------------------------------------------------------------------------------------
2823
   FUNCTION CenterGaussPoints( elm, EdgeBasis, PReferenceElement ) RESULT(CornerStuff)
2824
!---------------------------------------------------------------------------------------------
2825
     USE PElementMaps, ONLY : isActivePElement
2826
     TYPE( Element_t ) :: elm
2827
     LOGICAL, OPTIONAL :: EdgeBasis
2828
     LOGICAL, OPTIONAL :: PReferenceElement 
2829
     TYPE( GaussIntegrationPoints_t ) :: CornerStuff   
2830
!------------------------------------------------------------------------------
2831
     LOGICAL :: pElement
2832
     INTEGER :: n, i, j, k, t, ecode
2833
     TYPE( GaussIntegrationPoints_t), POINTER :: ip
2834
!------------------------------------------------------------------------------
2835
     ecode = elm % TYPE % ElementCode
2836

2837
     IF (PRESENT(PReferenceElement)) THEN
2838
       pElement = PReferenceElement
2839
     ELSE IF (PRESENT(EdgeBasis)) THEN
2840
       pElement = .TRUE.
2841
     ELSE
2842
       pElement = isActivePElement(elm)
2843
     END IF
2844

2845
     IF(.NOT. Ginit) CALL GaussPointsInit()
2846
     ip => IntegStuff
2847
     
2848
     ! Compute the number of corner nodes
2849
     n = ecode / 100
2850
     IF( n >= 5 .AND. n <= 7 ) n = n-1 
2851
     ip % n = 1
2852
     ip % s(1:n) = 1.0_dp
2853
     
2854
     SELECT CASE( ecode  / 100 )
2855

2856
     CASE( 3 )
2857
       ip % u(1) = 1.0_dp/3.0_dp; ip % v(1) = 1.0_dp/3.0_dp
2858
       
2859
     CASE( 4 )
2860
       ip % u(1) = 0.0_dp
2861
       ip % v(1) = 0.0_dp
2862

2863
     CASE( 5 )
2864
       ip % u(1) = 1.0_dp/4.0_dp; ip % v(1) = 1.0_dp/4.0_dp; ip % w(1) = 1.0_dp/4.0_dp
2865
       
2866
     CASE( 8 )
2867
       ip % u(1) = 0.0_dp
2868
       ip % v(1) = 0.0_dp
2869
       ip % w(1) = 0.0_dp
2870

2871
     CASE DEFAULT
2872
       WRITE(Message,'(A,I0)') 'Not implemented for element type: ', ecode
2873
       CALL Fatal('CornerGaussPoints',Message)
2874
     END SELECT
2875

2876
     IF( pElement ) THEN
2877
       DO i=1,n
2878
         CALL ConvertToPReference(ecode,ip % u(i),ip % v(i),ip % w(i))            
2879
       END DO
2880
     END IF
2881
     
2882
#if 0
2883
     PRINT *,'Corner Gauss Points:',n
2884
     PRINT *,'u:',IP % u(1:n)
2885
     PRINT *,'v:',IP % v(1:n)
2886
     IF(ecode > 500) PRINT *,'w:',IP % w(1:n)
2887
#endif
2888
     
2889
     CornerStuff = ip
2890
     
2891
   END FUNCTION CenterGaussPoints
2892

2893
   
2894
   
2895
!----------------------------------------------------------------------------------
2896
!>  Return a suitable version of the Gaussian numerical integration method for
2897
!>  H(curl)-conforming finite elements. The default here is that the edge basis
2898
!>  functions are obtained via calling the function EdgeElementInfo, so that
2899
!>  the reference element is chosen to be that used for p-approximation.
2900
!>  By giving the optional argument PiolaVersion = .FALSE. this function returns
2901
!>  an alternate (usually smaller) set of integration points on the classic
2902
!>  reference element. This option may be useful when the edge basis functions are
2903
!>  obtained via calling the alternate subroutine GetEdgeBasis.
2904
!----------------------------------------------------------------------------------------
2905
   FUNCTION EdgeElementGaussPoints(ElementFamily, PiolaVersion, BasisDegree) RESULT(IP)
2906
!---------------------------------------------------------------------------------------
2907
     INTEGER :: ElementFamily
2908
     LOGICAL, OPTIONAL :: PiolaVersion
2909
     INTEGER, OPTIONAL :: BasisDegree
2910
     TYPE(GaussIntegrationPoints_t) :: IP
2911
!------------------------------------------------------------------------------
2912
     LOGICAL :: PRefElement, SecondOrder
2913
!------------------------------------------------------------------------------
2914
     PRefElement = .TRUE.
2915
     SecondOrder = .FALSE.
2916
     IF ( PRESENT(PiolaVersion) ) PRefElement = PiolaVersion
2917
     IF ( PRESENT(BasisDegree) ) SecondOrder = BasisDegree > 1
2918

2919
     SELECT CASE(ElementFamily)
2920
     CASE (1)
2921
        IP = GaussPoints0D(1)
2922
     CASE (2)
2923
        IP = GaussPoints1D(2)
2924
     CASE(3)
2925
        IF (SecondOrder) THEN
2926
           IP = GaussPointsTriangle(6, PReferenceElement=PRefElement)
2927
        ELSE
2928
           IP = GaussPointsTriangle(3, PReferenceElement=PRefElement)
2929
        END IF
2930
     CASE(4)
2931
        IF (PRefElement) THEN
2932
           IP = GaussPointsQuad(9)
2933
        ELSE
2934
           IP = GaussPointsQuad(4)
2935
        END IF
2936
     CASE(5)
2937
        IF (SecondOrder) THEN
2938
           IP = GaussPointsTetra(11, PReferenceElement=PRefElement)
2939
        ELSE
2940
           IP = GaussPointsTetra(4, PReferenceElement=PRefElement)
2941
        END IF
2942
     CASE(6)
2943
        IF (PRefElement) THEN
2944
           IP = GaussPointsPPyramid(27)
2945
        ELSE
2946
           IP = GaussPointsPyramid(27)
2947
        END IF
2948
     CASE(7)
2949
        IF (PRefElement) THEN
2950
           IP = GaussPointsWedge2(6,3,PReferenceElement=PRefElement)
2951
        ELSE
2952
           IP = GaussPointsWedge2(3,2,PReferenceElement=PRefElement)
2953
        END IF
2954
     CASE(8)
2955
        IF (PRefElement) THEN
2956
           IP = GaussPointsPBrick(3,3,3)
2957
        ELSE
2958
           IP = GaussPointsBrick(8)
2959
        END IF
2960
     CASE DEFAULT
2961
        CALL Fatal('Integration::EdgeElementGaussPoints','Unsupported element type')
2962
     END SELECT
2963
!------------------------------------------------------------------------------
2964
   END FUNCTION EdgeElementGaussPoints
2965
!------------------------------------------------------------------------------
2966

2967
   SUBROUTINE ConvertToPReference(ElementCode,u,v,w,s)
2968
     INTEGER :: ElementCode
2969
     REAL(KIND=dp) :: u,v,w
2970
     REAL(KIND=dp), OPTIONAL :: s
2971
          
2972
     
2973
     SELECT CASE( ElementCode / 100 )
2974

2975
     CASE(1,2,4,8)
2976
       ! Nothing to do, reference element is the same
2977
       
2978
     CASE(3)
2979
       ! triangle
2980
       u = 2*u + v - 1
2981
       v = SQRT(3.0_dp)*v
2982
       IF( PRESENT(s) ) s = SQRT(3.0_dp)*2*s
2983
        
2984
     CASE(5)
2985
       ! tetrahedron
2986
       u = 2*u + v + w - 1
2987
       v = SQRT(3._dp)*v + 1/SQRT(3._dp)*w
2988
       w = 2*SQRT(2/3._dp)*w
2989
       IF(PRESENT(s)) s = 4*SQRT(2.0d0)*s       
2990

2991
     CASE(6) 
2992
       ! pyramid
2993
       w = SQRT(2._dp)*w
2994
       ! scaling of s??
2995
       
2996
     CASE(7)
2997
       ! wedge / prism
2998
       u = 2*u + v - 1
2999
       v = SQRT(3._dp)*v
3000
       IF(PRESENT(s)) s = 2*SQRT(3.0d0) * s       
3001

3002
     CASE DEFAULT
3003
       CALL Fatal('Integration::ConvertToPReference','Unsupported element type')
3004
     END SELECT
3005
            
3006
     
3007
   END SUBROUTINE ConvertToPReference
3008
   
3009

3010
!---------------------------------------------------------------------------
3011
END MODULE Integration
3012
!---------------------------------------------------------------------------
3013

3014
!> \} ElmerLib
3015

3016
Product

Resources

Company
