CoCalc -- H1Basis.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fem/src/H1Basis.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: Mikko Byckling, Juha Ruokolainen
27
! *  Web:     http://www.csc.fi/elmer
28
! *  Address: CSC - IT Center for Science Ltd.
29
! *           Keilaranta 14
30
! *           02101 Espoo, Finland 
31
! *
32
! *  Original Date: 31 May 2017
33
! *
34
! *****************************************************************************/
35

36
!> \ingroup ElmerLib
37
!> \{
38

39
!-----------------------------------------------------------------------------
40
!>  Module defining vectorized p element basis functions and mappings.
41
!-----------------------------------------------------------------------------
42

43
#include "../config.h"
44

45
MODULE H1Basis
46
  USE Messages
47
  USE Types, ONLY : dp, VECTOR_BLOCK_LENGTH
48
  
49
  ! Module contains vectorized version of FE basis
50
  ! functions for selected elements
51

52
  INTEGER, PARAMETER :: H1Basis_MaxPElementEdgeNodes=2, &
53
          H1Basis_MaxPElementEdges = 12, &
54
          H1Basis_MaxPElementFaceNodes = 4, &
55
          H1Basis_MaxPElementFaces = 6
56
CONTAINS
57

58
  SUBROUTINE H1Basis_GetEdgeDirection(ecode, nedges, globalind, direction)
59
    IMPLICIT NONE
60
    INTEGER, INTENT(IN) :: ecode, nedges
61
    INTEGER, DIMENSION(:), POINTER CONTIG, INTENT(IN) :: globalind
62
    INTEGER, DIMENSION(H1Basis_MaxPElementEdgeNodes, &
63
                       H1Basis_MaxPElementEdges), INTENT(INOUT) :: direction
64
    INTEGER :: i, tmp
65

66
    CALL H1Basis_GetEdgeMap(ecode, direction)
67

68
!DIR$ LOOP COUNT MAX=12
69
    DO i=1,nedges
70
      ! Swap local indices if needed to set the global direction (increasing)
71
      IF (globalind(direction(2,i))<globalind(direction(1,i))) THEN
72
        tmp = direction(1,i)
73
        direction(1,i)=direction(2,i)
74
        direction(2,i)=tmp
75
      END IF
76
    END DO
77
  END SUBROUTINE H1Basis_GetEdgeDirection
78
  
79
  SUBROUTINE H1Basis_GetTetraEdgeDirection(ttype, direction)
80
    IMPLICIT NONE
81
    INTEGER, INTENT(IN) :: ttype
82
    INTEGER, DIMENSION(H1Basis_MaxPElementEdgeNodes, &
83
                       H1Basis_MaxPElementEdges), TARGET, INTENT(INOUT) :: direction
84
    ! Local variables
85
    INTEGER, DIMENSION(:), POINTER CONTIG :: flatdir
86

87
    ! Remap array bounds
88
    flatdir(1:H1Basis_MaxPElementEdgeNodes*H1Basis_MaxPElementEdges) => direction
89
    SELECT CASE (ttype)
90
    CASE (1)
91
      flatdir(1:12) = [ 1,2, &
92
                        2,3, &
93
                        1,3, &
94
                        1,4, &
95
                        2,4, &
96
                        3,4 ]
97
    CASE (2)
98
      flatdir(1:12) = [ 1,2, &
99
                        3,2, &
100
                        1,3, &
101
                        1,4, &
102
                        2,4, &
103
                        3,4 ]      
104
    CASE DEFAULT
105
      CALL Fatal('H1Basis_GetTetraEdgeDirection','Unknown tetra type')
106
    END SELECT
107
  END SUBROUTINE H1Basis_GetTetraEdgeDirection
108
  
109
  SUBROUTINE H1Basis_GetFaceDirection(ecode, nfaces, globalind, direction)
110
    IMPLICIT NONE
111
    INTEGER, INTENT(IN) :: ecode, nfaces
112
    INTEGER, DIMENSION(:), POINTER CONTIG, INTENT(IN) :: globalind
113
    INTEGER, DIMENSION(H1Basis_MaxPElementFaceNodes, &
114
                       H1Basis_MaxPElementFaces), INTENT(INOUT) :: direction
115
    INTEGER :: face, tmp, tmp1, tmp2, i, minI
116

117
    CALL H1Basis_GetFaceMap(ecode, direction)
118

119
!DIR$ LOOP COUNT MAX=6
120
    DO face=1,nfaces
121
      IF (direction(H1Basis_MaxPElementFaceNodes, face) == 0) THEN
122
        ! Triangle face (last local index 0)
123
        
124
        ! Global face direction is consistent when the local indices are chosen such 
125
        ! that the global index is sorted
126

127
        ! Sort globalind(direction(:,face))
128
        ! (1,2) -> (1,2) or (2,1)
129
        IF (globalind(direction(1,face))>globalind(direction(2,face))) THEN
130
          tmp = direction(1,face)
131
          direction(1,face)=direction(2,face)
132
          direction(2,face)=tmp
133
        END IF
134
        ! (1,3) -> (1,3) or (3,1)
135
        IF (globalind(direction(1,face))>globalind(direction(3,face))) THEN
136
          tmp = direction(1,face)
137
          direction(1,face)=direction(3,face)
138
          direction(3,face)=tmp
139
        END IF
140
        ! (2,3) -> (2,3) or (3,2)
141
        IF (globalind(direction(2,face))>globalind(direction(3,face))) THEN
142
          tmp = direction(2,face)
143
          direction(2,face)=direction(3,face)
144
          direction(3,face)=tmp
145
        END IF
146
      ELSE
147
        ! Square face
148
        ! Find index of minimum global index (A)
149
        minI = 1
150
        DO i=2,4
151
          IF (globalind(direction(i,face))<globalind(direction(minI,face))) minI=i
152
        END DO
153
        
154
        ! In place rotate face to left or right cyclically to move minI as first element
155
        ! sqface(1:4)=cshift(direction(1:4,face),minI-1)
156
        SELECT CASE(minI-1)
157
        CASE(0)
158
          ! Nothing to do!
159
        CASE(1)
160
          tmp = direction(1,face)
161
          ! arr(i) -> arr(i+1)
162
          DO i=1,3
163
            direction(i,face)=direction(i+1,face)
164
          END DO
165
          direction(4,face)=tmp
166
        CASE(2)
167
          ! (1,3)->(3,1)
168
          tmp1 = direction(1,face)
169
          direction(1,face)=direction(3,face)
170
          direction(3,face)=tmp1
171
          ! (2,4)->(4,2)
172
          tmp2 = direction(2,face)
173
          direction(2,face)=direction(4,face)
174
          direction(4,face)=tmp2
175
        CASE(3)
176
          tmp = direction(4,face)
177
          ! arr(i+1) -> arr(i)
178
          DO i=3,1,-1
179
            direction(i+1,face)=direction(i,face)
180
          END DO
181
          direction(1,face)=tmp
182
        END SELECT
183

184
        ! Find second smallest index B next to global minimum (either 2 or 4)
185
        IF (globalind(direction(2,face))>globalind(direction(4,face))) THEN
186
          tmp=direction(2,face)
187
          direction(2,face)=direction(4,face)
188
          direction(4,face)=tmp
189
        END IF
190

191
        ! Let C be the remaining index next to the global minimum A and 
192
        ! D the index opposite of A -> [A,B,D,C] forms a globally consistent ordering
193
      END IF
194
    END DO
195
  END SUBROUTINE H1Basis_GetFaceDirection
196

197
  SUBROUTINE H1Basis_GetTetraFaceDirection(ttype, direction)
198
    IMPLICIT NONE
199
    INTEGER, INTENT(IN) :: ttype
200
    INTEGER, DIMENSION(H1Basis_MaxPElementFaceNodes, &
201
                       H1Basis_MaxPElementFaces), TARGET, INTENT(INOUT) :: direction
202
    ! Local variables
203
    INTEGER, DIMENSION(:), POINTER CONTIG :: flatdir
204

205
    ! Remap array bounds
206
    flatdir(1:H1Basis_MaxPElementFaceNodes*H1Basis_MaxPElementFaces) => direction
207
    SELECT CASE (ttype)
208
    CASE (1)
209
      flatdir(1:16) = [ 1,2,3,0, & 
210
                        1,2,4,0, &
211
                        2,3,4,0, &
212
                        1,3,4,0 ]
213
    CASE (2)
214
      flatdir(1:16) = [ 1,3,2,0, &
215
                        1,2,4,0, &
216
                        3,2,4,0, &
217
                        1,3,4,0 ]
218
    CASE DEFAULT
219
      CALL Fatal('H1Basis_GetTetraFaceDirection','Unknown tetra type')
220
    END SELECT
221
  END SUBROUTINE H1Basis_GetTetraFaceDirection
222

223
  SUBROUTINE H1Basis_GetEdgeMap(ecode, map)
224
    IMPLICIT NONE
225

226
    INTEGER, INTENT(IN) :: ecode
227
    INTEGER, DIMENSION(H1Basis_MaxPElementEdgeNodes, &
228
                       H1Basis_MaxPElementEdges), TARGET, INTENT(OUT) :: map
229
    ! Local variables
230
    INTEGER, DIMENSION(:), POINTER CONTIG :: flatmap
231

232
    ! Remap array bounds
233
    flatmap(1:H1Basis_MaxPElementEdgeNodes*H1Basis_MaxPElementEdges) => map
234

235
    SELECT CASE (ecode)
236
    CASE(202)
237
      ! Line edge mappings
238
      flatmap(1:2) = [ 1,2 ]
239
    CASE(303)
240
      ! Triangle edge mappings
241
      flatmap(1:6) = [ 1,2, &
242
                       2,3, &
243
                       3,1 ]
244
    CASE(404)
245
      ! Quad edge mappings
246
      flatmap(1:8) = [ 1,2,&
247
                       2,3,&
248
                       4,3,&
249
                       1,4 ]
250
    CASE(504)
251
      ! Tetra edge mapping (for completeness, not needed for enforcing parity)
252
      CALL H1Basis_GetTetraEdgeDirection(1, map)
253

254
    CASE(605)
255
      flatmap(1:16) = [ 1,2,&
256
                        2,3,&
257
                        4,3,&
258
                        1,4,&
259
                        1,5,&
260
                        2,5,&
261
                        3,5,&
262
                        4,5 ]
263

264
    CASE(706)
265
      flatmap(1:18) = [ 1,2,&
266
                        2,3,&
267
                        3,1,&
268
                        4,5,&
269
                        5,6,&
270
                        6,4,&
271
                        1,4,&
272
                        2,5,&
273
                        3,6 ]
274
    CASE(808)
275
      flatmap(1:24) = [ 1,2,&
276
                        2,3,&
277
                        4,3,&
278
                        1,4,&
279
                        5,6,&
280
                        6,7,&
281
                        8,7,&
282
                        5,8,&
283
                        1,5,&
284
                        2,6,&
285
                        3,7,&
286
                        4,8 ]
287

288
    CASE DEFAULT
289
      CALL Fatal('H1Basis_GetEdgeMap','Not fully implemented yet!')
290
    END SELECT
291
  END SUBROUTINE H1Basis_GetEdgeMap
292

293
  SUBROUTINE H1Basis_GetFaceMap(ecode, map)
294
    IMPLICIT NONE
295

296
    INTEGER, INTENT(IN) :: ecode
297
    INTEGER, DIMENSION(H1Basis_MaxPElementFaceNodes, &
298
                       H1Basis_MaxPElementFaces), TARGET, INTENT(OUT) :: map
299
    ! Local variables
300
    INTEGER, DIMENSION(:), POINTER CONTIG :: flatmap
301

302
    ! Remap array bounds
303
    flatmap(1:H1Basis_MaxPElementFaceNodes*H1Basis_MaxPElementFaces) => map
304

305
    SELECT CASE (ecode)
306
    CASE(303)
307
      ! 2D boundary element in a 3D mesh
308
      flatmap(1:4) = [ 1,2,3,0 ]
309
    CASE(404)
310
      ! 2D boundary element in a 3D mesh
311
      flatmap(1:4) = [ 1,2,3,4 ]
312
    CASE(504)
313
      ! Tetra face mapping (for completeness, not needed for enforcing parity)
314
      CALL H1Basis_GetTetraFaceDirection(1, map)
315

316
    CASE(605)
317
      ! Pyramid face mappings
318
      flatmap(1:20) = [ 1,2,3,4, &
319
                        1,2,5,0, &
320
                        2,3,5,0, &
321
                        3,4,5,0, &
322
                        4,1,5,0 ]
323

324
    CASE(706)
325
      flatmap(1:20) = [ 1,2,3,0, &
326
                        4,5,6,0, &
327
                        1,2,5,4, &
328
                        2,3,6,5, &
329
                        3,1,4,6 ]
330
    CASE(808)
331
      flatmap(1:24) = [ 1,2,3,4, &
332
                        5,6,7,8, &
333
                        1,2,6,5, &
334
                        2,3,7,6, &
335
                        4,3,7,8, &
336
                        1,4,8,5 ]
337
    CASE DEFAULT
338
      CALL Fatal('H1Basis_GetFaceMap','Not fully implemented yet!')
339
    END SELECT
340
  END SUBROUTINE H1Basis_GetFaceMap
341
  
342
  SUBROUTINE H1Basis_LineNodal(nvec, u, nbasismax, fval, nbasis)
343
    IMPLICIT NONE
344

345
    ! Parameters
346
    INTEGER, INTENT(IN) :: nvec
347
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u
348
    ! Variables
349
    INTEGER, INTENT(IN) :: nbasismax
350
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
351
    INTEGER, INTENT(INOUT) :: nbasis
352
    REAL(Kind=dp), PARAMETER :: c = 1D0/2D0
353
    INTEGER :: j
354
!DIR$ ASSUME_ALIGNED u:64, fval:64
355

356
    !_ELMER_OMP_SIMD
357
    DO j=1,nvec
358
      ! Node 1
359
      fval(j,nbasis+1) = c*(1-u(j))
360
      ! Node 2
361
      fval(j,nbasis+2) = c*(1+u(j))
362
    END DO
363
    nbasis = nbasis + 2
364
  END SUBROUTINE H1Basis_LineNodal
365

366
  SUBROUTINE H1Basis_dLineNodal(nvec, u, nbasismax, grad, nbasis)
367
    IMPLICIT NONE
368

369
    ! Parameters
370
    INTEGER, INTENT(IN) :: nvec
371
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u
372
    ! Variables
373
    INTEGER, INTENT(IN) :: nbasismax
374
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
375
    INTEGER, INTENT(INOUT) :: nbasis
376
    REAL(Kind=dp), PARAMETER :: c = 1D0/2D0
377
    INTEGER :: j
378
!DIR$ ASSUME_ALIGNED u:64, grad:64
379

380
    ! First coordinate (xi)
381
    !_ELMER_OMP_SIMD
382
    DO j=1,nvec
383
      grad(j,nbasis+1,1) = -c
384
      grad(j,nbasis+2,1) =  c
385
    END DO
386
    nbasis = nbasis + 2
387
  END SUBROUTINE H1Basis_dLineNodal
388

389
  SUBROUTINE H1Basis_LineBubbleP(nvec, u, pmax, nbasismax, fval, nbasis, invertEdge)
390
    IMPLICIT NONE
391

392
    INTEGER, INTENT(IN) :: nvec
393
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u
394
    INTEGER, INTENT(IN) :: pmax
395
    INTEGER, INTENT(IN) :: nbasismax
396
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
397
    INTEGER, INTENT(INOUT) :: nbasis
398
    LOGICAL, OPTIONAL :: invertEdge
399

400
    ! Local variables
401
    LOGICAL :: invert
402
    INTEGER :: p, j
403
!DIR$ ASSUME_ALIGNED u:64, fval:64
404

405
    ! Check if line basis has been inverted (not by default)
406
    invert = .FALSE.
407
    IF (PRESENT( invertEdge )) invert = invertEdge
408

409
    IF (invert) THEN
410
      DO p=1,pmax-1
411
        !_ELMER_OMP_SIMD 
412
        DO j=1,nvec
413
           fval(j,nbasis+p) = H1Basis_Phi(p+1,-u(j))
414
        END DO
415
      END DO
416
    ELSE
417
      DO p=1,pmax-1
418
        !_ELMER_OMP_SIMD 
419
        DO j=1,nvec
420
          fval(j,nbasis+p) = H1Basis_Phi(p+1,u(j))
421
        END DO
422
      END DO
423
    END IF
424

425
    ! nbasis = nbasis+(pmax-2)+1
426
    nbasis = nbasis+pmax-1
427
  END SUBROUTINE H1Basis_LineBubbleP
428

429
  SUBROUTINE H1Basis_dLineBubbleP(nvec, u, pmax, nbasismax, grad, nbasis, invertEdge)
430
    IMPLICIT NONE
431

432
    INTEGER, INTENT(IN) :: nvec
433
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u
434
    INTEGER, INTENT(IN) :: pmax
435
    INTEGER, INTENT(IN) :: nbasismax
436
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
437
    INTEGER, INTENT(INOUT) :: nbasis
438
    LOGICAL, OPTIONAL :: invertEdge
439

440
    ! Local variables
441
    LOGICAL :: invert
442
    INTEGER :: p, j
443
!DIR$ ASSUME_ALIGNED u:64, grad:64
444

445
    ! Check if line basis has been inverted (not by default)
446
    invert = .FALSE.
447
    IF (PRESENT( invertEdge )) invert = invertEdge
448

449
    IF (invert) THEN
450
      DO p=1,pmax-1
451
        ! First coordinate (xi)
452
        !_ELMER_OMP_SIMD
453
        DO j=1,nvec
454
          grad(j,nbasis+p,1) = H1Basis_dPhi(p+1,-u(j))
455
        END DO
456
      END DO
457
    ELSE
458
      DO p=1,pmax-1
459
        ! First coordinate (xi)
460
        !_ELMER_OMP_SIMD
461
        DO j=1,nvec
462
          grad(j,nbasis+p,1) = H1Basis_dPhi(p+1,u(j))
463
        END DO
464
      END DO
465
    END IF
466

467
    ! nbasis = nbasis+(pmax-2)+1
468
    nbasis = nbasis+pmax-1
469
  END SUBROUTINE H1Basis_dLineBubbleP
470

471
  SUBROUTINE H1Basis_TriangleNodalP(nvec, u, v, nbasismax, fval, nbasis)
472
    IMPLICIT NONE
473

474
    ! Parameters
475
    INTEGER, INTENT(IN) :: nvec
476
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v
477
    ! Variables
478
    INTEGER, INTENT(IN) :: nbasismax
479
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
480
    INTEGER, INTENT(INOUT) :: nbasis
481
    
482
    INTEGER :: j
483
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0, d = 1D0/SQRT(3D0)
484
!DIR$ ASSUME_ALIGNED u:64, v:64, fval:64
485

486
    !_ELMER_OMP_SIMD
487
    DO j=1,nvec
488
      fval(j,nbasis+1) = c*(1-u(j)-d*v(j))
489
      fval(j,nbasis+2) = c*(1+u(j)-d*v(j))
490
      fval(j,nbasis+3) = d*v(j)
491
    END DO
492
    nbasis = nbasis + 3
493
  END SUBROUTINE H1Basis_TriangleNodalP
494

495
  FUNCTION H1Basis_TriangleL(node, u, v) RESULT(fval)
496
    IMPLICIT NONE
497

498
    ! Parameters 
499
    INTEGER, INTENT(IN) :: node
500
    REAL(KIND=dp), INTENT(IN) :: u,v
501
    ! Result
502
    REAL(KIND=dp) :: fval
503
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0, d = 1D0/SQRT(3D0)
504
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) &
505
    !_ELMER_OMP _ELMER_LINEAR_REF(u) _ELMER_LINEAR_REF(v) NOTINBRANCH
506
    
507
    SELECT CASE(node)
508
    CASE(1)
509
      fval = c*(1-u-d*v)
510
    CASE(2)
511
      fval = c*(1+u-d*v)
512
    CASE(3)
513
      fval = d*v
514
    END SELECT
515
  END FUNCTION H1Basis_TriangleL
516

517
  SUBROUTINE H1Basis_dTriangleNodalP(nvec, u, v, nbasismax, grad, nbasis)
518
    IMPLICIT NONE
519

520
    ! Parameters
521
    INTEGER, INTENT(IN) :: nvec
522
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v
523
    ! Variables
524
    INTEGER, INTENT(IN) :: nbasismax
525
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
526
    INTEGER, INTENT(INOUT) :: nbasis
527
    INTEGER :: j
528
    REAL(KIND=dp), PARAMETER :: half=1._dp/2, a=1/SQRT(3._dp), b=-a/2
529
!DIR$ ASSUME_ALIGNED u:64, v:64, grad:64
530

531
    ! First coordinate (xi)
532
    !_ELMER_OMP_SIMD
533
    DO j=1,nvec
534
      grad(j,nbasis+1,1) = -half
535
    END DO
536
    !_ELMER_OMP_SIMD
537
    DO j=1,nvec
538
      grad(j,nbasis+2,1) =  half
539
    END DO
540
    !_ELMER_OMP_SIMD
541
    DO j=1,nvec
542
      grad(j,nbasis+3,1) = 0
543
    END DO
544
    ! Second coordinate (eta)
545
    !_ELMER_OMP_SIMD
546
    DO j=1,nvec
547
      grad(j,nbasis+1,2) = b
548
    END DO
549
    !_ELMER_OMP_SIMD
550
    DO j=1,nvec
551
      grad(j,nbasis+2,2) = b
552
    END DO
553
    !_ELMER_OMP_SIMD
554
    DO j=1,nvec
555
      grad(j,nbasis+3,2) = a
556
    END DO
557
    nbasis = nbasis + 3
558
  END SUBROUTINE H1Basis_dTriangleNodalP
559
  
560
  FUNCTION H1Basis_dTriangleL(node) RESULT(grad)
561
    IMPLICIT NONE
562

563
    ! Parameters 
564
    INTEGER, INTENT(IN) :: node
565
    ! REAL(KIND=dp), INTENT(IN) :: u,v
566
    ! Result
567
    REAL(KIND=dp) :: grad(2)
568
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0, d = 1D0/SQRT(3D0)
569
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) NOTINBRANCH
570
    
571
    SELECT CASE(node)
572
    CASE(1)
573
      grad = c*[REAL(-1,dp), REAL(-d,dp)]
574
    CASE(2)
575
      grad = c*[REAL(1,dp), REAL(-d,dp)]
576
    CASE(3)
577
      grad = [REAL(0,dp), REAL(d,dp)]
578
    END SELECT
579
  END FUNCTION H1Basis_dTriangleL
580
  
581
  SUBROUTINE H1Basis_TriangleEdgeP(nvec, u, v, pmax, nbasismax, fval, nbasis, edgedir)
582
    IMPLICIT NONE
583

584
    INTEGER, INTENT(IN) :: nvec
585
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
586
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
587
    INTEGER, INTENT(IN) :: nbasismax
588
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
589
    INTEGER, INTENT(INOUT) :: nbasis
590
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
591

592
    REAL(KIND=dp) :: La, Lb
593
    INTEGER :: i,j,k
594
!DIR$ ASSUME_ALIGNED u:64, v:64, fval:64
595

596
    ! For each edge
597
    DO i=1,3
598
      DO j=1,pmax(i)-1
599
        !_ELMER_OMP_SIMD PRIVATE(La, Lb)
600
        DO k=1,nvec
601
          La = H1Basis_TriangleL(edgedir(1,i),u(k),v(k))
602
          Lb = H1Basis_TriangleL(edgedir(2,i),u(k),v(k))
603
          fval(k, nbasis+j) = La*Lb*H1Basis_varPhi(j+1,Lb-La)
604
        END DO
605
      END DO
606
      ! nbasis = nbasis + (pmax(i)-2) + 1
607
      nbasis = nbasis + pmax(i) - 1
608
    END DO
609

610
  END SUBROUTINE H1Basis_TriangleEdgeP
611

612
  SUBROUTINE H1Basis_dTriangleEdgeP(nvec, u, v, pmax, nbasismax, grad, nbasis, edgedir)
613
    IMPLICIT NONE
614

615
    INTEGER, INTENT(IN) :: nvec
616
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
617
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
618
    INTEGER, INTENT(IN) :: nbasismax
619
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
620
    INTEGER, INTENT(INOUT) :: nbasis
621
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
622

623
    REAL(KIND=dp) :: La, Lb, vPhi, dVPhi, dLa(2), dLb(2)
624
    INTEGER :: i,j,k
625
!DIR$ ASSUME_ALIGNED u:64, v:64, grad:64
626

627
    ! For each edge
628
    DO i=1,3
629
      dLa = H1Basis_dTriangleL(edgedir(1,i))
630
      dLb = H1Basis_dTriangleL(edgedir(2,i))
631

632
      DO j=1,pmax(i)-1
633
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, vPhi, dVPhi)
634
        DO k=1,nvec
635
          La = H1Basis_TriangleL(edgedir(1,i),u(k),v(k))
636
          Lb = H1Basis_TriangleL(edgedir(2,i),u(k),v(k))
637
          vPhi = H1Basis_varPhi(j+1,Lb-La)
638
          dVPhi = H1Basis_dVarPhi(j+1,Lb-La)
639

640
          grad(k, nbasis+j, 1) = dLa(1)*Lb*vPhi+ &
641
                  La*dLb(1)*vPhi + La*Lb*dVPhi*(dLb(1)-dLa(1))
642
          grad(k, nbasis+j, 2) = dLa(2)*Lb*vPhi+ &
643
                  La*dLb(2)*vPhi + La*Lb*dVPhi*(dLb(2)-dLa(2))
644
        END DO
645
      END DO
646
      ! nbasis = nbasis + (pmax(i)-2) + 1
647
      nbasis = nbasis + pmax(i) - 1
648
    END DO
649

650
  END SUBROUTINE H1Basis_dTriangleEdgeP
651

652
  SUBROUTINE H1Basis_TriangleBubbleP(nvec, u, v, pmax, nbasismax, fval, nbasis, localnumbers)
653
    IMPLICIT NONE
654

655
    INTEGER, INTENT(IN) :: nvec
656
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
657
    INTEGER, INTENT(IN) :: pmax
658
    INTEGER, INTENT(IN) :: nbasismax
659
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
660
    INTEGER, INTENT(INOUT) :: nbasis
661
    INTEGER, INTENT(IN), DIMENSION(3), OPTIONAL :: localnumbers
662
      
663
    ! Variables
664
    INTEGER :: i,j,k
665
    REAL (KIND=dp) :: La, Lb, Lc
666
!DIR$ ASSUME_ALIGNED u:64, v:64, fval:64
667

668
    IF (.NOT. PRESENT(localnumbers)) THEN
669
      ! Triangle bubble basis
670
      DO i = 0,pmax-3
671
        DO j = 1,pmax-i-2
672
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc)
673
          DO k=1,nvec
674
            La = H1Basis_TriangleL(1,u(k),v(k))
675
            Lb = H1Basis_TriangleL(2,u(k),v(k))
676
            Lc = H1Basis_TriangleL(3,u(k),v(k))
677
          
678
            fval(k,nbasis+j) = La*Lb*Lc*(H1Basis_PowInt((Lb-La),i))*&
679
                    H1Basis_PowInt((2*Lc-1),j-1)
680
          END DO
681
        END DO
682
        ! nbasis = nbasis + (pmax-i-3) + 1
683
        nbasis = nbasis + MAX(pmax-i-2,0)
684
      END DO
685
    ELSE
686
      ! Triangular face basis
687
      DO i=0,pmax-3
688
        DO j=1,pmax-i-2
689
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc)
690
          DO k=1,nvec
691
            La=H1Basis_TriangleL(localnumbers(1),u(k),v(k))
692
            Lb=H1Basis_TriangleL(localnumbers(2),u(k),v(k))
693
            Lc=H1Basis_TriangleL(localnumbers(3),u(k),v(k))
694
            
695
            fval(k,nbasis+j) = La*Lb*Lc*H1Basis_LegendreP(i,Lb-La)* &
696
                                          H1Basis_LegendreP(j-1,2*Lc-1)
697
          END DO
698
        END DO
699
        ! nbasis = nbasis + (pmax-i-3 + 1)
700
        nbasis = nbasis + MAX(pmax-i-2,0)
701
      END DO
702
    END IF
703
  END SUBROUTINE H1Basis_TriangleBubbleP
704

705
  SUBROUTINE H1Basis_dTriangleBubbleP(nvec, u, v, pmax, nbasismax, grad, nbasis, localnumbers)
706
    IMPLICIT NONE
707

708
    INTEGER, INTENT(IN) :: nvec
709
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
710
    INTEGER, INTENT(IN) :: pmax
711
    INTEGER, INTENT(IN) :: nbasismax
712
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
713
    INTEGER, INTENT(INOUT) :: nbasis
714
    INTEGER, INTENT(IN), DIMENSION(3), OPTIONAL :: localnumbers
715

716
    ! Variables
717
    REAL (KIND=dp) :: La, Lb, Lc, Lc_1, Lb_Lai, Lc_1n, a, b, Lb_La,&
718
            dLa(2), dLb(2), dLc(2)
719
    INTEGER :: i,j,k
720
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0, d = -SQRT(3D0)/6, e = SQRT(3D0)/3
721
!DIR$ ASSUME_ALIGNED u:64, v:64, grad:64
722

723
    IF (.NOT. PRESENT(localnumbers)) THEN
724
      ! Triangle bubble basis
725
      DO i = 0,pmax-3
726
        DO j = 1,pmax-i-2
727
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Lb_Lai, Lc_1n)
728
          DO k=1,nvec
729
            La = H1Basis_TriangleL(1,u(k),v(k))
730
            Lb = H1Basis_TriangleL(2,u(k),v(k))
731
            Lc = H1Basis_TriangleL(3,u(k),v(k))
732
            
733
            Lb_Lai = H1Basis_PowInt((Lb-La), i)
734
            Lc_1n = H1Basis_PowInt((2*Lc-1), j-1)
735
            
736
            ! Calculate value of function from general form
737
            grad(k,nbasis+j,1) = -c*Lb*Lc*Lb_Lai*Lc_1n + La*c*Lc*Lb_Lai*Lc_1n + &
738
                    La*Lb*Lc*i*(H1Basis_PowInt((Lb-La),i-1))*Lc_1n
739
            grad(k,nbasis+j,2) = d*Lb*Lc*Lb_Lai*Lc_1n + La*d*Lc*Lb_Lai*Lc_1n + &
740
                    La*Lb*e*Lb_Lai*Lc_1n + &
741
                    La*Lb*Lc*Lb_Lai*(j-1)*(H1Basis_PowInt((2*Lc-1),j-2))*2*e
742
          END DO
743
        END DO
744
        ! nbasis = nbasis + (pmax-i-3) + 1
745
        nbasis = nbasis + MAX(pmax-i-2,0)
746
      END DO
747
    ELSE
748
      ! Triangular face basis
749
      dLa = H1Basis_dTriangleL(localnumbers(1))
750
      dLb = H1Basis_dTriangleL(localnumbers(2))
751
      dLc = H1Basis_dTriangleL(localnumbers(3))
752

753
      DO i=0,pmax-3
754
        DO j=1,pmax-i-2
755
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Lb_La, Lc_1, a, b)
756
          DO k=1,nvec
757
            La=H1Basis_TriangleL(localnumbers(1),u(k),v(k))
758
            Lb=H1Basis_TriangleL(localnumbers(2),u(k),v(k))
759
            Lc=H1Basis_TriangleL(localnumbers(3),u(k),v(k))
760
            Lb_La=Lb-La
761
            Lc_1=2*Lc-1
762
            a=H1Basis_LegendreP(i,Lb_La)
763
            b=H1Basis_LegendreP(j-1,Lc_1)
764

765
            grad(k,nbasis+j,1) = dLa(1)*Lb*Lc*a*b + La*dLb(1)*Lc*a*b + &
766
                    La*Lb*dLc(1)*a*b + &
767
                    La*Lb*Lc*H1Basis_dLegendreP(i,Lb_La)*(dLb(1)-dLa(1))*b + &
768
                    La*Lb*Lc*a*H1Basis_dLegendreP(j-1,Lc_1)*2*dLc(1)
769
            grad(k,nbasis+j,2) = dLa(2)*Lb*Lc*a*b + La*dLb(2)*Lc*a*b + &
770
                    La*Lb*dLc(2)*a*b + &
771
                    La*Lb*Lc*H1Basis_dLegendreP(i,Lb_La)*(dLb(2)-dLa(2))*b + &
772
                    La*Lb*Lc*a*H1Basis_dLegendreP(j-1,Lc_1)*2*dLc(2)
773
          END DO
774
        END DO
775
        ! nbasis = nbasis + (pmax-i-3 + 1)
776
        nbasis = nbasis + MAX(pmax-i-2,0)
777
      END DO
778
    END IF
779

780
  END SUBROUTINE H1Basis_dTriangleBubbleP
781

782
  FUNCTION H1Basis_PowInt(x,j) RESULT(powi)
783
    IMPLICIT NONE
784
    REAL(KIND=dp), INTENT(IN) :: x
785
    INTEGER, INTENT(IN) :: j
786
    
787
    INTEGER :: i
788
    REAL(KIND=dp) :: powi
789
    !_ELMER_OMP_DECLARE_SIMD _ELMER_LINEAR_REF(x) UNIFORM(j) NOTINBRANCH
790
    
791
    powi=REAL(1,dp)
792
    DO i=1,j
793
      powi=powi*x
794
    END DO
795

796
  END FUNCTION H1Basis_PowInt
797

798
  SUBROUTINE H1Basis_QuadNodal(nvec, u, v, nbasismax, fval, nbasis)
799
    IMPLICIT NONE
800

801
    ! Parameters
802
    INTEGER, INTENT(IN) :: nvec
803
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v
804
    ! Variables
805
    INTEGER, INTENT(IN) :: nbasismax
806
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
807
    INTEGER, INTENT(INOUT) :: nbasis
808
    REAL(Kind=dp), PARAMETER :: c = 1D0/4D0
809
    INTEGER :: j
810
!DIR$ ASSUME_ALIGNED u:64, v:64, fval:64
811

812
    !_ELMER_OMP_SIMD
813
    DO j=1,nvec
814
      ! Node 1
815
      fval(j,1) = c*(1-u(j))*(1-v(j))
816
      ! Node 2
817
      fval(j,2) = c*(1+u(j))*(1-v(j))
818
      ! Node 3
819
      fval(j,3) = c*(1+u(j))*(1+v(j))
820
      ! Node 4
821
      fval(j,4) = c*(1-u(j))*(1+v(j))
822
    END DO
823
    nbasis = nbasis + 4
824
  END SUBROUTINE H1Basis_QuadNodal
825

826
  SUBROUTINE H1Basis_dQuadNodal(nvec, u, v, nbasismax, grad, nbasis)
827
    IMPLICIT NONE
828

829
    ! Parameters
830
    INTEGER, INTENT(IN) :: nvec
831
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v
832
    ! Variables
833
    INTEGER, INTENT(IN) :: nbasismax
834
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
835
    INTEGER, INTENT(INOUT) :: nbasis
836
    REAL(Kind=dp), PARAMETER :: c = 1D0/4D0
837
    INTEGER :: j
838
!DIR$ ASSUME_ALIGNED u:64, v:64, grad:64
839

840
    ! First coordinate (xi)
841
    !_ELMER_OMP_SIMD
842
    DO j=1,nvec
843
      grad(j,nbasis+1,1) = -c*(1-v(j))
844
      grad(j,nbasis+2,1) =  c*(1-v(j))
845
      grad(j,nbasis+3,1) =  c*(1+v(j))
846
      grad(j,nbasis+4,1) = -c*(1+v(j))
847
    END DO
848
    
849
    ! Second coordinate (eta)
850
    !_ELMER_OMP_SIMD
851
    DO j=1,nvec
852
      grad(j,nbasis+1,2) = -c*(1-u(j))
853
      grad(j,nbasis+2,2) = -c*(1+u(j))
854
      grad(j,nbasis+3,2) =  c*(1+u(j))
855
      grad(j,nbasis+4,2) =  c*(1-u(j))
856
    END DO
857
    
858
    nbasis = nbasis + 4
859
  END SUBROUTINE H1Basis_dQuadNodal
860

861
! --- start serendipity quad
862

863

864
  SUBROUTINE H1Basis_SD_QuadEdgeP(nvec, u, v, pmax, nbasismax, fval, nbasis, edgedir)
865
    IMPLICIT NONE
866

867
    INTEGER, INTENT(IN) :: nvec
868
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
869
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
870
    INTEGER, INTENT(IN) :: nbasismax
871
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
872
    INTEGER, INTENT(INOUT) :: nbasis
873
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
874
    REAL(KIND=dp), PARAMETER :: c = 1/2D0
875

876
    REAL(KIND=dp) :: La, Lb
877
    INTEGER :: i,j,k
878
!DIR$ ASSUME_ALIGNED u:64, v:64, fval:64
879

880
    ! For each edge
881
    DO i=1,4
882
      DO j=1,pmax(i)-1
883
        !_ELMER_OMP_SIMD PRIVATE(La, Lb)
884
        DO k=1,nvec
885
          La = H1Basis_QuadL(edgedir(1,i), u(k), v(k))
886
          Lb = H1Basis_QuadL(edgedir(2,i), u(k), v(k))
887

888
          fval(k, nbasis+j) = c*(La+Lb-1)*H1Basis_Phi(j+1, Lb-La)
889
        END DO
890
      END DO
891
      ! nbasis = nbasis + (pmax(i)-2) + 1
892
      nbasis = nbasis + pmax(i) - 1
893
    END DO
894
  END SUBROUTINE H1Basis_SD_QuadEdgeP
895

896
  SUBROUTINE H1Basis_SD_dQuadEdgeP(nvec, u, v, pmax, nbasismax, grad, nbasis, edgedir)
897
    IMPLICIT NONE
898

899
    INTEGER, INTENT(IN) :: nvec
900
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
901
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
902
    INTEGER, INTENT(IN) :: nbasismax
903
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
904
    INTEGER, INTENT(INOUT) :: nbasis
905
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
906

907
    REAL(KIND=dp) :: La, Lb, Phi, dPhi, dLa(2), dLb(2)
908
    REAL(KIND=dp), PARAMETER :: c = 1/2D0
909
    INTEGER :: i,j,k
910
!DIR$ ASSUME_ALIGNED u:64, v:64, grad:64
911
    
912
    ! For each edge
913
    DO i=1,4
914
      dLa = H1Basis_dQuadL(edgedir(1,i))
915
      dLb = H1Basis_dQuadL(edgedir(2,i))
916

917
      DO j=1,pmax(i)-1
918
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Phi, dPhi)
919
        DO k=1,nvec
920
          La = H1Basis_QuadL(edgedir(1,i), u(k), v(k))
921
          Lb = H1Basis_QuadL(edgedir(2,i), u(k), v(k))
922
          
923
          Phi = H1Basis_Phi(j+1, Lb-La)
924
          dPhi = H1Basis_dPhi(j+1,Lb-La)
925

926
          grad(k, nbasis+j, 1) = c*((dLa(1)+dLb(1))*Phi + &
927
                  (La+Lb-1)*dPhi*(dLb(1)-dLa(1)))
928
          grad(k, nbasis+j, 2) = c*((dLa(2)+dLb(2))*Phi + &
929
                  (La+Lb-1)*dPhi*(dLb(2)-dLa(2)))
930
        END DO
931
      END DO
932
      ! nbasis = nbasis + (pmax(i)-2) + 1
933
      nbasis = nbasis + pmax(i) - 1
934
    END DO
935
  END SUBROUTINE H1Basis_SD_dQuadEdgeP
936

937
  SUBROUTINE H1Basis_SD_QuadBubbleP(nvec, u, v, pmax, nbasismax, fval, nbasis, localNumbers)
938
    IMPLICIT NONE
939

940
    INTEGER, INTENT(IN) :: nvec
941
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
942
    INTEGER, INTENT(IN) :: pmax
943
    INTEGER, INTENT(IN) :: nbasismax
944
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
945
    INTEGER, INTENT(INOUT) :: nbasis
946
    INTEGER, INTENT(IN), OPTIONAL :: localNumbers(4)
947

948
    INTEGER :: i,j,k
949
    REAL(KIND=dp) :: La, Lb, Lc
950
!DIR$ ASSUME_ALIGNED u:64, v:64, fval:64
951

952
    ! Calculate value of function without direction and return
953
    ! if local numbering not present
954
    IF (.NOT. PRESENT(localNumbers)) THEN
955
      DO i=2,(pmax-2)
956
        DO j=1,(pmax-i)-1
957
          !_ELMER_OMP_SIMD
958
          DO k=1,nvec
959
            fval(k,nbasis+j) = H1Basis_Phi(i,u(k))*H1Basis_Phi(j+1,v(k))
960
          END DO
961
        END DO
962
        nbasis = nbasis+MAX(pmax-i-1,0)
963
      END DO
964
    ELSE
965
      DO i=2,(pmax-2)
966
        DO j=1,(pmax-i)-1
967
          !_ELMER_OMP_SIMD PRIVATE(La,Lb,Lc)
968
          DO k=1,nvec
969
            ! Directed quad bubbles
970
            La = H1Basis_QuadL(localNumbers(1),u(k),v(k))
971
            Lb = H1Basis_QuadL(localNumbers(2),u(k),v(k))
972
            Lc = H1Basis_QuadL(localNumbers(4),u(k),v(k))
973

974
            ! Calculate value of function from the general form
975
            fval(k,nbasis+j) = H1Basis_Phi(i,Lb-La)*H1Basis_Phi(j+1,Lc-La)
976
          END DO
977
        END DO
978
        nbasis = nbasis+MAX(pmax-i-1,0)
979
      END DO
980
    END IF
981
  END SUBROUTINE H1Basis_SD_QuadBubbleP
982

983
  SUBROUTINE H1Basis_SD_dQuadBubbleP(nvec, u, v, pmax, nbasismax, grad, nbasis, localNumbers)
984
    IMPLICIT NONE
985

986
    INTEGER, INTENT(IN) :: nvec
987
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
988
    INTEGER, INTENT(IN) :: pmax
989
    INTEGER, INTENT(IN) :: nbasismax
990
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
991
    INTEGER, INTENT(INOUT) :: nbasis
992
    INTEGER, INTENT(IN), OPTIONAL :: localNumbers(4)
993

994
    INTEGER :: i,j,k
995
    REAL(KIND=dp) :: La, Lb, Lc
996
    REAL(Kind=dp), DIMENSION(2) :: dLa, dLb, dLc, dLbdLa, dLcdLa
997
!DIR$ ASSUME_ALIGNED u:64, v:64, grad:64
998

999
    IF (.NOT. PRESENT(localNumbers)) THEN
1000
      DO i=2,(pmax-2)
1001
        DO j=1,(pmax-i)-1
1002
          !_ELMER_OMP_SIMD
1003
          DO k=1,nvec
1004
            ! First coordinate (Xi)
1005
            grad(k,nbasis+j,1) = H1Basis_dPhi(i,u(k))*H1Basis_Phi(j+1,v(k))
1006
            ! Second coordinate (Eta)
1007
            grad(k,nbasis+j,2) = H1Basis_Phi(i,u(k))*H1Basis_dPhi(j+1,v(k))
1008
          END DO
1009
        END DO
1010
        nbasis = nbasis+MAX((pmax-i)-1,0)
1011
      END DO
1012
    ELSE
1013
      ! Numbering present, so use it
1014
      dLa = H1Basis_dQuadL(localNumbers(1))
1015
      dLb = H1Basis_dQuadL(localNumbers(2))
1016
      dLc = H1Basis_dQuadL(localNumbers(4))
1017

1018
      dLBdLa = dLb-dLa
1019
      dLcdLa = dLc-dLa
1020

1021
      DO i=2,(pmax-2)
1022
        DO j=1,(pmax-i)-1
1023
          !_ELMER_OMP_SIMD PRIVATE(La,Lb,Lc)
1024
          DO k=1,nvec
1025
            La = H1Basis_QuadL(localNumbers(1),u(k),v(k))
1026
            Lb = H1Basis_QuadL(localNumbers(2),u(k),v(k))
1027
            Lc = H1Basis_QuadL(localNumbers(4),u(k),v(k))
1028

1029
            grad(k,nbasis+j,1) = H1Basis_dPhi(i,Lb-La)*(dLbdLa(1))*H1Basis_Phi(j+1,Lc-La) + &
1030
                    H1Basis_Phi(i,Lb-La)*H1Basis_dPhi(j+1,Lc-La)*(dLcdLa(1))
1031
            grad(k,nbasis+j,2) = H1Basis_dPhi(i,Lb-La)*(dLbdLa(2))*H1Basis_Phi(j+1,Lc-La) + &
1032
                    H1Basis_Phi(i,Lb-La)*H1Basis_dPhi(j+1,Lc-La)*(dLcdLa(2))
1033
          END DO
1034
        END DO
1035
        nbasis = nbasis+MAX((pmax-i)-1,0)
1036
      END DO
1037
    END IF
1038
  END SUBROUTINE H1Basis_SD_dQuadBubbleP
1039

1040

1041
! --- end serendipity quad
1042

1043

1044
  SUBROUTINE H1Basis_QuadEdgeP(nvec, u, v, pmax, nbasismax, fval, nbasis, edgedir)
1045
    IMPLICIT NONE
1046

1047
    INTEGER, INTENT(IN) :: nvec
1048
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
1049
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
1050
    INTEGER, INTENT(IN) :: nbasismax
1051
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
1052
    INTEGER, INTENT(INOUT) :: nbasis
1053
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
1054
    REAL(KIND=dp), PARAMETER :: c = 1/2D0*2
1055

1056
    REAL(KIND=dp) :: La, Lb, N(VECTOR_BLOCK_LENGTH,nbasismax), Na, Nb
1057
    INTEGER :: i,j,k, nnb, node1, node2
1058
!DIR$ ASSUME_ALIGNED u:64, v:64, fval:64
1059

1060
    ! For each edge
1061
    DO i=1,4
1062
      node1 = edgedir(1,i)
1063
      node2 = edgedir(2,i)
1064
      DO j=1,pmax(i)-1
1065
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb)
1066
        DO k=1,nvec
1067
          La = H1Basis_QuadL(node1, u(k), v(k))
1068
          Lb = H1Basis_QuadL(node2, u(k), v(k))
1069

1070
          Na = fval(k,node1)
1071
          Nb = fval(k,node2)
1072

1073
          fval(k, nbasis+j) = c*Na*Nb*H1Basis_varPhi(j+1, Lb-La)
1074
        END DO
1075
      END DO
1076
      ! nbasis = nbasis + (pmax(i)-2) + 1
1077
      nbasis = nbasis + pmax(i) - 1
1078
    END DO
1079
  END SUBROUTINE H1Basis_QuadEdgeP
1080

1081
  SUBROUTINE H1Basis_dQuadEdgeP(nvec, u, v, pmax, nbasismax, grad, nbasis, edgedir)
1082
    IMPLICIT NONE
1083

1084
    INTEGER, INTENT(IN) :: nvec
1085
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
1086
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
1087
    INTEGER, INTENT(IN) :: nbasismax
1088
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
1089
    INTEGER, INTENT(INOUT) :: nbasis
1090
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
1091

1092
    REAL(KIND=dp) :: La, Lb, Phi, dPhi, dLa(2), dLb(2), N(VECTOR_BLOCK_LENGTH,nbasismax), &
1093
           dN(VECTOR_BLOCK_LENGTH,nbasismax,3), Na, Nb, dNa(2), dNb(2)
1094
    REAL(KIND=dp), PARAMETER :: c = 1/2D0*2
1095
    INTEGER :: i,j,k,nnb, node1, node2
1096
!DIR$ ASSUME_ALIGNED u:64, v:64, grad:64
1097
    
1098
    nnb = 0; CALL H1Basis_QuadNodal(nvec, u, v, nbasismax, n, nnb )
1099

1100
    ! For each edge
1101
    DO i=1,4
1102
      node1 = edgedir(1,i)
1103
      node2 = edgedir(2,i)
1104
      dLa = H1Basis_dQuadL(node1)
1105
      dLb = H1Basis_dQuadL(node2)
1106

1107
      DO j=1,pmax(i)-1
1108
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Phi, dPhi, Na, Nb, dNa, dNb)
1109
        DO k=1,nvec
1110
          La = H1Basis_QuadL(node1, u(k), v(k))
1111
          Lb = H1Basis_QuadL(node2, u(k), v(k))
1112
          
1113
          Na  = N(k,node1)
1114
          Nb  = N(k,node2)
1115
          dNa = grad(k,node1,1:2)
1116
          dNb = grad(k,node2,1:2)
1117

1118
          Phi = H1Basis_varPhi(j+1, Lb-La)
1119
          dPhi = H1Basis_dvarPhi(j+1,Lb-La)
1120

1121
          grad(k, nbasis+j, 1:2) = c*(dNa*Nb*Phi + Na*dNb*Phi + &
1122
                    Na*Nb*dPhi*(dLb-dLa))
1123
        END DO
1124
      END DO
1125
      ! nbasis = nbasis + (pmax(i)-2) + 1
1126
      nbasis = nbasis + pmax(i) - 1
1127
    END DO
1128
  END SUBROUTINE H1Basis_dQuadEdgeP
1129

1130
  SUBROUTINE H1Basis_QuadBubbleP(nvec, u, v, pmax, nbasismax, fval, nbasis, localNumbers)
1131
    IMPLICIT NONE
1132

1133
    INTEGER, INTENT(IN) :: nvec
1134
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
1135
    INTEGER, INTENT(IN) :: pmax
1136
    INTEGER, INTENT(IN) :: nbasismax
1137
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
1138
    INTEGER, INTENT(INOUT) :: nbasis
1139
    INTEGER, INTENT(IN), OPTIONAL :: localNumbers(4)
1140

1141
    INTEGER :: i,j,k
1142
    REAL(KIND=dp) :: La, Lb, Lc, Pa, Pb
1143
!DIR$ ASSUME_ALIGNED u:64, v:64, fval:64
1144

1145

1146
    ! Calculate value of function without direction and return
1147
    ! if local numbering not present
1148
    IF (.NOT. PRESENT(localNumbers)) THEN
1149
      DO i=2,pmax
1150
        DO j=1,pmax-1
1151
          !_ELMER_OMP_SIMD
1152
          DO k=1,nvec
1153
            fval(k,nbasis+j) = H1Basis_Phi(i,u(k))*H1Basis_Phi(j+1,v(k))
1154
          END DO
1155
        END DO
1156
        nbasis = nbasis + pmax - 1
1157
      END DO
1158
    ELSE
1159
      DO i=0,pmax-2
1160
        DO j=1,pmax-1
1161
          !_ELMER_OMP_SIMD PRIVATE(La,Lb,Lc,Pa,Pb)
1162
          DO k=1,nvec
1163
            ! Directed quad bubbles
1164
            La = H1Basis_QuadL(localNumbers(1),u(k),v(k))
1165
            Lb = H1Basis_QuadL(localNumbers(2),u(k),v(k))
1166
            Lc = H1Basis_QuadL(localNumbers(4),u(k),v(k))
1167

1168
            Pa = fval(k,localNumbers(1))
1169
            Pb = fval(k,localNumbers(3))
1170

1171
            ! Calculate value of function from the general form
1172
            fval(k,nbasis+j) = Pa*Pb*H1Basis_LegendreP(i,Lb-La)*H1Basis_LegendreP(j-1,Lc-La)
1173
          END DO
1174
        END DO
1175
        nbasis = nbasis + pmax - 1
1176
      END DO
1177
    END IF
1178
  END SUBROUTINE H1Basis_QuadBubbleP
1179

1180
  SUBROUTINE H1Basis_dQuadBubbleP(nvec, u, v, pmax, nbasismax, grad, nbasis, localNumbers)
1181
    IMPLICIT NONE
1182

1183
    INTEGER, INTENT(IN) :: nvec
1184
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v
1185
    INTEGER, INTENT(IN) :: pmax
1186
    INTEGER, INTENT(IN) :: nbasismax
1187
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
1188
    INTEGER, INTENT(INOUT) :: nbasis
1189
    INTEGER, INTENT(IN), OPTIONAL :: localNumbers(4)
1190

1191
    INTEGER :: i,j,k, nnb
1192
    REAL(KIND=dp) :: La, Lb, Lc, Legi, Legj, Pa, Pb
1193
    REAL(Kind=dp), DIMENSION(2) :: dLa, dLb, dLc, dLbdLa, dLcdLa, dLegi,dLegj, dPa,dPb
1194
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax) :: N
1195
!DIR$ ASSUME_ALIGNED u:64, v:64, grad:64
1196

1197
    IF (.NOT. PRESENT(localNumbers)) THEN
1198
      DO i=2,pmax
1199
        DO j=1,pmax-1
1200
          !_ELMER_OMP_SIMD
1201
          DO k=1,nvec
1202
            ! First coordinate (Xi)
1203
            grad(k,nbasis+j,1) = H1Basis_dPhi(i,u(k))*H1Basis_Phi(j+1,v(k))
1204
            ! Second coordinate (Eta)
1205
            grad(k,nbasis+j,2) = H1Basis_Phi(i,u(k))*H1Basis_dPhi(j+1,v(k))
1206
          END DO
1207
        END DO
1208
        nbasis = nbasis + pmax - 1
1209
      END DO
1210
    ELSE
1211
      nnb = 0; CALL H1Basis_QuadNodal(nvec, u, v, nbasismax, n, nnb )
1212

1213
      ! Numbering present, so use it
1214
      dLa = H1Basis_dQuadL(localNumbers(1))
1215
      dLb = H1Basis_dQuadL(localNumbers(2))
1216
      dLc = H1Basis_dQuadL(localNumbers(4))
1217

1218
      DO i=0,pmax-2
1219
        DO j=1,pmax-1
1220
          !_ELMER_OMP_SIMD PRIVATE(La,Lb,Lc,Legi,Legj,dLegi,dLegj,Pa,Pb,dPa,dPb)
1221
          DO k=1,nvec
1222
            La = H1Basis_QuadL(localNumbers(1),u(k),v(k))
1223
            Lb = H1Basis_QuadL(localNumbers(2),u(k),v(k))
1224
            Lc = H1Basis_QuadL(localNumbers(4),u(k),v(k))
1225

1226
            Pa = N(k,localNumbers(1))
1227
            Pb = N(k,localNumbers(3))
1228

1229
            dPa = grad(k,localNumbers(1),1:2)
1230
            dPb = grad(k,localNumbers(3),1:2)
1231

1232
            Legi = H1Basis_LegendreP(i,Lb-La)
1233
            Legj = H1Basis_LegendreP(j-1,Lc-La)
1234

1235
            dLegi = H1Basis_dLegendreP(i,Lb-La)*(dLb-dLa)
1236
            dLegj = H1Basis_dLegendreP(j-1,Lc-La)*(dLc-dLa)
1237

1238
            grad(k,nbasis+j,1:2) = dPa*Pb*Legi*Legj + Pa*dPb*Legi*Legj + &
1239
                                   Pa*Pb*dLegi*Legj + Pa*Pb*Legi*dLegj
1240
          END DO
1241
        END DO
1242
        nbasis = nbasis + pmax - 1
1243
      END DO
1244
    END IF
1245
  END SUBROUTINE H1Basis_dQuadBubbleP
1246

1247
  FUNCTION H1Basis_QuadL(node, u, v) RESULT(fval)
1248
    IMPLICIT NONE
1249

1250
    INTEGER, INTENT(IN) :: node
1251
    REAL(KIND=dp), INTENT(IN) :: u, v
1252
    REAL(KIND=dp) :: fval
1253
    REAL(KIND=dp), PARAMETER :: c = 1/2D0
1254
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) &
1255
    !_ELMER_OMP _ELMER_LINEAR_REF(u) _ELMER_LINEAR_REF(v) NOTINBRANCH
1256
    
1257
    SELECT CASE (node)
1258
    CASE (1)
1259
      fval = c*(2-u-v)
1260
    CASE (2)
1261
      fval = c*(2+u-v)
1262
    CASE (3)
1263
      fval = c*(2+u+v)
1264
    CASE (4)
1265
      fval = c*(2-u+v)
1266
    END SELECT
1267
  END FUNCTION H1Basis_QuadL
1268

1269
  FUNCTION H1Basis_dQuadL(node) RESULT(grad)
1270
    IMPLICIT NONE
1271

1272
    INTEGER, INTENT(IN) :: node
1273
    ! REAL(KIND=dp), INTENT(IN) :: u, v
1274
    REAL(KIND=dp) :: grad(2)
1275
    REAL(KIND=dp), PARAMETER :: c = 1/2D0
1276
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) NOTINBRANCH
1277
    
1278
    SELECT CASE(node)
1279
    CASE (1)
1280
      grad(1:2) = c*[-1,-1 ]
1281
    CASE (2)
1282
      grad(1:2) = c*[ 1,-1 ]
1283
    CASE (3)
1284
      grad(1:2) = c*[ 1, 1 ]
1285
    CASE (4)
1286
      grad(1:2) = c*[-1, 1 ]
1287
    END SELECT
1288
  END FUNCTION H1Basis_dQuadL
1289

1290
  SUBROUTINE H1Basis_TetraNodalP(nvec, u, v, w, nbasismax, fval, nbasis)
1291
    IMPLICIT NONE
1292

1293
    ! Parameters
1294
    INTEGER, INTENT(IN) :: nvec
1295
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
1296
    ! Variables
1297
    INTEGER, INTENT(IN) :: nbasismax
1298
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
1299
    INTEGER, INTENT(INOUT) :: nbasis
1300
    INTEGER :: j
1301

1302
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0, d = 1D0/SQRT(3D0), &
1303
            e = 1D0/SQRT(6D0), f = 1D0/SQRT(8D0), g = SQRT(3D0)*f
1304
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
1305

1306
    !_ELMER_OMP_SIMD
1307
    DO j=1,nvec
1308
      ! Node 1
1309
      fval(j,nbasis+1) = c*(1-u(j)-d*v(j)-e*w(j))
1310
      ! Node 2
1311
      fval(j,nbasis+2) = c*(1+u(j)-d*v(j)-e*w(j))
1312
      ! Node 3
1313
      fval(j,nbasis+3) = d*(v(j)-f*w(j))
1314
      ! Node 4
1315
      fval(j,nbasis+4) = g*w(j)
1316
    END DO
1317
    nbasis = nbasis + 4
1318
  END SUBROUTINE H1Basis_TetraNodalP
1319

1320
  FUNCTION H1Basis_TetraL(node, u, v, w) RESULT(fval)
1321
    IMPLICIT NONE
1322

1323
    ! Parameters 
1324
    INTEGER, INTENT(IN) :: node
1325
    REAL(KIND=dp), INTENT(IN) :: u,v,w
1326
    ! Result
1327
    REAL(KIND=dp) :: fval
1328
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0, d = 1D0/SQRT(3D0), &
1329
            e = 1D0/SQRT(6D0), f = 1D0/SQRT(8D0), g = SQRT(3D0)*f
1330
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) &
1331
    !_ELMER_OMP _ELMER_LINEAR_REF(u) _ELMER_LINEAR_REF(v) _ELMER_LINEAR_REF(w) NOTINBRANCH
1332
    
1333
    SELECT CASE(node)
1334
    CASE(1)
1335
      fval = c*(1-u-d*v-e*w)
1336
    CASE(2)
1337
      fval = c*(1+u-d*v-e*w)
1338
    CASE(3)
1339
      fval = d*(v-f*w)
1340
    CASE(4)
1341
      fval = g*w      
1342
    END SELECT
1343
  END FUNCTION H1Basis_TetraL
1344
  
1345
  SUBROUTINE H1Basis_dTetraNodalP(nvec, u, v, w, nbasismax, grad, nbasis)
1346
    IMPLICIT NONE
1347

1348
    ! Parameters
1349
    INTEGER, INTENT(IN) :: nvec
1350
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
1351
    ! Variables
1352
    INTEGER, INTENT(IN) :: nbasismax
1353
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
1354
    INTEGER, INTENT(INOUT) :: nbasis
1355
    INTEGER :: j
1356
    REAL(KIND=dp), PARAMETER :: half=1d0/2, a = 1/SQRT(3._dp), b=-a/2, c = -SQRT(6._dp)/12, &
1357
                  d = -3*c
1358
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
1359

1360
    ! First coordinate (xi)
1361
    !_ELMER_OMP_SIMD
1362
    DO j=1,nvec
1363
      grad(j,nbasis+1,1) =-half
1364
    END DO
1365
    !_ELMER_OMP_SIMD
1366
    DO j=1,nvec
1367
      grad(j,nbasis+2,1) = half
1368
    END DO
1369
    !_ELMER_OMP_SIMD
1370
    DO j=1,nvec
1371
      grad(j,nbasis+3,1) = 0
1372
    END DO
1373
    !_ELMER_OMP_SIMD
1374
    DO j=1,nvec
1375
      grad(j,nbasis+4,1) = 0
1376
    END DO
1377
    
1378
    ! Second coordinate (eta)
1379
    !_ELMER_OMP_SIMD
1380
    DO j=1,nvec
1381
      grad(j,nbasis+1,2) = b
1382
    END DO
1383
    !_ELMER_OMP_SIMD
1384
    DO j=1,nvec
1385
      grad(j,nbasis+2,2) = b
1386
    END DO
1387
    !_ELMER_OMP_SIMD
1388
    DO j=1,nvec
1389
      grad(j,nbasis+3,2) = a
1390
    END DO
1391
    !_ELMER_OMP_SIMD
1392
    DO j=1,nvec
1393
      grad(j,nbasis+4,2) = 0
1394
    END DO
1395
    
1396
    ! Third coordinate (zeta)
1397
    !_ELMER_OMP_SIMD
1398
    DO j=1,nvec
1399
      grad(j,nbasis+1,3) = c
1400
    END DO
1401
    !_ELMER_OMP_SIMD
1402
    DO j=1,nvec
1403
      grad(j,nbasis+2,3) = c
1404
    END DO
1405
    !_ELMER_OMP_SIMD
1406
    DO j=1,nvec
1407
      grad(j,nbasis+3,3) = c
1408
    END DO
1409
    !_ELMER_OMP_SIMD
1410
    DO j=1,nvec
1411
      grad(j,nbasis+4,3) = d
1412
    END DO
1413
    nbasis = nbasis + 4
1414
  END SUBROUTINE H1Basis_dTetraNodalP
1415
  
1416
  FUNCTION H1Basis_dTetraL(node) RESULT(grad)
1417
    IMPLICIT NONE
1418

1419
    ! Parameters 
1420
    INTEGER, INTENT(IN) :: node
1421
    ! Result
1422
    REAL(KIND=dp) :: grad(3)
1423
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0, d = 1D0/SQRT(3D0), &
1424
            e = 1D0/SQRT(6D0), f = 1D0/SQRT(8D0), g = SQRT(3._dp)*f
1425
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) NOTINBRANCH
1426
    
1427
    SELECT CASE(node)
1428
    CASE(1)
1429
      grad = c*[-1D0, -d, -e]
1430
    CASE(2)
1431
      grad = c*[1D0, -d, -e]
1432
    CASE(3)
1433
      grad = d*[0D0, 1D0, -f]
1434
    CASE(4)
1435
      grad = [0D0, 0D0, g]
1436
    END SELECT
1437
  END FUNCTION H1Basis_dTetraL
1438

1439
  SUBROUTINE H1Basis_TetraEdgeP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, edgedir)
1440
    IMPLICIT NONE
1441

1442
    INTEGER, INTENT(IN) :: nvec
1443
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
1444
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
1445
    INTEGER, INTENT(IN) :: nbasismax
1446
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
1447
    INTEGER, INTENT(INOUT) :: nbasis
1448
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
1449

1450
    REAL(KIND=dp) :: La, Lb
1451
    INTEGER :: i,j,k
1452
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
1453

1454
    ! For each edge
1455
    DO i=1,6
1456
      DO j=1,pmax(i)-1
1457
        !_ELMER_OMP_SIMD PRIVATE(La, Lb)
1458
        DO k=1,nvec
1459
          La = H1Basis_TetraL(edgedir(1,i), u(k), v(k), w(k))
1460
          Lb = H1Basis_TetraL(edgedir(2,i), u(k), v(k), w(k))
1461
          fval(k, nbasis+j) = La*Lb*H1Basis_varPhi(j+1,Lb-La)
1462
        END DO
1463
      END DO
1464
      ! nbasis = nbasis + (pmax(i)-2) + 1
1465
      nbasis = nbasis + pmax(i) - 1
1466
    END DO
1467
  END SUBROUTINE H1Basis_TetraEdgeP
1468
  
1469
  SUBROUTINE H1Basis_dTetraEdgeP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, edgedir)
1470
    IMPLICIT NONE
1471

1472
    INTEGER, INTENT(IN) :: nvec
1473
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
1474
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
1475
    INTEGER, INTENT(IN) :: nbasismax
1476
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
1477
    INTEGER, INTENT(INOUT) :: nbasis
1478
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
1479

1480
    REAL(KIND=dp) :: La, Lb, vPhi, dVPhi, dLa(3), dLb(3)
1481
    INTEGER :: i,j,k
1482
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
1483

1484
    ! For each edge
1485
    DO i=1,6
1486
      dLa = H1Basis_dTetraL(edgedir(1,i))
1487
      dLb = H1Basis_dTetraL(edgedir(2,i))
1488

1489
      DO j=1,pmax(i)-1
1490
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, vPhi, dVPhi)
1491
        DO k=1,nvec
1492
          La = H1Basis_TetraL(edgedir(1,i),u(k),v(k),w(k))
1493
          Lb = H1Basis_TetraL(edgedir(2,i),u(k),v(k),w(k))
1494
          vPhi = H1Basis_varPhi(j+1,Lb-La)
1495
          dVPhi = H1Basis_dVarPhi(j+1,Lb-La)
1496

1497
          grad(k, nbasis+j, 1) = dLa(1)*Lb*vPhi + &
1498
                  La*dLb(1)*vPhi + La*Lb*dVPhi*(dLb(1)-dLa(1))
1499
          grad(k, nbasis+j, 2) = dLa(2)*Lb*vPhi+ &
1500
                  La*dLb(2)*vPhi + La*Lb*dVPhi*(dLb(2)-dLa(2))
1501
          grad(k, nbasis+j, 3) = dLa(3)*Lb*vPhi+ &
1502
                  La*dLb(3)*vPhi + La*Lb*dVPhi*(dLb(3)-dLa(3))
1503
        END DO
1504
      END DO
1505
      ! nbasis = nbasis + (pmax(i)-2) + 1
1506
      nbasis = nbasis + pmax(i) - 1
1507
    END DO
1508
  END SUBROUTINE H1Basis_dTetraEdgeP
1509

1510
  SUBROUTINE H1Basis_TetraFaceP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, facedir)
1511
    IMPLICIT NONE
1512

1513
    INTEGER, INTENT(IN) :: nvec
1514
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
1515
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
1516
    INTEGER, INTENT(IN) :: nbasismax
1517
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
1518
    INTEGER, INTENT(INOUT) :: nbasis
1519
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
1520

1521
    REAL(KIND=dp) :: La, Lb, Lc
1522
    INTEGER :: i,j,k,l
1523
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
1524

1525
    ! For each face
1526
    DO i=1,4
1527
      DO j=0,pmax(i)-3
1528
        DO k=1,pmax(i)-j-2
1529
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc)
1530
          DO l=1,nvec
1531
            La=H1Basis_TetraL(facedir(1,i),u(l),v(l),w(l))
1532
            Lb=H1Basis_TetraL(facedir(2,i),u(l),v(l),w(l))
1533
            Lc=H1Basis_TetraL(facedir(3,i),u(l),v(l),w(l))
1534
    
1535
            fval(l,nbasis+k) = La*Lb*Lc*H1Basis_LegendreP(j,Lb-La)* &
1536
                                        H1Basis_LegendreP(k-1,2*Lc-1)
1537
          END DO
1538
        END DO
1539
        ! nbasis = nbasis + (pmax(i)-j-3 + 1)
1540
        nbasis = nbasis + MAX(pmax(i)-j-2,0)
1541
      END DO
1542
    END DO
1543
  END SUBROUTINE H1Basis_TetraFaceP
1544
  
1545
  SUBROUTINE H1Basis_dTetraFaceP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, facedir)
1546
    IMPLICIT NONE
1547

1548
    INTEGER, INTENT(IN) :: nvec
1549
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
1550
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
1551
    INTEGER, INTENT(IN) :: nbasismax
1552
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
1553
    INTEGER, INTENT(INOUT) :: nbasis
1554
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
1555

1556
    REAL(KIND=dp) :: La, Lb, Lc, dLa(3), dLb(3), dLc(3), Lb_La, Lc_1, a, b
1557
    INTEGER :: i,j,k,l
1558
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
1559

1560
    ! For each face
1561
    DO i=1,4
1562
      dLa = H1Basis_dTetraL(facedir(1,i))
1563
      dLb = H1Basis_dTetraL(facedir(2,i))
1564
      dLc = H1Basis_dTetraL(facedir(3,i))
1565

1566
      DO j=0,pmax(i)-3
1567
        DO k=1,pmax(i)-j-2
1568
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Lb_La, Lc_1, a, b)
1569
          DO l=1,nvec
1570
            La=H1Basis_TetraL(facedir(1,i),u(l),v(l),w(l))
1571
            Lb=H1Basis_TetraL(facedir(2,i),u(l),v(l),w(l))
1572
            Lc=H1Basis_TetraL(facedir(3,i),u(l),v(l),w(l))
1573
            Lb_La=Lb-La
1574
            Lc_1=2*Lc-1
1575
            a=H1Basis_LegendreP(j,Lb_La)
1576
            b=H1Basis_LegendreP(k-1,Lc_1)
1577

1578
            grad(l,nbasis+k,1) = dLa(1)*Lb*Lc*a*b + La*dLb(1)*Lc*a*b + &
1579
                    La*Lb*dLc(1)*a*b + &
1580
                    La*Lb*Lc*H1Basis_dLegendreP(j,Lb_La)*(dLb(1)-dLa(1))*b + &
1581
                    La*Lb*Lc*a*H1Basis_dLegendreP(k-1,Lc_1)*2*dLc(1)
1582
            grad(l,nbasis+k,2) = dLa(2)*Lb*Lc*a*b + La*dLb(2)*Lc*a*b + &
1583
                    La*Lb*dLc(2)*a*b + &
1584
                    La*Lb*Lc*H1Basis_dLegendreP(j,Lb_La)*(dLb(2)-dLa(2))*b + &
1585
                    La*Lb*Lc*a*H1Basis_dLegendreP(k-1,Lc_1)*2*dLc(2)
1586
            grad(l,nbasis+k,3) = dLa(3)*Lb*Lc*a*b + La*dLb(3)*Lc*a*b + &
1587
                    La*Lb*dLc(3)*a*b + &
1588
                    La*Lb*Lc*H1Basis_dLegendreP(j,Lb_La)*(dLb(3)-dLa(3))*b + &
1589
                    La*Lb*Lc*a*H1Basis_dLegendreP(k-1,Lc_1)*2*dLc(3)
1590
          END DO
1591
        END DO
1592
        ! nbasis = nbasis + (pmax(i)-j-3 + 1)
1593
        nbasis = nbasis + MAX(pmax(i)-j-2,0)
1594
      END DO
1595
    END DO
1596
  END SUBROUTINE H1Basis_dTetraFaceP
1597
  
1598
  SUBROUTINE H1Basis_TetraBubbleP(nvec, u, v, w, pmax, nbasismax, fval, nbasis)
1599
    IMPLICIT NONE
1600

1601
    INTEGER, INTENT(IN) :: nvec
1602
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
1603
    INTEGER, INTENT(IN) :: pmax
1604
    INTEGER, INTENT(IN) :: nbasismax
1605
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
1606
    INTEGER, INTENT(INOUT) :: nbasis
1607
    
1608
    ! Variables
1609
    INTEGER :: i, j, k, l
1610
    ! Variables
1611
    REAL(KIND=dp) :: L1, L2, L3, L4, L2_L1, L3_1, L4_1
1612
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
1613

1614
    DO i=0,pmax-4
1615
      DO j=0,pmax-i-4
1616
        DO k=1,pmax-i-j-3
1617
          !_ELMER_OMP_SIMD PRIVATE(L1,L2,L3,L4,L2_L1,L3_1,L4_1)
1618
          DO l=1,nvec
1619
            L1 = H1Basis_TetraL(1,u(l),v(l),w(l))
1620
            L2 = H1Basis_TetraL(2,u(l),v(l),w(l))
1621
            L3 = H1Basis_TetraL(3,u(l),v(l),w(l))
1622
            L4 = H1Basis_TetraL(4,u(l),v(l),w(l))
1623
            L2_L1 = L2-L1
1624
            L3_1 = 2*L3-1
1625
            L4_1 = 2*L4-1
1626

1627
            fval(l,nbasis+k) = L1*L2*L3*L4*&
1628
                    H1Basis_LegendreP(i,L2_L1)*&
1629
                    H1Basis_LegendreP(j,L3_1)*&
1630
                    H1Basis_LegendreP(k-1,L4_1)
1631
          END DO
1632
        END DO
1633
        ! nbasis = nbasis + (pmax-i-j-4) + 1
1634
        nbasis = nbasis + MAX(pmax-i-j-3,0)
1635
      END DO
1636
    END DO
1637
  END SUBROUTINE H1Basis_TetraBubbleP
1638

1639
  SUBROUTINE H1Basis_dTetraBubbleP(nvec, u, v, w, pmax, nbasismax, grad, nbasis)
1640
    IMPLICIT NONE
1641

1642
    INTEGER, INTENT(IN) :: nvec
1643
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
1644
    INTEGER, INTENT(IN) :: pmax
1645
    INTEGER, INTENT(IN) :: nbasismax
1646
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
1647
    INTEGER, INTENT(INOUT) :: nbasis
1648

1649
    ! Variables
1650
    INTEGER :: i, j, k, l
1651
    REAL(KIND=dp) :: L1, L2, L3, L4, L2_L1, L3_1, L4_1, a, b, c 
1652
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
1653

1654
    DO i=0,pmax-4
1655
      DO j=0,pmax-i-4
1656
        DO k=1,pmax-i-j-3
1657
          !_ELMER_OMP_SIMD PRIVATE(L1,L2,L3,L4,L2_L1,L3_1,L4_1,a,b,c)
1658
          DO l=1,nvec
1659
            L1 = H1Basis_TetraL(1,u(l),v(l),w(l))
1660
            L2 = H1Basis_TetraL(2,u(l),v(l),w(l))
1661
            L3 = H1Basis_TetraL(3,u(l),v(l),w(l))
1662
            L4 = H1Basis_TetraL(4,u(l),v(l),w(l))
1663
            L2_L1 = L2 - L1
1664
            L3_1 = 2*L3 - 1
1665
            L4_1 = 2*L4 - 1
1666
            a = H1Basis_LegendreP(i,L2_L1)
1667
            b = H1Basis_LegendreP(j,L3_1)
1668
            c = H1Basis_LegendreP(k-1,L4_1)
1669

1670
            ! Gradients of tetrahedral bubble basis functions 
1671
            grad(l,nbasis+k,1) = -1d0/2*L2*L3*L4*a*b*c + &
1672
                    1d0/2*L1*L3*L4*a*b*c + &
1673
                    L1*L2*L3*L4*H1Basis_dLegendreP(i,L2_L1)*b*c
1674
            grad(l,nbasis+k,2) = -SQRT(3d0)/6*L2*L3*L4*a*b*c - &
1675
                    SQRT(3d0)/6*L1*L3*L4*a*b*c + &
1676
                    SQRT(3d0)/3*L1*L2*L4*a*b*c &
1677
                    + 2*SQRT(3d0)/3*L1*L2*L3*L4*a*H1Basis_dLegendreP(j,L3_1)*c
1678
            grad(l,nbasis+k,3) = -SQRT(6d0)/12*L2*L3*L4*a*b*c - &
1679
                    SQRT(6d0)/12*L1*L3*L4*a*b*c - &
1680
                    SQRT(6d0)/12*L1*L2*L4*a*b*c &
1681
                    + SQRT(6d0)/4*L1*L2*L3*a*b*c - &
1682
                    SQRT(6d0)/6*L1*L2*L3*L4*a*H1Basis_dLegendreP(j,L3_1)*c &
1683
                    + SQRT(6d0)/2*L1*L2*L3*L4*a*b*H1Basis_dLegendreP(k-1,L4_1)
1684
          END DO
1685
        END DO
1686
        ! nbasis = nbasis + (pmax-i-j-4) + 1
1687
        nbasis = nbasis + MAX(pmax-i-j-3,0)
1688
      END DO
1689
    END DO
1690
  END SUBROUTINE H1Basis_dTetraBubbleP
1691

1692

1693
! --- start serendipity wedge
1694

1695
  
1696
  SUBROUTINE H1Basis_SD_WedgeEdgeP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, edgedir)
1697
    IMPLICIT NONE
1698

1699
    INTEGER, INTENT(IN) :: nvec
1700
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
1701
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
1702
    INTEGER, INTENT(IN) :: nbasismax
1703
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
1704
    INTEGER, INTENT(INOUT) :: nbasis
1705
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
1706

1707
    REAL(KIND=dp) :: La, Lb, Na, Nb
1708
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0
1709
    INTEGER :: i,j,k,l
1710
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
1711

1712
    ! For each triangle face edge
1713
    DO i=1,6
1714
      DO j=1,pmax(i)-1
1715
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb)
1716
        DO k=1,nvec
1717
          La = H1Basis_WedgeL(edgedir(1,i), u(k), v(k))
1718
          Lb = H1Basis_WedgeL(edgedir(2,i), u(k), v(k))
1719
          Na = H1Basis_WedgeH(edgedir(1,i), w(k))
1720
          Nb = H1Basis_WedgeH(edgedir(2,i), w(k))
1721
          fval(k, nbasis+j) = c*La*Lb*H1Basis_varPhi(j+1,Lb-La)*(1+Na+Nb)
1722
        END DO
1723
      END DO
1724
      ! nbasis = nbasis + (pmax(i)-2) + 1
1725
      nbasis = nbasis + pmax(i) - 1
1726
    END DO
1727
      
1728
    ! For each square face edge
1729
    DO i=7,9
1730
      DO j=1,pmax(i)-1
1731
        !_ELMER_OMP_SIMD PRIVATE(La, Na, Nb)
1732
        DO k=1,nvec
1733
          La = H1Basis_WedgeL(edgedir(1,i), u(k), v(k))
1734
          Na = H1Basis_WedgeH(edgedir(1,i), w(k))
1735
          Nb = H1Basis_WedgeH(edgedir(2,i), w(k))
1736
          fval(k, nbasis+j) = La*H1Basis_Phi(j+1, Nb-Na)
1737
        END DO
1738
      END DO
1739
      ! nbasis = nbasis + (pmax(i)-2) + 1
1740
      nbasis = nbasis + pmax(i) - 1
1741
    END DO
1742
  END SUBROUTINE H1Basis_SD_WedgeEdgeP
1743
  
1744
  SUBROUTINE H1Basis_SD_dWedgeEdgeP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, edgedir)
1745
    IMPLICIT NONE
1746

1747
    INTEGER, INTENT(IN) :: nvec
1748
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
1749
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
1750
    INTEGER, INTENT(IN) :: nbasismax
1751
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
1752
    INTEGER, INTENT(INOUT) :: nbasis
1753
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
1754

1755
    REAL(KIND=dp) :: La, Lb, Na, Nb, Phi, dPhi, vPhi, dVPhi, &
1756
          dLa(3), dLb(3), dNa(3), dNb(3), NaNb
1757
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0
1758
    INTEGER :: i,j,k
1759
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
1760

1761
    ! For each triangle face edge
1762
    DO i=1,6
1763
      dLa = H1Basis_dWedgeL(edgedir(1,i))
1764
      dLb = H1Basis_dWedgeL(edgedir(2,i))
1765
      dNa = H1Basis_dWedgeH(edgedir(1,i))
1766
      dNb = H1Basis_dWedgeH(edgedir(2,i))
1767
      DO j=1,pmax(i)-1
1768
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb, vPhi, dVPhi, NaNb)
1769
        DO k=1,nvec
1770
          La = H1Basis_WedgeL(edgedir(1,i), u(k), v(k))
1771
          Lb = H1Basis_WedgeL(edgedir(2,i), u(k), v(k))
1772
          Na = H1Basis_WedgeH(edgedir(1,i), w(k))
1773
          Nb = H1Basis_WedgeH(edgedir(2,i), w(k))
1774
          vPhi=H1Basis_varPhi(j+1,Lb-La)
1775
          dVPhi=H1Basis_dVarPhi(j+1,Lb-La)
1776
          NaNb=1+Na+Nb
1777
          grad(k, nbasis+j,1) = c*dLa(1)*Lb*vPhi*NaNb + &
1778
                c*La*dLb(1)*vPhi*NaNb + c*La*Lb*dVPhi*(dLb(1)-dLa(1))*NaNb
1779
          grad(k, nbasis+j,2) = c*dLa(2)*Lb*vPhi*NaNb + &
1780
                c*La*dLb(2)*vPhi*NaNb + c*La*Lb*dVPhi*(dLb(2)-dLa(2))*NaNb
1781
          grad(k, nbasis+j,3) = c*La*Lb*vPhi*(dNb(3)+dNa(3))
1782
        END DO
1783
      END DO
1784
      ! nbasis = nbasis + (pmax(i)-2) + 1
1785
      nbasis = nbasis + pmax(i) - 1
1786
    END DO
1787
      
1788
    ! For each square face edge
1789
    DO i=7,9
1790
      dLa = H1Basis_dWedgeL(edgedir(1,i))
1791
      dNa = H1Basis_dWedgeH(edgedir(1,i))
1792
      dNb = H1Basis_dWedgeH(edgedir(2,i))
1793
      DO j=1,pmax(i)-1
1794
        !_ELMER_OMP_SIMD PRIVATE(La, Na, Nb, Phi, dPhi)
1795
        DO k=1,nvec
1796
          La = H1Basis_WedgeL(edgedir(1,i), u(k), v(k))
1797
          Na = H1Basis_WedgeH(edgedir(1,i), w(k))
1798
          Nb = H1Basis_WedgeH(edgedir(2,i), w(k))
1799
          Phi = H1Basis_Phi(j+1, Nb-Na)
1800
          dPhi = H1Basis_dPhi(j+1,Nb-Na)
1801
          
1802
          grad(k, nbasis+j,1) = dLa(1)*Phi+La*dPhi*(dNb(1)-dNa(1))
1803
          grad(k, nbasis+j,2) = dLa(2)*Phi+La*dPhi*(dNb(2)-dNa(2))
1804
          grad(k, nbasis+j,3) = dLa(3)*Phi+La*dPhi*(dNb(3)-dNa(3)) 
1805
        END DO
1806
      END DO
1807
      ! nbasis = nbasis + (pmax(i)-2) + 1
1808
      nbasis = nbasis + pmax(i) - 1
1809
    END DO
1810
  END SUBROUTINE H1Basis_SD_dWedgeEdgeP
1811
  
1812
  SUBROUTINE H1Basis_SD_WedgeFaceP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, facedir)
1813
    IMPLICIT NONE
1814

1815
    INTEGER, INTENT(IN) :: nvec
1816
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
1817
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
1818
    INTEGER, INTENT(IN) :: nbasismax
1819
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
1820
    INTEGER, INTENT(INOUT) :: nbasis
1821
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
1822

1823
    REAL(KIND=dp) :: La, Lb, Lc, Na, Nc
1824
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0
1825
    INTEGER :: i,j,k,l
1826
    LOGICAL :: nonpermuted
1827
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
1828

1829
    ! Triangle faces
1830
    DO i=1,2      
1831
      DO j=0,pmax(i)-3
1832
        DO k=1,pmax(i)-j-2
1833
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Na)
1834
          DO l=1,nvec
1835
            La = H1Basis_WedgeL(facedir(1,i), u(l), v(l))
1836
            Lb = H1Basis_WedgeL(facedir(2,i), u(l), v(l))
1837
            Lc = H1Basis_WedgeL(facedir(3,i), u(l), v(l))
1838
            Na = H1Basis_WedgeH(facedir(1,i), w(l))
1839

1840
            fval(l,nbasis+k) = c*(1+2*Na)*La*Lb*Lc* &
1841
                                            H1Basis_LegendreP(j, Lb-La)* &
1842
                                            H1Basis_LegendreP(k-1, 2*Lc-1)
1843
          END DO
1844
        END DO
1845
        ! nbasis = nbasis + (p-j-3) + 1
1846
        nbasis = nbasis + MAX(pmax(i)-j-2,0)
1847
      END DO
1848
    END DO
1849
    ! Square faces
1850
    DO i=3,5
1851
      ! First and second node must form an edge in upper or lower triangle
1852
      nonpermuted = (facedir(1,i) >= 1 .AND. facedir(1,i) <= 3 .AND.&
1853
                     facedir(2,i) >= 1 .AND. facedir(2,i) <= 3) .OR. &
1854
                    (facedir(1,i) >= 4 .AND. facedir(1,i) <= 6 .AND.&
1855
                     facedir(2,i) >= 4 .AND. facedir(2,i) <= 6)
1856
      
1857
      DO j=2,pmax(i)-2
1858
        DO k=1,pmax(i)-j-1
1859
          IF (nonpermuted) THEN
1860
            !_ELMER_OMP_SIMD PRIVATE(La,Lb,Na,Nc)
1861
            DO l=1,nvec
1862
              La = H1Basis_WedgeL(facedir(1,i), u(l), v(l))
1863
              Lb = H1Basis_WedgeL(facedir(2,i), u(l), v(l))
1864
              Na = H1Basis_WedgeH(facedir(1,i), w(l))
1865
              Nc = H1Basis_WedgeH(facedir(4,i), w(l))
1866
              fval(l,nbasis+k) = La*Lb*H1Basis_varPhi(j, Lb-La)* &
1867
                                         H1Basis_Phi(k+1, Nc-Na)
1868
            END DO
1869
          ELSE
1870
            !_ELMER_OMP_SIMD PRIVATE(La,Lb,Na,Nc)
1871
            DO l=1,nvec
1872
              ! Numbering permuted, switch indices j,k and nodes 2 and 4
1873
              La = H1Basis_WedgeL(facedir(1,i), u(l), v(l))
1874
              Lb = H1Basis_WedgeL(facedir(4,i), u(l), v(l))
1875
              Na = H1Basis_WedgeH(facedir(1,i), w(l))
1876
              Nc = H1Basis_WedgeH(facedir(2,i), w(l))
1877
              fval(l,nbasis+k) = La*Lb*H1Basis_varPhi(k+1, Lb-La)* &
1878
                                         H1Basis_Phi(j, Nc-Na)
1879
            END DO
1880
          END IF
1881
        END DO
1882
        ! nbasis = nbasis + (pmax(i)-j-2) + 1
1883
        nbasis = nbasis + MAX(pmax(i)-j-1,0)
1884
      END DO
1885
    END DO
1886
  END SUBROUTINE H1Basis_SD_WedgeFaceP
1887
  
1888
  SUBROUTINE H1Basis_SD_dWedgeFaceP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, facedir)
1889
    IMPLICIT NONE
1890

1891
    INTEGER, INTENT(IN) :: nvec
1892
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
1893
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
1894
    INTEGER, INTENT(IN) :: nbasismax
1895
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
1896
    INTEGER, INTENT(INOUT) :: nbasis
1897
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
1898

1899
    REAL(KIND=dp) :: La, Lb, Lc, Na, Nc, LegPLbLa, LegP2Lc1, &
1900
            cNa, vPhi, Phi, dLa(3), dLb(3), dLc(3), dNa(3), dNc(3)
1901
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0
1902
    INTEGER :: i,j,k,l
1903
    LOGICAL :: nonpermuted
1904
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
1905

1906
    ! Triangle faces
1907
    DO i=1,2
1908
      dLa = H1Basis_dWedgeL(facedir(1,i))
1909
      dLb = H1Basis_dWedgeL(facedir(2,i))
1910
      dLc = H1Basis_dWedgeL(facedir(3,i))
1911
      dNa = H1Basis_dWedgeH(facedir(1,i))
1912
      DO j=0,pmax(i)-3
1913
        DO k=1,pmax(i)-j-2
1914
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Na, LegPLbLa, LegP2Lc1, cNa)
1915
          DO l=1,nvec
1916
            La = H1Basis_WedgeL(facedir(1,i), u(l), v(l))
1917
            Lb = H1Basis_WedgeL(facedir(2,i), u(l), v(l))
1918
            Lc = H1Basis_WedgeL(facedir(3,i), u(l), v(l))
1919
            Na = H1Basis_WedgeH(facedir(1,i), w(l))
1920

1921
            LegPLbLa = H1Basis_LegendreP(j, Lb-La)
1922
            LegP2Lc1 = H1Basis_LegendreP(k-1, 2*Lc-1)
1923
            cNa = c*(1+2*Na)
1924

1925
            grad(l,nbasis+k,1) = cNa*dLa(1)*Lb*Lc*LegPLbLa*LegP2Lc1 + &
1926
                    cNa*La*dLb(1)*Lc*LegPLbLa*LegP2Lc1 + &
1927
                    cNa*La*Lb*dLc(1)*LegPLbLa*LegP2Lc1 + &
1928
                    cNa*La*Lb*Lc*H1Basis_dLegendreP(j, Lb-La)*(dLb(1)-dLa(1))*LegP2Lc1 + &
1929
                    cNa*La*Lb*Lc*LegPLbLa*H1Basis_dLegendreP(k-1, 2*Lc-1)*(2*dLc(1))
1930
            grad(l,nbasis+k,2) = cNa*dLa(2)*Lb*Lc*LegPLbLa*LegP2Lc1 + &
1931
                    cNa*La*dLb(2)*Lc*LegPLbLa*LegP2Lc1 + &
1932
                    cNa*La*Lb*dLc(2)*LegPLbLa*LegP2Lc1 + &
1933
                    cNa*La*Lb*Lc*H1Basis_dLegendreP(j, Lb-La)*(dLb(2)-dLa(2))*LegP2Lc1 + &
1934
                    cNa*La*Lb*Lc*LegPLbLa*H1Basis_dLegendreP(k-1, 2*Lc-1)*(2*dLc(2))
1935
            grad(l,nbasis+k,3) = 2*c*dNa(3)*La*Lb*Lc*LegPLbLa*LegP2Lc1
1936
          END DO
1937
        END DO
1938
        ! nbasis = nbasis + (p-j-3) + 1
1939
        nbasis = nbasis + MAX(pmax(i)-j-2,0)
1940
      END DO
1941
    END DO
1942
    ! Square faces
1943
    DO i=3,5
1944
      ! First and second node must form an edge in upper or lower triangle
1945
      nonpermuted = (facedir(1,i) >= 1 .AND. facedir(1,i) <= 3 .AND.&
1946
                     facedir(2,i) >= 1 .AND. facedir(2,i) <= 3) .OR. &
1947
                    (facedir(1,i) >= 4 .AND. facedir(1,i) <= 6 .AND.&
1948
                     facedir(2,i) >= 4 .AND. facedir(2,i) <= 6)
1949
      
1950
      IF (nonpermuted) THEN
1951
        dLa = H1Basis_dWedgeL(facedir(1,i))
1952
        dLb = H1Basis_dWedgeL(facedir(2,i))
1953
        dNa = H1Basis_dWedgeH(facedir(1,i))
1954
        dNc = H1Basis_dWedgeH(facedir(4,i))
1955
      ELSE
1956
        dLa = H1Basis_dWedgeL(facedir(1,i))
1957
        dLb = H1Basis_dWedgeL(facedir(4,i))
1958
        dNa = H1Basis_dWedgeH(facedir(1,i))
1959
        dNc = H1Basis_dWedgeH(facedir(2,i))
1960
      END IF
1961

1962
      DO j=2,pmax(i)-2
1963
        DO k=1,pmax(i)-j-1
1964
          IF (nonpermuted) THEN
1965
            !_ELMER_OMP_SIMD PRIVATE(La,Lb,Na,Nc,vPhi,Phi)
1966
            DO l=1,nvec
1967
              La = H1Basis_WedgeL(facedir(1,i), u(l), v(l))
1968
              Lb = H1Basis_WedgeL(facedir(2,i), u(l), v(l))
1969
              Na = H1Basis_WedgeH(facedir(1,i), w(l))
1970
              Nc = H1Basis_WedgeH(facedir(4,i), w(l))
1971
              
1972
              vPhi = H1Basis_varPhi(j, Lb-La)
1973
              Phi = H1Basis_Phi(k+1, Nc-Na)
1974

1975
              grad(l,nbasis+k,1) = dLa(1)*Lb*vPhi*Phi + La*dLb(1)*vPhi*Phi + &
1976
                      La*Lb*H1Basis_dVarPhi(j, Lb-La)*(dLb(1)-dLa(1))*Phi
1977
              grad(l,nbasis+k,2) = dLa(2)*Lb*vPhi*Phi + La*dLb(2)*vPhi*Phi + &
1978
                      La*Lb*H1Basis_dVarPhi(j, Lb-La)*(dLb(2)-dLa(2))*Phi
1979
              grad(l,nbasis+k,3) = La*Lb*vPhi*H1Basis_dPhi(k+1, Nc-Na)*(dNc(3)-dNa(3))
1980
            END DO
1981
          ELSE
1982
            !_ELMER_OMP_SIMD PRIVATE(La,Lb,Na,Nc,vPhi,Phi)
1983
            DO l=1,nvec
1984
              ! Numbering permuted, switch indices j,k and nodes 2 and 4
1985
              La = H1Basis_WedgeL(facedir(1,i), u(l), v(l))
1986
              Lb = H1Basis_WedgeL(facedir(4,i), u(l), v(l))
1987
              Na = H1Basis_WedgeH(facedir(1,i), w(l))
1988
              Nc = H1Basis_WedgeH(facedir(2,i), w(l))
1989
              
1990
              vPhi = H1Basis_varPhi(k+1, Lb-La)
1991
              Phi = H1Basis_Phi(j, Nc-Na)
1992

1993
              grad(l,nbasis+k,1) = dLa(1)*Lb*vPhi*Phi + La*dLb(1)*vPhi*Phi + &
1994
                      La*Lb*H1Basis_dVarPhi(k+1, Lb-La)*(dLb(1)-dLa(1))*Phi
1995
              grad(l,nbasis+k,2) = dLa(2)*Lb*vPhi*Phi + La*dLb(2)*vPhi*Phi + &
1996
                      La*Lb*H1Basis_dVarPhi(k+1, Lb-La)*(dLb(2)-dLa(2))*Phi
1997
              grad(l,nbasis+k,3) = La*Lb*vPhi*H1Basis_dPhi(j, Nc-Na)*(dNc(3)-dNa(3))
1998
            END DO
1999
          END IF
2000
        END DO
2001
        ! nbasis = nbasis + (pmax(i)-j-2) + 1
2002
        nbasis = nbasis + MAX(pmax(i)-j-1,0)
2003
      END DO
2004
    END DO
2005
  END SUBROUTINE H1Basis_SD_dWedgeFaceP
2006
  
2007
  SUBROUTINE H1Basis_SD_WedgeBubbleP(nvec, u, v, w, pmax, nbasismax, fval, nbasis)
2008
    IMPLICIT NONE
2009

2010
    INTEGER, INTENT(IN) :: nvec
2011
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
2012
    INTEGER, INTENT(IN) :: pmax
2013
    INTEGER, INTENT(IN) :: nbasismax
2014
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
2015
    INTEGER, INTENT(INOUT) :: nbasis
2016
    
2017
    INTEGER :: i, j, k, l
2018
    REAL(KIND=dp) :: L1, L2, L3, L2_L1, L3_1
2019
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
2020

2021
    DO i=0,pmax-5
2022
      DO j=0,pmax-5-i
2023
        DO k=1,pmax-4-i-j
2024
          !_ELMER_OMP_SIMD PRIVATE(L1, L2, L3, L2_L1, L3_1)
2025
          DO l=1,nvec
2026
            L1 = H1Basis_WedgeL(1,u(l),v(l))
2027
            L2 = H1Basis_WedgeL(2,u(l),v(l))
2028
            L3 = H1Basis_WedgeL(3,u(l),v(l))
2029
            L2_L1 = L2-L1
2030
            L3_1 = 2*L3-1
2031

2032
            ! Get value of bubble function
2033
            fval(l,nbasis+k) = L1*L2*L3*H1Basis_LegendreP(i,L2_L1)*&
2034
                    H1Basis_LegendreP(j,L3_1)*H1Basis_Phi(k+1,w(l))
2035
          END DO
2036
        END DO
2037
        ! nbasis = nbasis + (pmax-3-i-j)-2+1
2038
        nbasis = nbasis + MAX(pmax-4-i-j,0)
2039
      END DO
2040
    END DO
2041
  END SUBROUTINE H1Basis_SD_WedgeBubbleP
2042

2043
  SUBROUTINE H1Basis_SD_dWedgeBubbleP(nvec, u, v, w, pmax, nbasismax, grad, nbasis)
2044
    IMPLICIT NONE
2045

2046
    INTEGER, INTENT(IN) :: nvec
2047
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
2048
    INTEGER, INTENT(IN) :: pmax
2049
    INTEGER, INTENT(IN) :: nbasismax
2050
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
2051
    INTEGER, INTENT(INOUT) :: nbasis
2052
    
2053
    ! Parameters
2054
    INTEGER :: i,j,k,l
2055
    ! Variables
2056
    REAL(KIND=dp) :: L1,L2,L3,Legi,Legj,phiW,L2_L1,L3_1
2057
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64      
2058

2059
    DO i=0,pmax-5
2060
      DO j=0,pmax-5-i
2061
        DO k=1,pmax-4-i-j
2062
          !_ELMER_OMP_SIMD PRIVATE(L1, L2, L3, Legi, Legj, phiW, L2_L1, L3_1)
2063
          DO l=1,nvec
2064

2065
            ! Values of function L
2066
            L1 = H1Basis_WedgeL(1,u(l),v(l))
2067
            L2 = H1Basis_WedgeL(2,u(l),v(l))
2068
            L3 = H1Basis_WedgeL(3,u(l),v(l))
2069
            
2070
            L2_L1 = L2-L1
2071
            L3_1 = 2d0*L3-1
2072

2073
            Legi = H1Basis_LegendreP(i,L2_L1)
2074
            Legj = H1Basis_LegendreP(j,L3_1)
2075
            phiW = H1Basis_Phi(k+1,w(l))
2076
            
2077
            grad(l,nbasis+k,1) = ((-1d0/2)*L2*L3*Legi*Legj + &
2078
                    L1*(1d0/2)*L3*Legi*Legj +&
2079
                    L1*L2*L3*H1Basis_dLegendreP(i,L2_L1)*Legj)*phiW
2080
            grad(l,nbasis+k,2) = ((-SQRT(3d0)/6)*L2*L3*Legi*Legj + &
2081
                    L1*(-SQRT(3d0)/6)*L3*Legi*Legj +&
2082
                    L1*L2*(SQRT(3D0)/3)*Legi*Legj + &
2083
                    L1*L2*L3*Legi*H1Basis_dLegendreP(j,L3_1)*(2d0*SQRT(3D0)/3))*phiW 
2084
            grad(l,nbasis+k,3) = L1*L2*L3*Legi*Legj*H1Basis_dPhi(k+1,w(l))
2085
          END DO
2086
        END DO
2087
        ! nbasis = nbasis + (pmax-3-i-j)-2+1
2088
        nbasis = nbasis + MAX(pmax-4-i-j,0)
2089
      END DO
2090
    END DO
2091
  END SUBROUTINE H1Basis_SD_dWedgeBubbleP
2092

2093

2094
! --- end serendipity wedge
2095

2096

2097
  SUBROUTINE H1Basis_WedgeNodalP(nvec, u, v, w, nbasismax, fval, nbasis)
2098
    IMPLICIT NONE
2099

2100
    ! Parameters
2101
    INTEGER, INTENT(IN) :: nvec
2102
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
2103
    ! Result
2104
    INTEGER, INTENT(IN) :: nbasismax
2105
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
2106
    INTEGER, INTENT(INOUT) :: nbasis
2107
    INTEGER :: j
2108

2109
    REAL(KIND=dp), PARAMETER :: c = 1D0/4D0, d = 1D0/SQRT(3D0), &
2110
            e = SQRT(3d0)/6D0
2111
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
2112

2113
    !_ELMER_OMP_SIMD
2114
    DO j=1,nvec
2115
      ! Node 1
2116
      fval(j,nbasis+1) = c*(1-u(j)-d*v(j))*(1-w(j))
2117
      ! Node 2
2118
      fval(j,nbasis+2) = c*(1+u(j)-d*v(j))*(1-w(j))
2119
      ! Node 3
2120
      fval(j,nbasis+3) = e*v(j)*(1-w(j))
2121
      ! Node 4
2122
      fval(j,nbasis+4) = c*(1-u(j)-d*v(j))*(1+w(j))
2123
      ! Node 5
2124
      fval(j,nbasis+5) = c*(1+u(j)-d*v(j))*(1+w(j))
2125
      ! Node 6
2126
      fval(j,nbasis+6) = e*v(j)*(1+w(j))
2127
    END DO
2128
    nbasis = nbasis + 6
2129
  END SUBROUTINE H1Basis_WedgeNodalP
2130

2131
  SUBROUTINE H1Basis_dWedgeNodalP(nvec, u, v, w, nbasismax, grad, nbasis)
2132
    IMPLICIT NONE
2133

2134
    ! Parameters
2135
    INTEGER, INTENT(IN) :: nvec
2136
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
2137
    INTEGER, INTENT(IN) :: nbasismax
2138
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
2139
    INTEGER, INTENT(INOUT) :: nbasis
2140
    
2141
    ! Variables
2142
    INTEGER :: j
2143
    REAL(KIND=dp), PARAMETER :: c = 1D0/4D0, d = 1D0/SQRT(3D0), &
2144
            e = SQRT(3d0)/6D0, f = SQRT(3d0)/12D0
2145
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
2146

2147
    ! First coordinate (xi)
2148
    !_ELMER_OMP_SIMD
2149
    DO j=1,nvec
2150
      grad(j,nbasis+1,1) = -c*(1-w(j))
2151
      grad(j,nbasis+2,1) =  c*(1-w(j))  
2152
      grad(j,nbasis+3,1) =  0
2153
      grad(j,nbasis+4,1) = -c*(1+w(j))
2154
      grad(j,nbasis+5,1) =  c*(1+w(j))
2155
      grad(j,nbasis+6,1) =  0
2156
    END DO
2157
    
2158
    ! Second coordinate (eta)
2159
    !_ELMER_OMP_SIMD
2160
    DO j=1,nvec
2161
      grad(j,nbasis+1,2)=  -f*(1-w(j))
2162
      grad(j,nbasis+2,2) = -f*(1-w(j))
2163
      grad(j,nbasis+3,2) =  e*(1-w(j))
2164
      grad(j,nbasis+4,2) = -f*(1+w(j))
2165
      grad(j,nbasis+5,2) = -f*(1+w(j))
2166
      grad(j,nbasis+6,2) =  e*(1+w(j))
2167
    END DO
2168
    
2169
    ! Third coordinate (zeta)
2170
    !_ELMER_OMP_SIMD
2171
    DO j=1,nvec
2172
      grad(j,nbasis+1,3) = -c*(1-u(j)-d*v(j))
2173
      grad(j,nbasis+2,3) = -c*(1+u(j)-d*v(j))
2174
      grad(j,nbasis+3,3) = -e*v(j)
2175
      grad(j,nbasis+4,3) =  c*(1-u(j)-d*v(j))
2176
      grad(j,nbasis+5,3) =  c*(1+u(j)-d*v(j))
2177
      grad(j,nbasis+6,3) =  e*v(j)
2178
    END DO
2179
        nbasis = nbasis + 6
2180
  END SUBROUTINE H1Basis_dWedgeNodalP
2181

2182
  FUNCTION H1Basis_WedgeL(node, u, v) RESULT(fval)
2183
    IMPLICIT NONE
2184

2185
    ! Parameters
2186
    INTEGER, INTENT(IN) :: node
2187
    REAL(KIND=dp), INTENT(IN) :: u,v
2188
    ! Result
2189
    REAL(KIND=dp) :: fval
2190
    REAL(KIND=dp), PARAMETER :: a=1/SQRT(3.0_dp)
2191
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) &
2192
    !_ELMER_OMP _ELMER_LINEAR_REF(u) _ELMER_LINEAR_REF(v) NOTINBRANCH
2193
    
2194
    SELECT CASE(node)
2195
    CASE (1,4)
2196
      fval = (1-u-a*v)/2
2197
    CASE (2,5)
2198
      fval = (1+u-a*v)/2
2199
    CASE (3,6)
2200
      fval = a*v
2201
    END SELECT
2202
  END FUNCTION H1Basis_WedgeL
2203

2204
  FUNCTION H1Basis_dWedgeL(node) RESULT(grad)
2205
    IMPLICIT NONE
2206

2207
    ! Parameters
2208
    INTEGER, INTENT(IN) :: node
2209
    ! Result
2210
    REAL(KIND=dp) :: grad(3)
2211
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0, d = 1D0/SQRT(3D0)
2212
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) NOTINBRANCH
2213
    
2214
    SELECT CASE(node)
2215
    CASE(1,4)
2216
      grad = c*[REAL(-1,dp), REAL(-d,dp), 0D0]
2217
    CASE(2,5)
2218
      grad = c*[REAL(1,dp), REAL(-d,dp), 0D0]
2219
    CASE(3,6)
2220
      grad =   [REAL(0,dp), REAL(d,dp), 0D0]
2221
    END SELECT    
2222
  END FUNCTION H1Basis_dWedgeL
2223

2224
  FUNCTION H1Basis_WedgeH(node, w) RESULT(fval)
2225
    IMPLICIT NONE
2226

2227
    ! Parameters
2228
    INTEGER, INTENT(IN) :: node
2229
    REAL(KIND=dp), INTENT(IN) :: w
2230
    ! Result
2231
    REAL(KIND=dp) :: fval
2232
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) _ELMER_LINEAR_REF(w) NOTINBRANCH
2233
    
2234
    SELECT CASE(node)
2235
    CASE (1,2,3)
2236
      fval = -w/2
2237
    CASE (4,5,6)
2238
      fval =  w/2
2239
    END SELECT
2240
  END FUNCTION H1Basis_WedgeH
2241

2242
  FUNCTION H1Basis_dWedgeH(node) RESULT(grad)
2243
    IMPLICIT NONE
2244

2245
    ! Parameters
2246
    INTEGER, INTENT(IN) :: node
2247
    ! Result
2248
    REAL(KIND=dp) :: grad(3)
2249
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0
2250
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) NOTINBRANCH
2251
    
2252
    SELECT CASE(node)
2253
    CASE (1,2,3)
2254
      grad = [0D0, 0D0, -c]
2255
    CASE (4,5,6)
2256
      grad = [0D0, 0D0,  c]
2257
    END SELECT
2258
  END FUNCTION H1Basis_dWedgeH
2259
  
2260
  SUBROUTINE H1Basis_WedgeEdgeP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, edgedir)
2261
    IMPLICIT NONE
2262

2263
    INTEGER, INTENT(IN) :: nvec
2264
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
2265
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
2266
    INTEGER, INTENT(IN) :: nbasismax
2267
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
2268
    INTEGER, INTENT(INOUT) :: nbasis
2269
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
2270

2271
    REAL(KIND=dp) :: La, Lb, Na, Nb
2272
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0
2273
    INTEGER :: i,j,k,l, node1, node2
2274
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
2275

2276
    ! For each triangle face edge
2277
    DO i=1,6
2278
      node1 = edgedir(1,i)
2279
      node2 = edgedir(2,i)
2280
      DO j=1,pmax(i)-1
2281
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb)
2282
        DO k=1,nvec
2283
          La = H1Basis_WedgeL(node1, u(k), v(k))
2284
          Lb = H1Basis_WedgeL(node2, u(k), v(k))
2285
          Na = fval(k,node1)
2286
          Nb = fval(k,node2)
2287
          fval(k, nbasis+j) = Na*Nb*H1Basis_varPhi(j+1,Lb-La)
2288
        END DO
2289
      END DO
2290
      ! nbasis = nbasis + (pmax(i)-2) + 1
2291
      nbasis = nbasis + pmax(i) - 1
2292
    END DO
2293
      
2294
    ! For each square face edge
2295
    DO i=7,9
2296
      node1 = edgedir(1,i)
2297
      node2 = edgedir(2,i)
2298
      DO j=1,pmax(i)-1
2299
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb)
2300
        DO k=1,nvec
2301
          La = H1Basis_WedgeH(node1, w(k))
2302
          Lb = H1Basis_WedgeH(node2, w(k))
2303
          Na = fval(k,node1)
2304
          Nb = fval(k,node2)
2305
          fval(k, nbasis+j) = Na*Nb*H1Basis_varPhi(j+1, Lb-La)
2306
        END DO
2307
      END DO
2308
      ! nbasis = nbasis + (pmax(i)-2) + 1
2309
      nbasis = nbasis + pmax(i) - 1
2310
    END DO
2311
  END SUBROUTINE H1Basis_WedgeEdgeP
2312
  
2313
  SUBROUTINE H1Basis_dWedgeEdgeP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, edgedir)
2314
    IMPLICIT NONE
2315

2316
    INTEGER, INTENT(IN) :: nvec
2317
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
2318
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
2319
    INTEGER, INTENT(IN) :: nbasismax
2320
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
2321
    INTEGER, INTENT(INOUT) :: nbasis
2322
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
2323

2324
    REAL(KIND=dp) :: La, Lb, Na, Nb, Phi, dPhi, vPhi, dVPhi, &
2325
          dLa(3), dLb(3), dNa(3), dNb(3), N(VECTOR_BLOCK_LENGTH,nbasismax)
2326
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0
2327
    INTEGER :: i,j,k, node1, node2, nnb
2328
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
2329

2330
    nnb = 0; CALL H1Basis_WedgeNodalP(nvec, u, v, w, nbasismax, n, nnb )
2331

2332
    ! For each triangle face edge
2333
    DO i=1,6
2334
      node1 = edgedir(1,i)
2335
      node2 = edgedir(2,i)
2336
      dLa = H1Basis_dWedgeL(node1)
2337
      dLb = H1Basis_dWedgeL(node2)
2338

2339
      DO j=1,pmax(i)-1
2340
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb, dNa, dNb, vPhi, dVPhi)
2341
        DO k=1,nvec
2342
          La = H1Basis_WedgeL(node1, u(k), v(k))
2343
          Lb = H1Basis_WedgeL(node2, u(k), v(k))
2344

2345
          Na = N(k,node1)
2346
          Nb = N(k,node2)
2347
          dNa = grad(k,node1,:)
2348
          dNb = grad(k,node2,:)
2349

2350
          vPhi = H1Basis_varPhi(j+1,Lb-La)
2351
          dVPhi = H1Basis_dVarPhi(j+1,Lb-La)
2352
          ! fval(k, nbasis+j-1) = c*La*Lb*H1Basis_varPhi(j,Lb-La)*(1+Na+Nb)
2353

2354
          grad(k, nbasis+j,:) = dNa*Nb*vPhi + Na*dNb*vPhi + Na*Nb*dvPhi*(dLb-dLa)
2355
        END DO
2356
      END DO
2357
      ! nbasis = nbasis + (pmax(i)-2) + 1
2358
      nbasis = nbasis + pmax(i) - 1
2359
    END DO
2360
      
2361
    ! For each square face edge
2362
    DO i=7,9
2363
      node1 = edgedir(1,i)
2364
      node2 = edgedir(2,i)
2365
      dLa = H1Basis_dWedgeH(node1)
2366
      dLb = H1Basis_dWedgeH(node2)
2367
      DO j=1,pmax(i)-1
2368
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb, dNa, dNb, vPhi, dvPhi)
2369
        DO k=1,nvec
2370
          La = H1Basis_WedgeH(node1, w(k))
2371
          Lb = H1Basis_WedgeH(node2, w(k))
2372

2373
          vPhi = H1Basis_varPhi(j+1, Lb-La)
2374
          dvPhi = H1Basis_dvarPhi(j+1, Lb-La)
2375

2376
          Na = N(k,node1)
2377
          Nb = N(k,node2)
2378
          dNa = grad(k,node1,:)
2379
          dNb = grad(k,node2,:)
2380

2381
          grad(k, nbasis+j,:) = dNa*Nb*vPhi + Na*dNb*vPhi + Na*Nb*dvPhi*(dLb-dLa)
2382
        END DO
2383
      END DO
2384
      ! nbasis = nbasis + (pmax(i)-2) + 1
2385
      nbasis = nbasis + pmax(i) - 1
2386
    END DO
2387
  END SUBROUTINE H1Basis_dWedgeEdgeP
2388
  
2389
  SUBROUTINE H1Basis_WedgeFaceP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, facedir)
2390
    IMPLICIT NONE
2391

2392
    INTEGER, INTENT(IN) :: nvec
2393
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
2394
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
2395
    INTEGER, INTENT(IN) :: nbasismax
2396
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
2397
    INTEGER, INTENT(INOUT) :: nbasis
2398
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
2399

2400
    REAL(KIND=dp) :: La, Lb, Lc, Na, Nb, Lha, Lhc
2401
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0
2402
    INTEGER :: i,j,k,l, node1, node2, node3, node4
2403
    LOGICAL :: nonpermuted
2404
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
2405

2406
    ! Triangle faces
2407
    DO i=1,2      
2408
      DO j=0,pmax(i)-3
2409
        DO k=1,pmax(i)-j-2
2410
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Na)
2411
          DO l=1,nvec
2412
            La = H1Basis_WedgeL(facedir(1,i), u(l), v(l))
2413
            Lb = H1Basis_WedgeL(facedir(2,i), u(l), v(l))
2414
            Lc = H1Basis_WedgeL(facedir(3,i), u(l), v(l))
2415
            Na = H1Basis_WedgeH(facedir(1,i), w(l))
2416

2417
            fval(l,nbasis+k) = c*(1+2*Na)*La*Lb*Lc* &
2418
                                            H1Basis_LegendreP(j, Lb-La)* &
2419
                                            H1Basis_LegendreP(k-1, 2*Lc-1)
2420
          END DO
2421
        END DO
2422
        ! nbasis = nbasis + (p-j-3) + 1
2423
        nbasis = nbasis + MAX(pmax(i)-j-2,0)
2424
      END DO
2425
    END DO
2426
    ! Square faces
2427
    DO i=3,5
2428
      ! First and second node must form an edge in upper or lower triangle
2429
      node1 = facedir(1,i)
2430
      node2 = facedir(2,i)
2431
      node3 = facedir(3,i)
2432
      node4 = facedir(4,i)
2433

2434
      nonpermuted = (node1 >= 1 .AND. node1 <= 3 .AND.&
2435
                     node2 >= 1 .AND. node2 <= 3) .OR. &
2436
                    (node1 >= 4 .AND. node1 <= 6 .AND.&
2437
                     node2 >= 4 .AND. node2 <= 6)
2438

2439
      
2440
      DO j=0,pmax(i)-2
2441
        DO k=1,pmax(i)-1
2442
          IF (nonpermuted) THEN
2443
            !_ELMER_OMP_SIMD PRIVATE(La,Lb,Na,Nb,Lha,Lhc)
2444
            DO l=1,nvec
2445
              La = H1Basis_WedgeL(node1, u(l), v(l))
2446
              Lb = H1Basis_WedgeL(node2, u(l), v(l))
2447

2448
              Lha = H1Basis_WedgeH(node1, w(l))
2449
              Lhc = H1Basis_WedgeH(node4, w(l))
2450

2451
              Na = fval(l,node1)
2452
              Nb = fval(l,node3)
2453

2454
              fval(l,nbasis+k) = Na*Nb*H1Basis_LegendreP(j, Lb-La)*H1Basis_LegendreP(k-1,Lhc-Lha)
2455
            END DO
2456
          ELSE
2457
            !_ELMER_OMP_SIMD PRIVATE(La,Lb,Na,Nb,Lha,Lhc)
2458
            DO l=1,nvec
2459
              ! Numbering permuted, switch indices j,k and nodes 2 and 4
2460
              La = H1Basis_WedgeL(node1, u(l), v(l))
2461
              Lb = H1Basis_WedgeL(node4, u(l), v(l))
2462

2463
              Lha = H1Basis_WedgeH(node1, w(l))
2464
              Lhc = H1Basis_WedgeH(node2, w(l))
2465

2466
              Na = fval(l,node1)
2467
              Nb = fval(l,node3)
2468

2469
              fval(l,nbasis+k) = Na*Nb*H1Basis_LegendreP(k-1, Lb-La)*H1Basis_LegendreP(j,Lhc-Lha)
2470
            END DO
2471
          END IF
2472
        END DO
2473
        nbasis = nbasis + MAX(pmax(i)-1,0)
2474
      END DO
2475
    END DO
2476
  END SUBROUTINE H1Basis_WedgeFaceP
2477
  
2478
  SUBROUTINE H1Basis_dWedgeFaceP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, facedir)
2479
    IMPLICIT NONE
2480

2481
    INTEGER, INTENT(IN) :: nvec
2482
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
2483
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
2484
    INTEGER, INTENT(IN) :: nbasismax
2485
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
2486
    INTEGER, INTENT(INOUT) :: nbasis
2487
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
2488

2489
    REAL(KIND=dp) :: La, Lb, Lc, Lha, Lhc, Na, Nb, Nc, LegPLbLa, LegP2Lc1, &
2490
            cNa, vPhi, Phi, dLa(3), dLb(3), dLc(3), dLhc(3), dLha(3), dNc(3), &
2491
                 Legj, Legk, dLegj(3), dLegk(3), dNa(3), dNb(3), N(VECTOR_BLOCK_LENGTH, nbasismax) 
2492
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0
2493
    INTEGER :: i,j,k,l, node1, node2, node3, node4, nnb
2494
    LOGICAL :: nonpermuted
2495
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
2496

2497
    nnb = 0; CALL H1Basis_WedgeNodalP(nvec, u, v, w, nbasismax, n, nnb)
2498

2499
    ! Triangle faces
2500
    DO i=1,2
2501
      dLa = H1Basis_dWedgeL(facedir(1,i))
2502
      dLb = H1Basis_dWedgeL(facedir(2,i))
2503
      dLc = H1Basis_dWedgeL(facedir(3,i))
2504
      dNa = H1Basis_dWedgeH(facedir(1,i))
2505
      DO j=0,pmax(i)-3
2506
        DO k=1,pmax(i)-j-2
2507
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Na, LegPLbLa, LegP2Lc1, cNa)
2508
          DO l=1,nvec
2509
            La = H1Basis_WedgeL(facedir(1,i), u(l), v(l))
2510
            Lb = H1Basis_WedgeL(facedir(2,i), u(l), v(l))
2511
            Lc = H1Basis_WedgeL(facedir(3,i), u(l), v(l))
2512
            Na = H1Basis_WedgeH(facedir(1,i), w(l))
2513

2514
            LegPLbLa = H1Basis_LegendreP(j, Lb-La)
2515
            LegP2Lc1 = H1Basis_LegendreP(k-1, 2*Lc-1)
2516
            cNa = c*(1+2*Na)
2517

2518
            grad(l,nbasis+k,1) = cNa*dLa(1)*Lb*Lc*LegPLbLa*LegP2Lc1 + &
2519
                    cNa*La*dLb(1)*Lc*LegPLbLa*LegP2Lc1 + &
2520
                    cNa*La*Lb*dLc(1)*LegPLbLa*LegP2Lc1 + &
2521
                    cNa*La*Lb*Lc*H1Basis_dLegendreP(j, Lb-La)*(dLb(1)-dLa(1))*LegP2Lc1 + &
2522
                    cNa*La*Lb*Lc*LegPLbLa*H1Basis_dLegendreP(k-1, 2*Lc-1)*(2*dLc(1))
2523
            ! + 2*c*dNa(1)*La*Lb*Lc*LegPLbLa*LegP2Lc1 ! Term always zero
2524
            grad(l,nbasis+k,2) = cNa*dLa(2)*Lb*Lc*LegPLbLa*LegP2Lc1 + &
2525
                    cNa*La*dLb(2)*Lc*LegPLbLa*LegP2Lc1 + &
2526
                    cNa*La*Lb*dLc(2)*LegPLbLa*LegP2Lc1 + &
2527
                    cNa*La*Lb*Lc*H1Basis_dLegendreP(j, Lb-La)*(dLb(2)-dLa(2))*LegP2Lc1 + &
2528
                    cNa*La*Lb*Lc*LegPLbLa*H1Basis_dLegendreP(k-1, 2*Lc-1)*(2*dLc(2))
2529
            ! + 2*c*dNa(2)*La*Lb*Lc*LegPLbLa*LegP2Lc1 ! Term always zero
2530
            grad(l,nbasis+k,3) = 2*c*dNa(3)*La*Lb*Lc*LegPLbLa*LegP2Lc1
2531
            ! Rest of the terms always zero
2532
            ! cNa*dLa(3)*Lb*Lc*LegPLbLa*LegP2Lc1 + &
2533
            ! cNa*La*dLb(3)*Lc*LegPLbLa*LegP2Lc1 + &
2534
            ! cNa*La*Lb*dLc(3)*LegPLbLa*LegP2Lc1 + &
2535
            ! cNa*La*Lb*Lc*H1Basis_dLegendreP(j, Lb-La)*(dLb(3)-dLa(3))*LegP2Lc1 + &
2536
            ! cNa*La*Lb*Lc*LegPLbLa*H1Basis_dLegendreP(k, 2*Lc-1)*(2*dLc(3))
2537
          END DO
2538
        END DO
2539
        ! nbasis = nbasis + (p-j-3) + 1
2540
        nbasis = nbasis + MAX(pmax(i)-j-2,0)
2541
      END DO
2542
    END DO
2543

2544
    ! Square faces
2545
    DO i=3,5
2546

2547
      node1 = facedir(1,i)
2548
      node2 = facedir(2,i)
2549
      node3 = facedir(3,i)
2550
      node4 = facedir(4,i)
2551

2552
      ! First and second node must form an edge in upper or lower triangle
2553
      nonpermuted = (node1 >= 1 .AND. node1 <= 3 .AND.&
2554
                     node2 >= 1 .AND. node2 <= 3) .OR. &
2555
                    (node1 >= 4 .AND. node1 <= 6 .AND.&
2556
                     node2 >= 4 .AND. node2 <= 6)
2557
      
2558
      IF (nonpermuted) THEN
2559
        dLa  = H1Basis_dWedgeL(node1)
2560
        dLb  = H1Basis_dWedgeL(node2)
2561
        dLha = H1Basis_dWedgeH(node1)
2562
        dLhc = H1Basis_dWedgeH(node4)
2563
      ELSE
2564
        dLa  = H1Basis_dWedgeL(node1)
2565
        dLb  = H1Basis_dWedgeL(node4)
2566
        dLha = H1Basis_dWedgeH(node1)
2567
        dLhc = H1Basis_dWedgeH(node2)
2568
      END IF
2569

2570
      DO j=0,pmax(i)-2
2571
        DO k=1,pmax(i)-1
2572
          IF (nonpermuted) THEN
2573
            !_ELMER_OMP_SIMD PRIVATE(La,Lb,Lha,Lhc,Na,Nb,dNa,dNb,Legj,Legk,dLegj,dLegk)
2574
            DO l=1,nvec
2575
              La = H1Basis_WedgeL(node1, u(l), v(l))
2576
              Lb = H1Basis_WedgeL(node2, u(l), v(l))
2577

2578
              Lha = H1Basis_WedgeH(node1, w(l))
2579
              Lhc = H1Basis_WedgeH(node4, w(l))
2580

2581
              Na = N(l,node1)
2582
              Nb = N(l,node3)
2583

2584
              dNa = grad(l,node1,:)
2585
              dNb = grad(l,node3,:)
2586

2587
              Legj = H1Basis_LegendreP(j,Lb-La)
2588
              Legk = H1Basis_LegendreP(k-1,Lhc-Lha)
2589

2590
              dLegj = H1Basis_dLegendreP(j,Lb-La)*(dLb-dLa)
2591
              dLegk = H1Basis_dLegendreP(k-1,Lhc-Lha)*(dLhc-dLha)
2592

2593
              grad(l,nbasis+k,:) = dNa*Nb*Legj*Legk + Na*dNb*Legj*Legk + &
2594
                     Na*Nb*dLegj*Legk + Na*Nb*Legj*dLegk
2595
            END DO
2596
          ELSE
2597
            !_ELMER_OMP_SIMD PRIVATE(La,Lb,Lha,Lhc,Na,Nb,dNa,dNb,Legj,Legk,dLegj,dLegk)
2598
            DO l=1,nvec
2599
              La = H1Basis_WedgeL(node1, u(l), v(l))
2600
              Lb = H1Basis_WedgeL(node4, u(l), v(l))
2601

2602
              Lha = H1Basis_WedgeH(node1, w(l))
2603
              Lhc = H1Basis_WedgeH(node2, w(l))
2604

2605
              Na = N(l,node1)
2606
              Nb = N(l,node3)
2607

2608
              dNa = grad(l,node1,:)
2609
              dNb = grad(l,node3,:)
2610

2611
              Legj = H1Basis_LegendreP(k-1,Lb-La)
2612
              Legk = H1Basis_LegendreP(j,Lhc-Lha)
2613

2614
              dLegj = H1Basis_dLegendreP(k-1,Lb-La)*(dLb-dLa)
2615
              dLegk = H1Basis_dLegendreP(j,Lhc-Lha)*(dLhc-dLha)
2616

2617
              grad(l,nbasis+k,:) = dNa*Nb*Legj*Legk + Na*dNb*Legj*Legk + &
2618
                     Na*Nb*dLegj*Legk + Na*Nb*Legj*dLegk
2619
            END DO
2620
          END IF
2621
        END DO
2622
        nbasis = nbasis + MAX(pmax(i)-1,0)
2623
      END DO
2624
    END DO
2625
  END SUBROUTINE H1Basis_dWedgeFaceP
2626
  
2627

2628
  SUBROUTINE H1Basis_WedgeBubbleP(nvec, u, v, w, pmax, nbasismax, fval, nbasis)
2629
    IMPLICIT NONE
2630

2631
    INTEGER, INTENT(IN) :: nvec
2632
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
2633
    INTEGER, INTENT(IN) :: pmax
2634
    INTEGER, INTENT(IN) :: nbasismax
2635
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
2636
    INTEGER, INTENT(INOUT) :: nbasis
2637
    
2638
    INTEGER :: i, j, k, l
2639
    REAL(KIND=dp) :: L1, L2, L3, s,t
2640
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
2641

2642
    DO i=0,pmax-3
2643
      DO j=0,pmax-3-i
2644
        DO k=1,pmax-1
2645
          !_ELMER_OMP_SIMD PRIVATE(L1, L2, L3, s,t)
2646
          DO l=1,nvec
2647
            L1 = H1Basis_WedgeL(1,u(l),v(l))
2648
            L2 = H1Basis_WedgeL(2,u(l),v(l))
2649
            L3 = H1Basis_WedgeL(3,u(l),v(l))
2650

2651
            s = L2-L1
2652
            t = 2*L3-1
2653

2654
            ! Get value of bubble function
2655
            fval(l,nbasis+k) = L1*L2*L3*H1Basis_LegendreP(i,s)*&
2656
                    H1Basis_LegendreP(j,t)*H1Basis_Phi(k+1,w(l))
2657
          END DO
2658
        END DO
2659
        ! nbasis = nbasis + (pmax-3-i-j)-2+1
2660
        nbasis = nbasis + MAX(pmax-1,0)
2661
      END DO
2662
    END DO
2663
  END SUBROUTINE H1Basis_WedgeBubbleP
2664

2665
  SUBROUTINE H1Basis_dWedgeBubbleP(nvec, u, v, w, pmax, nbasismax, grad, nbasis)
2666
    IMPLICIT NONE
2667

2668
    INTEGER, INTENT(IN) :: nvec
2669
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
2670
    INTEGER, INTENT(IN) :: pmax
2671
    INTEGER, INTENT(IN) :: nbasismax
2672
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
2673
    INTEGER, INTENT(INOUT) :: nbasis
2674
    
2675
    ! Parameters
2676
    INTEGER :: i,j,k,l
2677
    ! Variables
2678
    REAL(KIND=dp) :: L1,L2,L3,Legi,Legj, Legk, dLegi(3), dLegj(3), dLegk(3), &
2679
                  s,t,ds(3),dt(3), dL1(3),dL2(3),dL3(3)
2680
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64      
2681

2682
    dL1 = H1Basis_dWedgeL(1)
2683
    dL2 = H1Basis_dWedgeL(2)
2684
    dL3 = H1Basis_dWedgeL(3)
2685

2686
    DO i=0,pmax-3
2687
      DO j=0,pmax-3-i
2688
        DO k=1,pmax-1
2689
          !_ELMER_OMP_SIMD PRIVATE(L1, L2, L3, Legi, Legj, Legk, dLegi, dLegj, dLegk, s,t,ds,dt)
2690
          DO l=1,nvec
2691

2692
            ! Values of function L
2693
            L1 = H1Basis_WedgeL(1,u(l),v(l))
2694
            L2 = H1Basis_WedgeL(2,u(l),v(l))
2695
            L3 = H1Basis_WedgeL(3,u(l),v(l))
2696
            
2697
            s = L2-L1
2698
            t = 2*L3-1
2699
            ds = dL2-dL1
2700
            dt = 2*dL3
2701

2702
            Legi = H1Basis_LegendreP(i,s)
2703
            Legj = H1Basis_LegendreP(j,t)
2704
            Legk = H1Basis_Phi(k+1,w(l))
2705

2706
            dLegi = H1Basis_dLegendreP(i,s)*ds
2707
            dLegj = H1Basis_dLegendreP(j,t)*dt
2708
            dLegk = H1Basis_dPhi(k+1,w(l))*[0,0,1]
2709

2710
            grad(l,nbasis+k,:) = dL1*L2*L3*Legi*Legj*Legk + L1*dL2*L3*Legi*Legj*Legk + &
2711
                    L1*L2*dL3*Legi*Legj*Legk + L1*L2*L3*dLegi*Legj*Legk + &
2712
                    L1*L2*L3*Legi*dLegj*Legk + L1*L2*L3*Legi*Legj*dLegk
2713
          END DO
2714
        END DO
2715
        ! nbasis = nbasis + (pmax-3-i-j)-2+1
2716
        nbasis = nbasis + MAX(pmax-1,0)
2717
      END DO
2718
    END DO
2719
  END SUBROUTINE H1Basis_dWedgeBubbleP
2720

2721
!------------------------------------------------------------------------------
2722
!>     Pyramid nodal basis at points (u,v,w)
2723
!------------------------------------------------------------------------------
2724
  SUBROUTINE H1Basis_PyramidNodalP(nvec, u, v, w, nbasismax, fval, nbasis)
2725
    IMPLICIT NONE
2726

2727
    ! Parameters
2728
    INTEGER, INTENT(IN) :: nvec
2729
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
2730
    ! Result
2731
    INTEGER, INTENT(IN) :: nbasismax
2732
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
2733
    INTEGER, INTENT(INOUT) :: nbasis
2734
    INTEGER :: j
2735
!------------------------------------------------------------------------------
2736
    REAL(KIND=dp) :: s, sq2 = SQRT(2.0_dp)
2737
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
2738

2739
      !_ELMER_OMP_SIMD
2740
      DO j=1,nvec
2741
        s = w(j) / sq2
2742
        fval(j,nbasis+1) = (1-u(j)-v(j)-s+u(j)*v(j)/(1-s))/4
2743
        fval(j,nbasis+2) = (1+u(j)-v(j)-s-u(j)*v(j)/(1-s))/4
2744
        fval(j,nbasis+3) = (1+u(j)+v(j)-s+u(j)*v(j)/(1-s))/4
2745
        fval(j,nbasis+4) = (1-u(j)+v(j)-s-u(j)*v(j)/(1-s))/4
2746
        fval(j,nbasis+5) =  s
2747
      END DO
2748
      nbasis = nbasis + 5
2749
    END SUBROUTINE H1Basis_PyramidNodalP
2750

2751

2752
!------------------------------------------------------------------------------
2753
!>     Gradient of pyramids nodal basis at point (u,v,w)
2754
!------------------------------------------------------------------------------
2755
  SUBROUTINE H1Basis_dPyramidNodalP(nvec, u, v, w, nbasismax, grad, nbasis)
2756
    IMPLICIT NONE
2757

2758
    ! Parameters
2759
    INTEGER, INTENT(IN) :: nvec
2760
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
2761
    ! Result
2762
    INTEGER, INTENT(IN) :: nbasismax
2763
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
2764
    INTEGER, INTENT(INOUT) :: nbasis
2765
    INTEGER :: j
2766
    REAL(KIND=dp) :: s, sq2=SQRT(2.0_dp), sq42=4*SQRT(2._dp)
2767
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64
2768

2769
      !_ELMER_OMP_SIMD
2770
      DO j=1,nvec
2771
        s = v(j)/(1-w(j)/sq2)
2772
        grad(j,nbasis+1,1) = (-1 + s)/4
2773
        grad(j,nbasis+2,1) = ( 1 - s)/4
2774
        grad(j,nbasis+3,1) = ( 1 + s)/4
2775
        grad(j,nbasis+4,1) = (-1 - s)/4
2776
        grad(j,nbasis+5,1) = 0
2777
      END DO
2778

2779
      !_ELMER_OMP_SIMD
2780
      DO j=1,nvec
2781
        s = u(j)/(1-w(j)/sq2)
2782
        grad(j,nbasis+1,2) = (-1 + s)/4
2783
        grad(j,nbasis+2,2) = (-1 - s)/4
2784
        grad(j,nbasis+3,2) = ( 1 + s)/4
2785
        grad(j,nbasis+4,2) = ( 1 - s)/4
2786
        grad(j,nbasis+5,2) = 0
2787
      END DO
2788

2789
      !_ELMER_OMP_SIMD
2790
      DO j=1,nvec
2791
        s = u(j)*v(j)/(1-w(j)/sq2)**2
2792
        grad(j,nbasis+1,3) = (-1 + s)/sq42
2793
        grad(j,nbasis+2,3) = (-1 - s)/sq42
2794
        grad(j,nbasis+3,3) = (-1 + s)/sq42
2795
        grad(j,nbasis+4,3) = (-1 - s)/sq42
2796
        grad(j,nbasis+5,3) = 1/sq2
2797
      END DO
2798
      nbasis = nbasis + 5
2799
  END SUBROUTINE H1Basis_dPyramidnodalP
2800

2801
    ! Define affine coordinates for pyramid square face
2802
    FUNCTION H1Basis_PyramidL(which, u, v) RESULT(value)
2803
      IMPLICIT NONE
2804

2805
      ! Parameters
2806
      INTEGER, INTENT(IN) :: which
2807
      REAL(KIND=dp), INTENT(IN) :: u,v
2808
      ! Variables
2809
      REAL(KIND=dp) :: value
2810
      
2811
      SELECT CASE (which)
2812
      CASE (1)
2813
         value = ((1-u)+(1-v))/2 
2814
      CASE (2)
2815
         value = ((1+u)+(1-v))/2
2816
      CASE (3)
2817
         value = ((1+u)+(1+v))/2
2818
      CASE (4)
2819
         value = ((1-u)+(1+v))/2
2820
      CASE DEFAULT
2821
         CALL Fatal('PElementBase::PyramidL','Unknown affine coordinate for square face')
2822
      END SELECT
2823
    END FUNCTION H1Basis_PyramidL
2824

2825
    FUNCTION H1Basis_dPyramidL(which) RESULT(grad)
2826
      IMPLICIT NONE
2827
      
2828
      ! Parameters
2829
      INTEGER, INTENT(IN) :: which
2830
      ! Variables
2831
      REAL(KIND=dp) :: grad(3)
2832

2833
      SELECT CASE (which)
2834
      CASE (1)
2835
         grad = [-1d0/2,-1d0/2, 0d0 ]
2836
      CASE (2)
2837
         grad = [ 1d0/2,-1d0/2, 0d0 ]
2838
      CASE (3)
2839
         grad = [ 1d0/2, 1d0/2, 0d0 ]
2840
      CASE (4)
2841
         grad = [-1d0/2, 1d0/2, 0d0 ]
2842
      CASE DEFAULT
2843
         CALL Fatal('PElementBase::dPyramidL','Unknown affine coordinate for square face')
2844
      END SELECT
2845
    END FUNCTION H1Basis_dPyramidL
2846

2847

2848
    FUNCTION H1Basis_PyramidTL(which, u, v, w) RESULT(value)
2849
      IMPLICIT NONE
2850
      
2851
      ! Parameters
2852
      INTEGER, INTENT(IN) :: which
2853
      REAL(KIND=dp), INTENT(IN) :: u,v,w 
2854
      ! Variables
2855
      REAL(KIND=dp) :: value,s
2856
      
2857
      value = 0
2858
      s = w/SQRT(2.0_dp)
2859
      SELECT CASE(which)
2860
      CASE (1)
2861
         value = (2-u-v-s)/2
2862
      CASE (2)
2863
         value = (2+u-v-s)/2
2864
      CASE (3)
2865
         value = (2+u+v-s)/2
2866
      CASE (4)
2867
         value = (2-u+v-s)/2
2868
      CASE (5)
2869
         value = s
2870
      CASE DEFAULT
2871
         CALL Fatal('PElementBase::PyramidTL','Unknown function L for brick')
2872
      END SELECT
2873
    END FUNCTION H1Basis_PyramidTL
2874

2875
    FUNCTION H1Basis_dPyramidTL(which) RESULT(grad)
2876
      IMPLICIT NONE
2877
      
2878
      ! Parameters
2879
      INTEGER, INTENT(IN) :: which
2880
      ! Variables
2881
      REAL(KIND=dp) :: grad(3),s
2882
      
2883
      grad = 0
2884
      SELECT CASE(which)
2885
      CASE (1)
2886
         grad(1) = -1
2887
         grad(2) = -1
2888
         grad(3) = -1
2889
      CASE (2)
2890
         grad(1) =  1
2891
         grad(2) = -1
2892
         grad(3) = -1
2893
      CASE (3)
2894
         grad(1) =  1
2895
         grad(2) =  1
2896
         grad(3) = -1
2897
      CASE (4)
2898
         grad(1) = -1
2899
         grad(2) =  1
2900
         grad(3) = -1
2901
      CASE (5)
2902
         grad(3) =  2
2903
      CASE DEFAULT
2904
         CALL Fatal('PElementBase::PyramidTL','Unknown function L for brick')
2905
      END SELECT
2906
      grad = grad/2
2907
      grad(3) = grad(3)/SQRT(2._dp)
2908
    END FUNCTION H1Basis_dPyramidTL
2909

2910

2911

2912
  SUBROUTINE H1Basis_PyramidEdgeP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, edgedir)
2913
    IMPLICIT NONE
2914

2915
    INTEGER, INTENT(IN) :: nvec
2916
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
2917
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
2918
    INTEGER, INTENT(IN) :: nbasismax
2919
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
2920
    INTEGER, INTENT(INOUT) :: nbasis
2921
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
2922

2923
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax) :: n
2924
    REAL(KIND=dp) :: La, Lb, Na, Nb
2925
    REAL(KIND=dp), PARAMETER :: c = 1D0/2D0
2926
    INTEGER :: i,j,k,l, node1, node2, nnb
2927
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
2928

2929
    DO i=1,4 ! square face
2930
      node1 = edgedir(1,i)
2931
      node2 = edgedir(2,i)
2932

2933
      DO j=1,pmax(i)-1
2934
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb)
2935
        DO k=1,nvec
2936
          Na = fval(k,node1)
2937
          Nb = fval(k,node2)
2938
          La = H1Basis_PyramidL(node1,u(k),v(k))
2939
          Lb = H1Basis_PyramidL(node2,u(k),v(k))
2940

2941
          fval(k,nbasis+j) = Na*Nb*H1Basis_varPhi(j+1,Lb-La)
2942
        END DO
2943
      END DO
2944
      nbasis = nbasis + MAX(pmax(i)-1,0)
2945
    END DO
2946

2947
    DO i=5,8 ! triangular faces
2948
      Node1 = edgedir(1,i)
2949
      Node2 = edgedir(2,i)
2950

2951
      DO j=1,pmax(i)-1
2952
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb)
2953
        DO k=1,nvec
2954
          Na = fval(k,node1)
2955
          Nb = fval(k,node2)
2956
          La = H1Basis_PyramidTL(node1,u(k),v(k),w(k))
2957
          Lb = H1Basis_PyramidTL(node2,u(k),v(k),w(k))
2958
          fval(k,nbasis+j) = Na*Nb*H1Basis_varPhi(j+1,Lb-La)
2959
        END DO
2960
      END DO
2961
      nbasis = nbasis + MAX(pmax(i)-1,0)
2962
    END DO
2963
  END SUBROUTINE H1Basis_PyramidEdgeP
2964

2965

2966
  SUBROUTINE H1Basis_dPyramidEdgeP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, edgedir)
2967
    IMPLICIT NONE
2968

2969
    INTEGER, INTENT(IN) :: nvec
2970
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
2971
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
2972
    INTEGER, INTENT(IN) :: nbasismax
2973
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
2974
    INTEGER, INTENT(INOUT) :: nbasis
2975
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
2976

2977
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax) :: n
2978
    REAL(KIND=dp) :: La, Lb, dLa(3), dLb(3), Na, Nb, dNa(3), dNb(3), Phi, dPhi(3)
2979
    INTEGER :: i,j,k,l, node1, node2, nnb
2980
!!!!DIR$ ASSUME_ALIGNED u:64, v:64, w:64
2981

2982
    REAL(KIND=dp) :: s, sq2=SQRT(2.0_dp)
2983

2984
    nnb=0; CALL H1Basis_PyramidNodalP(nvec, u, v, w, nbasismax, n, nnb)
2985

2986
    DO i=1,4 ! square face
2987
      node1 = edgedir(1,i)
2988
      node2 = edgedir(2,i)
2989

2990
      dLa = H1Basis_dPyramidL(node1)
2991
      dLb = H1Basis_dPyramidL(node2)
2992

2993
      DO j=1,pmax(i)-1
2994
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb, dNa, dNb, Phi, dPhi)
2995
        DO k=1,nvec
2996
          Na = N(k,node1)
2997
          Nb = N(k,node2)
2998

2999
          dNa = grad(k,node1,:)
3000
          dNb = grad(k,node2,:)
3001

3002
          La = H1Basis_PyramidL(node1,u(k),v(k))
3003
          Lb = H1Basis_PyramidL(node2,u(k),v(k))
3004

3005
          Phi = H1Basis_varPhi(j+1,Lb-La)
3006
          dPhi = H1Basis_dvarPhi(j+1,Lb-La)*(dLb-dLa)
3007

3008
          grad(k,nbasis+j,:) = dNa*Nb*Phi + Na*dNb*Phi + Na*Nb*dPhi
3009
        END DO
3010
      END DO
3011
      nbasis = nbasis + MAX(pmax(i)-1,0)
3012
    END DO
3013

3014
    DO i=5,8 ! triangular faces
3015
      Node1 = edgedir(1,i)
3016
      Node2 = edgedir(2,i)
3017

3018
      dLa = H1Basis_dPyramidTL(node1)
3019
      dLb = H1Basis_dPyramidTL(node2)
3020

3021
      DO j=1,pmax(i)-1
3022
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb, dNa, dNb, Phi, dPhi)
3023
        DO k=1,nvec
3024
          Na = N(k,node1)
3025
          Nb = N(k,node2)
3026
          dNa = grad(k,node1,:)
3027
          dNb = grad(k,node2,:)
3028

3029
          La = H1Basis_PyramidTL(node1,u(k),v(k),w(k))
3030
          Lb = H1Basis_PyramidTL(node2,u(k),v(k),w(k))
3031

3032
          Phi = H1Basis_varPhi(j+1,Lb-La)
3033
          dPhi = H1Basis_dvarPhi(j+1,Lb-La)*(dLb-dLa)
3034

3035
          grad(k,nbasis+j,:) = dNa*Nb*Phi + Na*dNb*Phi + Na*Nb*dPhi
3036
        END DO
3037
      END DO
3038
      nbasis = nbasis + MAX(pmax(i)-1,0)
3039
    END DO
3040
  END SUBROUTINE H1Basis_dPyramidEdgeP
3041

3042

3043

3044
  SUBROUTINE H1Basis_PyramidFaceP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, facedir)
3045
    IMPLICIT NONE
3046

3047
    INTEGER, INTENT(IN) :: nvec
3048
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3049
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
3050
    INTEGER, INTENT(IN) :: nbasismax
3051
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
3052
    INTEGER, INTENT(INOUT) :: nbasis
3053
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
3054

3055
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax) :: n
3056
    REAL(KIND=dp) :: La, Lb, Lc, Pa, Pb, Pc
3057
    INTEGER :: i,j,k,l, node1, node2, node3, node4, nnb
3058
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
3059

3060
    DO i=1,1 ! square face
3061
      node1 = facedir(1,i)
3062
      node2 = facedir(2,i)
3063
      node3 = facedir(3,i)
3064
      node4 = facedir(4,i)
3065

3066
      DO j=0,pmax(i)-2
3067
        DO k=1,pmax(i)-1
3068
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Pa, Pb)
3069
          DO l=1,nvec
3070
            Pa = fval(l,node1)
3071
            Pb = fval(l,node3)
3072

3073
            La = H1Basis_PyramidL(node1,u(l),v(l))
3074
            Lb = H1Basis_PyramidL(node2,u(l),v(l))
3075
            Lc = H1Basis_PyramidL(node4,u(l),v(l))
3076
            fval(l,nbasis+k) = Pa*Pb*H1Basis_LegendreP(j,Lb-La)*H1Basis_LegendreP(k-1,Lc-La)
3077
          END DO
3078
        END DO
3079
        nbasis = nbasis + pmax(i)-1
3080
      END DO
3081
    END DO
3082

3083
    DO i=2,5 ! tri face
3084
      node1 = facedir(1,i)
3085
      node2 = facedir(2,i)
3086
      node3 = facedir(3,i)
3087
      DO j=0,pmax(i)-3
3088
        DO k=1,pmax(i)-2-j
3089
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Pa, Pb, Pc)
3090
          DO l=1,nvec
3091
            Pa = fval(l,node1)
3092
            Pb = fval(l,node2)
3093
            Pc = fval(l,node3)
3094

3095
            La = H1Basis_PyramidTL(node1,u(l),v(l),w(l))
3096
            Lb = H1Basis_PyramidTL(node2,u(l),v(l),w(l))
3097
            Lc = H1Basis_PyramidTL(node3,u(l),v(l),w(l))
3098
            fval(l,nbasis+k) = Pa*Pb*Pc*H1Basis_LegendreP(j,Lb-La)*H1Basis_LegendreP(k-1,2*Lc-1)
3099
          END DO
3100
        END DO
3101
        nbasis = nbasis + MAX(pmax(i)-j-2,0)
3102
      END DO
3103
    END DO
3104
  END SUBROUTINE H1Basis_PyramidFaceP
3105

3106

3107
  SUBROUTINE H1Basis_dPyramidFaceP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, facedir)
3108
    IMPLICIT NONE
3109

3110
    INTEGER, INTENT(IN) :: nvec
3111
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3112
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
3113
    INTEGER, INTENT(IN) :: nbasismax
3114
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
3115
    INTEGER, INTENT(INOUT) :: nbasis
3116
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
3117

3118
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax) :: n
3119
    REAL(KIND=dp) :: La, Lb, Lc, Pa, Pb, Pc, Legi, Legj
3120
    REAL(KIND=dp), DIMENSION(3) :: dLa, dLb, dLc, dPa, dPb, dPc, dLegi, dLegj
3121
    INTEGER :: i,j,k,l, node1, node2, node3, node4, nnb
3122
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64
3123

3124
    nnb = 0; CALL H1Basis_PyramidNodalP(nvec,u,v,w,nbasismax,n,nnb)
3125

3126
    DO i=1,1 ! square face
3127
      node1 = facedir(1,i)
3128
      node2 = facedir(2,i)
3129
      node3 = facedir(3,i)
3130
      node4 = facedir(4,i)
3131

3132
      dLa = H1Basis_dPyramidL(node1)
3133
      dLb = H1Basis_dPyramidL(node2)
3134
      dLc = H1Basis_dPyramidL(node4)
3135

3136
      DO j=0,pmax(i)-2
3137
        DO k=1,pmax(i)-1
3138
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Pa, Pb, dPa, dPb, Legi, Legj, dLegi, dLegj)
3139
          DO l=1,nvec
3140
            Pa = n(l,node1)
3141
            Pb = n(l,node3)
3142

3143
            dPa = grad(l,node1,:)
3144
            dPb = grad(l,node3,:)
3145

3146
            La = H1Basis_PyramidL(node1,u(l),v(l))
3147
            Lb = H1Basis_PyramidL(node2,u(l),v(l))
3148
            Lc = H1Basis_PyramidL(node4,u(l),v(l))
3149

3150
            Legi = H1Basis_LegendreP(j,Lb-La)
3151
            Legj = H1Basis_LegendreP(k-1,Lc-La)
3152

3153
            dLegi = H1Basis_dLegendreP(j,Lb-La)*(dLb-dLa)
3154
            dLegj = H1Basis_dLegendreP(k-1,Lc-La)*(dLc-dLa)
3155
 
3156
            grad(l,nbasis+k,:) = dPa*Pb*Legi*Legj + Pa*dPb*Legi*Legj + &
3157
                       Pa*Pb*dLegi*Legj + Pa*Pb*Legi*dLegj
3158
          END DO
3159
        END DO
3160
        nbasis = nbasis + pmax(i)-1
3161
      END DO
3162
    END DO
3163

3164
    DO i=2,5 ! tri face
3165
      node1 = facedir(1,i)
3166
      node2 = facedir(2,i)
3167
      node3 = facedir(3,i)
3168
 
3169
      dLa = H1Basis_dPyramidTL(node1)
3170
      dLb = H1Basis_dPyramidTL(node2)
3171
      dLc = H1Basis_dPyramidTL(node3)
3172

3173
      DO j=0,pmax(i)-3
3174
        DO k=1,pmax(i)-2-j
3175
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Pa, Pb, Pc, dPa, dPb, dPc, Legi, Legj, dLegi, dLegj)
3176
          DO l=1,nvec
3177
            Pa = n(l,node1)
3178
            Pb = n(l,node2)
3179
            Pc = n(l,node3)
3180

3181
            dPa = grad(l,node1,:)
3182
            dPb = grad(l,node2,:)
3183
            dPc = grad(l,node3,:)
3184

3185
            La = H1Basis_PyramidTL(node1,u(l),v(l),w(l))
3186
            Lb = H1Basis_PyramidTL(node2,u(l),v(l),w(l))
3187
            Lc = H1Basis_PyramidTL(node3,u(l),v(l),w(l))
3188

3189
            Legi = H1Basis_LegendreP(j,Lb-La)
3190
            Legj = H1Basis_LegendreP(k-1,2*Lc-1)
3191

3192
            dLegi = H1Basis_dLegendreP(j,Lb-La)*(dLb-dLa)
3193
            dLegj = H1Basis_dLegendreP(k-1,2*Lc-1)*2*dLc
3194
 
3195
            grad(l,nbasis+k,:) = dPa*Pb*Pc*Legi*Legj + Pa*dPb*Pc*Legi*Legj + &
3196
               Pa*Pb*dPc*Legi*Legj + Pa*Pb*Pc*dLegi*Legj + Pa*Pb*Pc*Legi*dLegj
3197
          END DO
3198
        END DO
3199
        nbasis = nbasis + MAX(pmax(i)-j-2,0)
3200
      END DO
3201
    END DO
3202
  END SUBROUTINE H1Basis_dPyramidFaceP
3203

3204
  SUBROUTINE H1Basis_PyramidBubbleP(nvec, u, v, w, pmax, nbasismax, fval, nbasis)
3205
    IMPLICIT NONE
3206

3207
    INTEGER, INTENT(IN) :: nvec
3208
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3209
    INTEGER, INTENT(IN) :: pmax
3210
    INTEGER, INTENT(IN) :: nbasismax
3211
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
3212
    INTEGER, INTENT(INOUT) :: nbasis
3213

3214
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax) :: n
3215
    REAL(KIND=dp) :: La, Lb, Lc, Pa, Pb, Pc, Legi, Legj, Legk, s
3216
    REAL(KIND=dp), DIMENSION(3) :: dLa, dLb, dLc, dPa, dPb, dPc, dLegi, dLegj
3217
    INTEGER :: i,j,k,l, node1, node2, node3, node4, nnb
3218
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
3219

3220
    ! Calculate value of function
3221
    DO i=0,pmax-3
3222
      DO j=0,pmax-i-3
3223
        DO k=1,pmax-i-j-2
3224
          nbasis = nbasis + 1
3225
          !_ELMER_OMP_SIMD PRIVATE(Pa, Pb, Pc, Legi, Legj, Legk,s)
3226
          DO l=1,nvec
3227
             s = w(l) / SQRT(2._dp)
3228
             Pa = fval(l,1)
3229
             Pb = fval(l,3)
3230
             Pc = fval(l,5)
3231
             Legi = H1Basis_LegendreP(i,u(l))
3232
             Legj = H1Basis_LegendreP(j,v(l))
3233
             Legk = H1Basis_LegendreP(k-1,2*s-1)
3234
             fval(l,nbasis) = Pa*Pb*Pc*Legi*Legj*Legk
3235
          END DO
3236
        END DO
3237
      END DO
3238
    END DO
3239
  END SUBROUTINE H1Basis_PyramidBubbleP
3240

3241
  SUBROUTINE H1Basis_dPyramidBubbleP(nvec, u, v, w, pmax, nbasismax, grad, nbasis)
3242
    IMPLICIT NONE
3243

3244
    INTEGER, INTENT(IN) :: nvec
3245
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3246
    INTEGER, INTENT(IN) :: pmax
3247
    INTEGER, INTENT(IN) :: nbasismax
3248
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
3249
    INTEGER, INTENT(INOUT) :: nbasis
3250

3251
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax) :: n
3252
    REAL(KIND=dp) :: La, Lb, Lc, Pa, Pb, Pc, Legi, Legj, Legk, s
3253
    REAL(KIND=dp), DIMENSION(3) :: dLa, dLb, dLc, dPa, dPb, dPc, dLegi, dLegj, dLegk
3254
    INTEGER :: i,j,k,l, node1, node2, node3, node4, nnb
3255
    REAL(KIND=dp), PARAMETER :: a = 1/SQRT(2._dp)
3256
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64
3257

3258

3259
    nnb = 0; CALL H1Basis_PyramidNodalP(nvec,u,v,w,nbasismax,n,nnb)
3260
    ! Calculate value of function
3261
    DO i=0,pmax-3
3262
      DO j=0,pmax-i-3
3263
        DO k=1,pmax-i-j-2
3264
          nbasis = nbasis + 1
3265
          !_ELMER_OMP_SIMD PRIVATE(Pa, Pb, Pc, Legi, Legj, Legk, dPa, dPb, dPc, dLegi, dLegj, dLegk, s)
3266
          DO l=1,nvec
3267
             s = a * w(l)
3268
             Pa = N(l,1)
3269
             Pb = N(l,3)
3270
             Pc = N(l,5)
3271

3272
             dPa = grad(l,1,:)
3273
             dPb = grad(l,3,:)
3274
             dPc = grad(l,5,:)
3275

3276
             Legi = H1Basis_LegendreP(i,u(l))
3277
             Legj = H1Basis_LegendreP(j,v(l))
3278
             Legk = H1Basis_LegendreP(k-1,2*s-1)
3279

3280
             dLegi = 0; dLegj=0; dLegk=0
3281
             dLegi(1) = H1Basis_dLegendreP(i,u(l))
3282
             dLegj(2) = H1Basis_dLegendreP(j,v(l))
3283
             dLegk(3) = H1Basis_dLegendreP(k-1,2*s-1)*2*a
3284

3285
             grad(l,nbasis,:) = dPa*Pb*Pc*Legi*Legj*Legk + Pa*dPb*Pc*Legi*Legj*Legk + &
3286
                                Pa*Pb*dPc*Legi*Legj*Legk + Pa*Pb*Pc*dLegi*Legj*Legk + &
3287
                                Pa*Pb*Pc*Legi*dLegj*Legk + Pa*Pb*Pc*Legi*Legj*dLegk
3288
          END DO
3289
        END DO
3290
      END DO
3291
    END DO
3292
  END SUBROUTINE H1Basis_dPyramidBubbleP
3293

3294

3295
  SUBROUTINE H1Basis_BrickNodal(nvec, u, v, w, nbasismax, fval, nbasis)
3296
    IMPLICIT NONE
3297

3298
    ! Parameters
3299
    INTEGER, INTENT(IN) :: nvec
3300
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
3301
    INTEGER, INTENT(IN) :: nbasismax
3302
    ! Variables
3303
    REAL(KIND=dp) :: fval(VECTOR_BLOCK_LENGTH,nbasismax)
3304
    INTEGER, INTENT(INOUT) :: nbasis
3305
    
3306
    REAL(Kind=dp), PARAMETER :: c = 1D0/8D0
3307
    INTEGER :: j
3308
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
3309

3310
    !_ELMER_OMP_SIMD
3311
    DO j=1,nvec
3312
      ! Node 1
3313
      fval(j,nbasis+1) = c*(1-u(j))*(1-v(j))*(1-w(j))
3314
      ! Node 2
3315
      fval(j,nbasis+2) = c*(1+u(j))*(1-v(j))*(1-w(j))
3316
      ! Node 3
3317
      fval(j,nbasis+3) = c*(1+u(j))*(1+v(j))*(1-w(j))
3318
      ! Node 4
3319
      fval(j,nbasis+4) = c*(1-u(j))*(1+v(j))*(1-w(j))
3320
      ! Node 5
3321
      fval(j,nbasis+5) = c*(1-u(j))*(1-v(j))*(1+w(j))
3322
      ! Node 6
3323
      fval(j,nbasis+6) = c*(1+u(j))*(1-v(j))*(1+w(j))
3324
      ! Node 7
3325
      fval(j,nbasis+7) = c*(1+u(j))*(1+v(j))*(1+w(j))
3326
      ! Node 8
3327
      fval(j,nbasis+8) = c*(1-u(j))*(1+v(j))*(1+w(j))
3328
    END DO
3329
    nbasis = nbasis + 8
3330
  END SUBROUTINE H1Basis_BrickNodal
3331

3332
  SUBROUTINE H1Basis_dBrickNodal(nvec, u, v, w, nbasismax, grad, nbasis)
3333
    IMPLICIT NONE
3334

3335
    ! Parameters
3336
    INTEGER, INTENT(IN) :: nvec
3337
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
3338
    INTEGER, INTENT(IN) :: nbasismax
3339
    ! Variables
3340
    REAL(Kind=dp) :: grad(VECTOR_BLOCK_LENGTH,nbasismax,3)
3341
    INTEGER, INTENT(INOUT) :: nbasis
3342
    REAL(Kind=dp), PARAMETER :: c = 1D0/8D0
3343
    INTEGER :: j
3344
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
3345

3346
    ! First coordinate (xi)
3347
    !_ELMER_OMP_SIMD
3348
    DO j=1,nvec
3349
      grad(j,nbasis+1,1) = -c*(1-v(j))*(1-w(j))
3350
      grad(j,nbasis+2,1) =  c*(1-v(j))*(1-w(j))
3351
      grad(j,nbasis+3,1) =  c*(1+v(j))*(1-w(j))
3352
      grad(j,nbasis+4,1) = -c*(1+v(j))*(1-w(j))
3353
      grad(j,nbasis+5,1) = -c*(1-v(j))*(1+w(j))
3354
      grad(j,nbasis+6,1) =  c*(1-v(j))*(1+w(j))
3355
      grad(j,nbasis+7,1) =  c*(1+v(j))*(1+w(j))
3356
      grad(j,nbasis+8,1) = -c*(1+v(j))*(1+w(j))
3357
    END DO
3358
    
3359
    ! Second coordinate (eta)
3360
    !_ELMER_OMP_SIMD
3361
    DO j=1,nvec
3362
      grad(j,nbasis+1,2) = -c*(1-u(j))*(1-w(j))
3363
      grad(j,nbasis+2,2) = -c*(1+u(j))*(1-w(j))
3364
      grad(j,nbasis+3,2) =  c*(1+u(j))*(1-w(j))
3365
      grad(j,nbasis+4,2) =  c*(1-u(j))*(1-w(j))
3366
      grad(j,nbasis+5,2) = -c*(1-u(j))*(1+w(j))
3367
      grad(j,nbasis+6,2) = -c*(1+u(j))*(1+w(j))
3368
      grad(j,nbasis+7,2) =  c*(1+u(j))*(1+w(j))
3369
      grad(j,nbasis+8,2) =  c*(1-u(j))*(1+w(j))
3370
    END DO
3371
    
3372
    ! Third coordinate (zeta)
3373
    !_ELMER_OMP_SIMD
3374
    DO j=1,nvec
3375
      grad(j,nbasis+1,3) = -c*(1-u(j))*(1-v(j))
3376
      grad(j,nbasis+2,3) = -c*(1+u(j))*(1-v(j))
3377
      grad(j,nbasis+3,3) = -c*(1+u(j))*(1+v(j))
3378
      grad(j,nbasis+4,3) = -c*(1-u(j))*(1+v(j))
3379
      grad(j,nbasis+5,3) =  c*(1-u(j))*(1-v(j))
3380
      grad(j,nbasis+6,3) =  c*(1+u(j))*(1-v(j))
3381
      grad(j,nbasis+7,3) =  c*(1+u(j))*(1+v(j))
3382
      grad(j,nbasis+8,3) =  c*(1-u(j))*(1+v(j))
3383
    END DO
3384
        nbasis = nbasis + 8
3385
  END SUBROUTINE H1Basis_dBrickNodal
3386

3387
  FUNCTION H1Basis_BrickL(node, u, v, w) RESULT(fval)
3388
    IMPLICIT NONE
3389

3390
    INTEGER, INTENT(IN) :: node
3391
    REAL(KIND=dp), INTENT(IN) :: u, v, w
3392
    REAL(KIND=dp) :: fval
3393
    REAL(KIND=dp), PARAMETER :: c = 1/2D0
3394
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) _ELMER_LINEAR_REF(u) &
3395
    !_ELMER_OMP _ELMER_LINEAR_REF(v) _ELMER_LINEAR_REF(w) NOTINBRANCH
3396
    
3397
    SELECT CASE (node)
3398
    CASE (1)
3399
      fval = c*(3-u-v-w)
3400
    CASE (2)
3401
      fval = c*(3+u-v-w)
3402
    CASE (3)
3403
      fval = c*(3+u+v-w)
3404
    CASE (4)
3405
      fval = c*(3-u+v-w)
3406
    CASE (5)
3407
      fval = c*(3-u-v+w)
3408
    CASE (6)
3409
      fval = c*(3+u-v+w)
3410
    CASE (7)
3411
      fval = c*(3+u+v+w)
3412
    CASE (8)
3413
      fval = c*(3-u+v+w)
3414
    END SELECT
3415
  END FUNCTION H1Basis_BrickL
3416

3417
  FUNCTION H1Basis_dBrickL(node) RESULT(grad)
3418
    IMPLICIT NONE
3419

3420
    INTEGER, INTENT(IN) :: node
3421
    ! REAL(KIND=dp), INTENT(IN) :: u, v
3422
    REAL(KIND=dp) :: grad(3)
3423
    REAL(KIND=dp), PARAMETER :: c = 1/2D0
3424
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(node) NOTINBRANCH
3425
    
3426
    SELECT CASE(node)
3427
    CASE (1)
3428
      grad(1:3) = c*[-1,-1,-1 ]
3429
    CASE (2)
3430
      grad(1:3) = c*[ 1,-1,-1 ]
3431
    CASE (3)
3432
      grad(1:3) = c*[ 1, 1,-1 ]
3433
    CASE (4)
3434
      grad(1:3) = c*[-1, 1,-1 ]
3435
    CASE (5)
3436
      grad(1:3) = c*[-1,-1, 1 ]
3437
    CASE (6)
3438
      grad(1:3) = c*[ 1,-1, 1 ]
3439
    CASE (7)
3440
      grad(1:3) = c*[ 1, 1, 1 ]
3441
    CASE (8)
3442
      grad(1:3) = c*[-1, 1, 1 ]
3443
    END SELECT
3444
  END FUNCTION H1Basis_dBrickL
3445

3446
  SUBROUTINE H1Basis_BrickEdgeL(edge, u, v, w, La, Lb)
3447
    IMPLICIT NONE
3448

3449
    INTEGER, INTENT(IN) :: edge
3450
    REAL(KIND=dp), INTENT(IN) :: u, v, w
3451
    REAL(KIND=dp), INTENT(OUT) :: La, Lb
3452
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(edge) _ELMER_LINEAR_REF(u) _ELMER_LINEAR_REF(v) &
3453
    !_ELMER_OMP _ELMER_LINEAR_REF(w) _ELMER_LINEAR_REF(La) _ELMER_LINEAR_REF(Lb) NOTINBRANCH
3454
    
3455
    SELECT CASE(edge)
3456
    CASE (1)
3457
      La=1-v
3458
      Lb=1-w
3459
    CASE (2)
3460
      La=1+u
3461
      Lb=1-w
3462
    CASE (3)
3463
      La=1+v
3464
      Lb=1-w
3465
    CASE (4)
3466
      La=1-u
3467
      Lb=1-w
3468
    CASE (5)
3469
      La=1-v
3470
      Lb=1+w
3471
    CASE (6)
3472
      La=1+u
3473
      Lb=1+w
3474
    CASE (7)
3475
      La=1+v
3476
      Lb=1+w
3477
    CASE (8)
3478
      La=1-u
3479
      Lb=1+w
3480
    CASE (9)
3481
      La=1-u
3482
      Lb=1-v
3483
    CASE (10)
3484
      La=1+u
3485
      Lb=1-v
3486
    CASE (11)
3487
      La=1+u
3488
      Lb=1+v
3489
    CASE (12)
3490
      La=1-u
3491
      Lb=1+v
3492
    END SELECT
3493
  END SUBROUTINE H1Basis_BrickEdgeL
3494

3495
  SUBROUTINE H1Basis_dBrickEdgeL(edge, dLa, dLb)
3496
    IMPLICIT NONE
3497

3498
    INTEGER, INTENT(IN) :: edge
3499
    REAL(KIND=dp), INTENT(OUT) :: dLa(3), dLb(3)
3500
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(edge) NOTINBRANCH
3501
    
3502
    SELECT CASE(edge)
3503
    CASE (1)
3504
      dLa=[ 0,-1, 0 ] ! 1-v
3505
      dLb=[ 0, 0,-1 ] ! 1-w
3506
    CASE (2)
3507
      dLa=[ 1, 0, 0 ] ! 1+u
3508
      dLb=[ 0, 0,-1 ] ! 1-w
3509
    CASE (3)
3510
      dLa=[ 0, 1, 0 ] ! 1+v
3511
      dLb=[ 0, 0,-1 ] ! 1-w
3512
    CASE (4)
3513
      dLa=[-1, 0, 0 ] ! 1-u
3514
      dLb=[ 0, 0,-1 ] ! 1-w
3515
    CASE (5)
3516
      dLa=[ 0,-1, 0 ] ! 1-v
3517
      dLb=[ 0, 0, 1 ] ! 1+w
3518
    CASE (6)
3519
      dLa=[ 1, 0, 0 ] ! 1+u
3520
      dLb=[ 0, 0, 1 ] ! 1+w
3521
    CASE (7)
3522
      dLa=[ 0, 1, 0 ] ! 1+v
3523
      dLb=[ 0, 0, 1 ] ! 1+w
3524
    CASE (8)
3525
      dLa=[-1, 0, 0 ] ! 1-u
3526
      dLb=[ 0, 0, 1 ] ! 1+w
3527
    CASE (9)
3528
      dLa=[-1, 0, 0 ] ! 1-u
3529
      dLb=[ 0,-1, 0 ] ! 1-v
3530
    CASE (10)
3531
      dLa=[ 1, 0, 0 ] ! 1+u
3532
      dLb=[ 0,-1, 0 ] ! 1-v
3533
    CASE (11)
3534
      dLa=[ 1, 0, 0 ] ! 1+u
3535
      dLb=[ 0, 1, 0 ] ! 1+v
3536
    CASE (12)
3537
      dLa=[-1, 0, 0 ] ! 1-u
3538
      dLb=[ 0, 1, 0 ] ! 1+v
3539
    END SELECT
3540
  END SUBROUTINE H1Basis_dBrickEdgeL
3541

3542
! --- start serendipity brick
3543

3544
  
3545
  SUBROUTINE H1Basis_SD_BrickEdgeP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, edgedir)
3546
    IMPLICIT NONE
3547

3548
    INTEGER, INTENT(IN) :: nvec
3549
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3550
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
3551
    INTEGER, INTENT(IN) :: nbasismax
3552
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
3553
    INTEGER, INTENT(INOUT) :: nbasis
3554
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
3555

3556
    REAL(KIND=dp), PARAMETER :: c = 1/4D0
3557
    REAL(KIND=dp) :: La, Lb, Aa, Ba
3558
    INTEGER :: i,j,k,l
3559
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
3560

3561
    ! For each edge
3562
    DO i=1,12
3563
      DO j=1,pmax(i)-1
3564
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Aa, Ba)
3565
        DO k=1,nvec
3566
          La=H1Basis_BrickL(edgedir(1,i), u(k), v(k), w(k))
3567
          Lb=H1Basis_BrickL(edgedir(2,i), u(k), v(k), w(k))
3568
          CALL H1Basis_BrickEdgeL(i, u(k), v(k), w(k), Aa, Ba)
3569
          fval(k, nbasis+j) = c*H1Basis_Phi(j+1, Lb-La)*Aa*Ba
3570
        END DO
3571
      END DO
3572
      ! nbasis = nbasis + (pmax(i)-2) + 1
3573
      nbasis = nbasis + pmax(i) - 1
3574
    END DO
3575
  END SUBROUTINE H1Basis_SD_BrickEdgeP
3576
  
3577
  SUBROUTINE H1Basis_SD_dBrickEdgeP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, edgedir)
3578
    IMPLICIT NONE
3579

3580
    INTEGER, INTENT(IN) :: nvec
3581
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3582
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
3583
    INTEGER, INTENT(IN) :: nbasismax
3584
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
3585
    INTEGER, INTENT(INOUT) :: nbasis
3586
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
3587

3588
    REAL(KIND=dp), PARAMETER :: c = 1/4D0
3589
    REAL(KIND=dp) :: La, Lb, Aa, Ba, Phi, dPhi, dLa(3), &
3590
            dLb(3), dAa(3), dBa(3)
3591
    INTEGER :: i,j,k
3592
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
3593

3594
    ! For each edge
3595
    DO i=1,12
3596
      dLa = H1Basis_dBrickL(edgedir(1,i))
3597
      dLb = H1Basis_dBrickL(edgedir(2,i))
3598
      CALL H1Basis_dBrickEdgeL(i, dAa, dBa)
3599

3600
      DO j=1,pmax(i)-1
3601
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Aa, Ba, Phi, dPhi)
3602
        DO k=1,nvec
3603
          La=H1Basis_BrickL(edgedir(1,i), u(k), v(k), w(k))
3604
          Lb=H1Basis_BrickL(edgedir(2,i), u(k), v(k), w(k))
3605
          CALL H1Basis_BrickEdgeL(i, u(k), v(k), w(k), Aa, Ba)
3606

3607
          Phi = H1Basis_Phi(j+1, Lb-La)
3608
          dPhi = H1Basis_dPhi(j+1, Lb-La)
3609

3610
          grad(k,nbasis+j,1) = c*dPhi*(dLb(1)-dLa(1))*Aa*Ba +&
3611
                  c*Phi*dAa(1)*Ba + c*Phi*Aa*dBa(1)
3612
          grad(k,nbasis+j,2) = c*dPhi*(dLb(2)-dLa(2))*Aa*Ba +&
3613
                  c*Phi*dAa(2)*Ba + c*Phi*Aa*dBa(2)
3614
          grad(k,nbasis+j,3) = c*dPhi*(dLb(3)-dLa(3))*Aa*Ba +&
3615
                  c*Phi*dAa(3)*Ba + c*Phi*Aa*dBa(3)
3616
        END DO
3617
      END DO
3618
      ! nbasis = nbasis + (pmax(i)-2) + 1
3619
      nbasis = nbasis + pmax(i) - 1
3620
    END DO
3621
  END SUBROUTINE H1Basis_SD_dBrickEdgeP
3622
  
3623
  SUBROUTINE H1Basis_SD_BrickFaceP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, facedir)
3624
    IMPLICIT NONE
3625

3626
    INTEGER, INTENT(IN) :: nvec
3627
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3628
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
3629
    INTEGER, INTENT(IN) :: nbasismax
3630
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
3631
    INTEGER, INTENT(INOUT) :: nbasis
3632
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
3633

3634
    REAL(KIND=dp) :: La, Lb, Lc, Ld
3635
    REAL(KIND=dp), PARAMETER :: a = 1D0/4
3636
    INTEGER :: i,j,k,l
3637
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
3638

3639
    ! For each face
3640
    DO i=1,6
3641
      DO j=2,pmax(i)
3642
        DO k=1,pmax(i)-j-1
3643
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Ld)
3644
          DO l=1,nvec
3645
            La=H1Basis_BrickL(facedir(1,i), u(l), v(l), w(l))
3646
            Lb=H1Basis_BrickL(facedir(2,i), u(l), v(l), w(l))
3647
            Lc=H1Basis_BrickL(facedir(3,i), u(l), v(l), w(l))
3648
            Ld=H1Basis_BrickL(facedir(4,i), u(l), v(l), w(l))
3649
            
3650
            fval(l, nbasis+k) = (a*(La+Lb+Lc+Ld)-1)*H1Basis_Phi(j, Lb-La) &
3651
                               *H1Basis_Phi(k+1, Ld-La)
3652
          END DO
3653
        END DO
3654
        ! nbasis = nbasis + (pmax(i)-j-2) + 1
3655
        nbasis = nbasis + MAX(pmax(i)-j-1,0)
3656
      END DO
3657
    END DO
3658
  END SUBROUTINE H1Basis_SD_BrickFaceP
3659
  
3660
  SUBROUTINE H1Basis_SD_dBrickFaceP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, facedir)
3661
    IMPLICIT NONE
3662

3663
    INTEGER, INTENT(IN) :: nvec
3664
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3665
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
3666
    INTEGER, INTENT(IN) :: nbasismax
3667
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
3668
    INTEGER, INTENT(INOUT) :: nbasis
3669
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
3670

3671
    REAL(KIND=dp) :: La, Lb, Lc, Ld, dLa(3), dLb(3), dLc(3), dLd(3), &
3672
          PhiLbLa, PhiLdLa, LaLbLcLd
3673
    REAL(KIND=dp), PARAMETER :: a=1D0/4D0
3674
    INTEGER :: i,j,k,l
3675
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
3676

3677
    ! For each face
3678
    DO i=1,6
3679
      dLa=H1Basis_dBrickL(facedir(1,i))
3680
      dLb=H1Basis_dBrickL(facedir(2,i))
3681
      dLc=H1Basis_dBrickL(facedir(3,i))
3682
      dLd=H1Basis_dBrickL(facedir(4,i))
3683
      DO j=2,pmax(i)
3684
        DO k=1,pmax(i)-j-1
3685
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Lc, Ld, PhiLbLa, PhiLdLa, LaLbLcLd)
3686
          DO l=1,nvec
3687
            La=H1Basis_BrickL(facedir(1,i), u(l), v(l), w(l))
3688
            Lb=H1Basis_BrickL(facedir(2,i), u(l), v(l), w(l))
3689
            Lc=H1Basis_BrickL(facedir(3,i), u(l), v(l), w(l))
3690
            Ld=H1Basis_BrickL(facedir(4,i), u(l), v(l), w(l))
3691
            
3692
            PhiLbLa=H1Basis_Phi(j, Lb-La)
3693
            PhiLdLa=H1Basis_Phi(k+1, Ld-La)
3694
            LaLbLcLd=a*(La+Lb+Lc+Ld)-1
3695

3696
            grad(l, nbasis+k, 1)=a*(dLa(1)+dLb(1)+dLc(1)+dLd(1))*PhiLbLa*PhiLdLa + &
3697
                  LaLbLcLd*H1Basis_dPhi(j,Lb-La)*(dLb(1)-dLa(1))*PhiLdLa + &
3698
                  LaLbLcLd*PhiLbLa*H1Basis_dPhi(k+1, Ld-La)*(dLd(1)-dLa(1))
3699
            grad(l, nbasis+k, 2)=a*(dLa(2)+dLb(2)+dLc(2)+dLd(2))*PhiLbLa*PhiLdLa + &
3700
                  LaLbLcLd*H1Basis_dPhi(j,Lb-La)*(dLb(2)-dLa(2))*PhiLdLa + &
3701
                  LaLbLcLd*PhiLbLa*H1Basis_dPhi(k+1, Ld-La)*(dLd(2)-dLa(2))
3702
            grad(l, nbasis+k, 3)=a*(dLa(3)+dLb(3)+dLc(3)+dLd(3))*PhiLbLa*PhiLdLa + &
3703
                  LaLbLcLd*H1Basis_dPhi(j,Lb-La)*(dLb(3)-dLa(3))*PhiLdLa + &
3704
                  LaLbLcLd*PhiLbLa*H1Basis_dPhi(k+1, Ld-La)*(dLd(3)-dLa(3))
3705
          END DO
3706
        END DO
3707
        ! nbasis = nbasis + (pmax(i)-j-2) + 1
3708
        nbasis = nbasis + MAX(pmax(i)-j-1,0)
3709
      END DO
3710
    END DO
3711
  END SUBROUTINE H1Basis_SD_dBrickFaceP
3712
  
3713
  SUBROUTINE H1Basis_SD_BrickBubbleP(nvec, u, v, w, pmax, nbasismax, fval, nbasis)
3714
    IMPLICIT NONE
3715

3716
    INTEGER, INTENT(IN) :: nvec
3717
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3718
    INTEGER, INTENT(IN) :: pmax
3719
    INTEGER, INTENT(IN) :: nbasismax
3720
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
3721
    INTEGER, INTENT(INOUT) :: nbasis
3722

3723
    INTEGER :: i, j, k, l
3724
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
3725

3726
    ! For each bubble calculate value of basis function and its derivative
3727
    ! for index pairs i,j,k=2,..,p-4, i+j+k=6,..,pmax
3728
    DO i=2,pmax-4
3729
      DO j=2,pmax-i-2
3730
        DO k=1,pmax-i-j-1
3731
          !_ELMER_OMP_SIMD
3732
          DO l=1,nvec
3733
            fval(l,nbasis+k) = H1Basis_Phi(i,u(l))*H1Basis_Phi(j,v(l))*H1Basis_Phi(k+1,w(l))
3734
          END DO
3735
        END DO
3736
        nbasis = nbasis+MAX(pmax-i-j-1,0)
3737
      END DO
3738
    END DO
3739
  END SUBROUTINE H1Basis_SD_BrickBubbleP
3740

3741
  SUBROUTINE H1Basis_SD_dBrickBubbleP(nvec, u, v, w, pmax, nbasismax, grad, nbasis)
3742
    IMPLICIT NONE
3743

3744
    INTEGER, INTENT(IN) :: nvec
3745
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3746
    INTEGER, INTENT(IN) :: pmax
3747
    INTEGER, INTENT(IN) :: nbasismax
3748
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
3749
    INTEGER, INTENT(INOUT) :: nbasis
3750

3751
    ! Variables
3752
    INTEGER :: i, j, k, l
3753
    REAL(KIND=dp) :: phiU, phiV, phiW
3754
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
3755

3756
    ! For each bubble calculate value of basis function and its derivative
3757
    ! for index pairs i,j,k=2,..,p-4, i+j+k=6,..,pmax
3758
    DO i=2,pmax-4
3759
      DO j=2,pmax-i-2
3760
        DO k=1,pmax-i-j-1
3761
          !_ELMER_OMP_SIMD PRIVATE(phiU, phiV, phiW)
3762
          DO l=1,nvec
3763
            phiU = H1Basis_Phi(i,u(l))
3764
            phiV = H1Basis_Phi(j,v(l))
3765
            phiW = H1Basis_Phi(k+1,w(l))
3766
            grad(l,nbasis+k,1) = H1Basis_dPhi(i,u(l))*phiV*phiW
3767
            grad(l,nbasis+k,2) = phiU*H1Basis_dPhi(j,v(l))*phiW
3768
            grad(l,nbasis+k,3) = phiU*phiV*H1Basis_dPhi(k+1,w(l))
3769
          END DO
3770
        END DO
3771
        nbasis = nbasis+MAX(pmax-i-j-1,0)
3772
      END DO
3773
    END DO
3774
  END SUBROUTINE H1Basis_SD_dBrickBubbleP
3775

3776

3777
! --- end serendipity brick
3778
  
3779
  SUBROUTINE H1Basis_BrickEdgeP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, edgedir)
3780
    IMPLICIT NONE
3781

3782
    INTEGER, INTENT(IN) :: nvec
3783
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3784
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
3785
    INTEGER, INTENT(IN) :: nbasismax
3786
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
3787
    INTEGER, INTENT(INOUT) :: nbasis
3788
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
3789

3790
    REAL(KIND=dp) :: La, Lb, Na, Nb, Lp(nvec,8)
3791
    INTEGER :: i,j,k,l, node1, node2, nnb
3792
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
3793

3794
    DO i=1,8
3795
      !_ELMER_OMP_SIMD
3796
      DO k=1,nvec
3797
        Lp(k,i) = H1Basis_BrickL(i, u(k), v(k), w(k))
3798
      END DO
3799
    END DO
3800

3801
    ! For each edge
3802
    DO i=1,12
3803
      node1 = edgedir(1,i)
3804
      node2 = edgedir(2,i)
3805
      DO j=1,pmax(i)-1
3806
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb)
3807
        DO k=1,nvec
3808
          Na = fval(k, node1)
3809
          Nb = fval(k, node2)
3810
          La = Lp(k, node1)
3811
          Lb = Lp(k, node2)
3812
          fval(k, nbasis+j) = Na*Nb*H1Basis_varPhi(j+1, Lb-La)
3813
        END DO
3814
      END DO
3815
      nbasis = nbasis + pmax(i) - 1
3816
    END DO
3817
  END SUBROUTINE H1Basis_BrickEdgeP
3818
  
3819
  SUBROUTINE H1Basis_dBrickEdgeP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, edgedir)
3820
    IMPLICIT NONE
3821

3822
    INTEGER, INTENT(IN) :: nvec
3823
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3824
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
3825
    INTEGER, INTENT(IN) :: nbasismax
3826
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
3827
    INTEGER, INTENT(INOUT) :: nbasis
3828
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: edgedir
3829
    REAL(KIND=dp) :: La, Lb, Na, Nb, Phi, dPhi, dLa(3), Lp(nvec,8), &
3830
            dLb(3), dNa(3), dNb(3), N(VECTOR_BLOCK_LENGTH, nbasismax)
3831
    INTEGER :: i,j,k, node1, node2, nnb
3832
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
3833

3834
    nnb = 0; CALL H1Basis_BrickNodal(nvec, u, v, w, nbasismax, n, nnb )
3835

3836
    DO i=1,8
3837
      !_ELMER_OMP_SIMD
3838
      DO k=1,nvec
3839
        Lp(k,i) = H1Basis_BrickL(i, u(k), v(k), w(k))
3840
      END DO
3841
    END DO
3842

3843
    ! For each edge
3844
    DO i=1,12
3845
      node1 = edgedir(1,i)
3846
      node2 = edgedir(2,i)
3847

3848
      dLa = H1Basis_dBrickL(node1)
3849
      dLb = H1Basis_dBrickL(node2)
3850

3851
      DO j=1,pmax(i)-1
3852
        !_ELMER_OMP_SIMD PRIVATE(La, Lb, Na, Nb, Phi, dPhi)
3853
        DO k=1,nvec
3854
          La = Lp(k,node1)
3855
          Lb = Lp(k,node2)
3856

3857
          Na = N(k,node1)
3858
          Nb = N(k,node2)
3859
          dNa = grad(k,node1,:)
3860
          dNb = grad(k,node2,:)
3861

3862
          Phi = H1Basis_varPhi(j+1, Lb-La)
3863
          dPhi = H1Basis_dvarPhi(j+1, Lb-La)
3864

3865
          grad(k,nbasis+j,1) = dPhi*(dLb(1)-dLa(1))*Na*Nb + &
3866
                  Phi*dNa(1)*Nb + Phi*Na*dNb(1)
3867
          grad(k,nbasis+j,2) = dPhi*(dLb(2)-dLa(2))*Na*Nb + &
3868
                  Phi*dNa(2)*Nb + Phi*Na*dNb(2)
3869
          grad(k,nbasis+j,3) = dPhi*(dLb(3)-dLa(3))*Na*Nb + &
3870
                  Phi*dNa(3)*Nb + Phi*Na*dNb(3)
3871
        END DO
3872
      END DO
3873
      ! nbasis = nbasis + (pmax(i)-2) + 1
3874
      nbasis = nbasis + pmax(i) - 1
3875
    END DO
3876
  END SUBROUTINE H1Basis_dBrickEdgeP
3877
  
3878
  SUBROUTINE H1Basis_BrickFaceP(nvec, u, v, w, pmax, nbasismax, fval, nbasis, facedir)
3879
    IMPLICIT NONE
3880

3881
    INTEGER, INTENT(IN) :: nvec
3882
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3883
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
3884
    INTEGER, INTENT(IN) :: nbasismax
3885
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
3886
    INTEGER, INTENT(INOUT) :: nbasis
3887
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
3888

3889
    REAL(KIND=dp) :: La, Lb, Lc, Ld, Na, Nb, Lp(nvec,8)
3890
    INTEGER :: i,j,k,l, node1, node2, node3, node4
3891
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
3892

3893
    DO i=1,8
3894
      !_ELMER_OMP_SIMD
3895
      DO k=1,nvec
3896
        Lp(k,i) = H1Basis_BrickL(i, u(k), v(k), w(k))
3897
      END DO
3898
    END DO
3899

3900
    ! For each face
3901
    DO i=1,6
3902
      node1 = facedir(1,i)
3903
      node2 = facedir(2,i)
3904
      node3 = facedir(3,i)
3905
      node4 = facedir(4,i)
3906
      DO j=0,pmax(i)-2
3907
        DO k=1,pmax(i)-1
3908
          !_ELMER_OMP_SIMD PRIVATE(Na,Nb,La, Lb, Ld)
3909
          DO l=1,nvec
3910
            Na = fval(l,node1)
3911
            Nb = fval(l,node3)
3912
            La = Lp(l,node1)
3913
            Lb = Lp(l,node2)
3914
            Ld = Lp(l,node4)
3915
            fval(l, nbasis+k) = Na*Nb*H1Basis_LegendreP(j, Lb-La)*H1Basis_LegendreP(k-1, Ld-La)
3916
          END DO
3917
        END DO
3918
        ! nbasis = nbasis + (pmax(i)-j-2) + 1
3919
        nbasis = nbasis + pmax(i)-1
3920
      END DO
3921
    END DO
3922
  END SUBROUTINE H1Basis_BrickFaceP
3923
  
3924
  SUBROUTINE H1Basis_dBrickFaceP(nvec, u, v, w, pmax, nbasismax, grad, nbasis, facedir)
3925
    IMPLICIT NONE
3926

3927
    INTEGER, INTENT(IN) :: nvec
3928
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3929
    INTEGER, DIMENSION(:) CONTIG, INTENT(IN) :: pmax
3930
    INTEGER, INTENT(IN) :: nbasismax
3931
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
3932
    INTEGER, INTENT(INOUT) :: nbasis
3933
    INTEGER, DIMENSION(:,:) CONTIG, INTENT(IN) :: facedir
3934

3935
    REAL(KIND=dp) :: La, Lb, Lc, Ld, dLa(3), dLb(3), dLc(3), dLd(3), Lp(nvec,8), &
3936
          Na,Nb,dNa(3),dNb(3), N(VECTOR_BLOCK_LENGTH, nbasismax), PhiU, PhiV, dPhiU, dPhiV
3937
    INTEGER :: i,j,k,l, nnb, node1, node2, node3, node4
3938
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
3939

3940
    nnb = 0; CALL H1Basis_BrickNodal(nvec, u, v, w, nbasismax, n, nnb )
3941
    DO i=1,8
3942
      !_ELMER_OMP_SIMD
3943
      DO k=1,nvec
3944
        Lp(k,i) = H1Basis_BrickL(i, u(k), v(k), w(k))
3945
      END DO
3946
    END DO
3947

3948
    ! For each face
3949
    DO i=1,6
3950
      node1 = facedir(1,i)
3951
      node2 = facedir(2,i)
3952
      node3 = facedir(3,i)
3953
      node4 = facedir(4,i)
3954
      dLa = H1Basis_dBrickL(node1)
3955
      dLb = H1Basis_dBrickL(node2)
3956
      dLc = H1Basis_dBrickL(node3)
3957
      dLd = H1Basis_dBrickL(node4)
3958

3959
      DO j=0,pmax(i)-2
3960
        DO k=1,pmax(i)-1
3961
          !_ELMER_OMP_SIMD PRIVATE(La, Lb, Ld, PhiU, PhiV, Na, Nb, dNa, dNb)
3962
          DO l=1,nvec
3963

3964
            La = Lp(l,node1)
3965
            Lb = Lp(l,node2)
3966
            Ld = Lp(l,node4)
3967
            
3968
            PhiU = H1Basis_LegendreP(j, Lb-La)
3969
            PhiV = H1Basis_LegendreP(k-1, Ld-La)
3970

3971
            dPhiU = H1Basis_dLegendreP(j, Lb-La)
3972
            dPhiV = H1Basis_dLegendreP(k-1, Ld-La)
3973

3974
            Na = N(l,node1)
3975
            Nb = N(l,node3)
3976

3977
            dNa = grad(l,node1,:)
3978
            dNb = grad(l,node3,:)
3979

3980
            grad(l, nbasis+k, :) = dNa*Nb*PhiU*PhiV + Na*dNb*PhiU*PhiV + &
3981
                Na*Nb*dPhiU*(dLb-dLa)*PhiV + Na*Nb*PhiU*dPhiV*(dLd-dLa)
3982
          END DO
3983
        END DO
3984
        ! nbasis = nbasis + (pmax(i)-j-2) + 1
3985
        nbasis = nbasis + pmax(i)-1
3986
      END DO
3987
    END DO
3988
  END SUBROUTINE H1Basis_dBrickFaceP
3989
  
3990
  SUBROUTINE H1Basis_BrickBubbleP(nvec, u, v, w, pmax, nbasismax, fval, nbasis)
3991
    IMPLICIT NONE
3992

3993
    INTEGER, INTENT(IN) :: nvec
3994
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
3995
    INTEGER, INTENT(IN) :: pmax
3996
    INTEGER, INTENT(IN) :: nbasismax
3997
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
3998
    INTEGER, INTENT(INOUT) :: nbasis
3999

4000
    INTEGER :: i, j, k, l
4001
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, fval:64
4002

4003
    ! For each bubble calculate value of basis function and its derivative
4004
    ! for index pairs i,j,k=2,..,p-4, i+j+k=6,..,pmax
4005
    DO i=2,pmax
4006
      DO j=2,pmax
4007
        DO k=1,pmax-1
4008
          !_ELMER_OMP_SIMD
4009
          DO l=1,nvec
4010
            fval(l,nbasis+k) = H1Basis_Phi(i,u(l))*H1Basis_Phi(j,v(l))*H1Basis_Phi(k+1,w(l))
4011
          END DO
4012
        END DO
4013
        nbasis = nbasis + pmax-1
4014
      END DO
4015
    END DO
4016
  END SUBROUTINE H1Basis_BrickBubbleP
4017

4018
  SUBROUTINE H1Basis_dBrickBubbleP(nvec, u, v, w, pmax, nbasismax, grad, nbasis)
4019
    IMPLICIT NONE
4020

4021
    INTEGER, INTENT(IN) :: nvec
4022
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u, v, w
4023
    INTEGER, INTENT(IN) :: pmax
4024
    INTEGER, INTENT(IN) :: nbasismax
4025
    REAL(KIND=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
4026
    INTEGER, INTENT(INOUT) :: nbasis
4027

4028
    ! Variables
4029
    INTEGER :: i, j, k, l
4030
    REAL(KIND=dp) :: phiU, phiV, phiW
4031
!DIR$ ASSUME_ALIGNED u:64, v:64, w:64, grad:64
4032

4033
    ! For each bubble calculate value of basis function and its derivative
4034
    ! for index pairs i,j,k=2,..,p-4, i+j+k=6,..,pmax
4035
    DO i=2,pmax
4036
      DO j=2,pmax
4037
        DO k=1,pmax-1
4038
          !_ELMER_OMP_SIMD PRIVATE(phiU, phiV, phiW)
4039
          DO l=1,nvec
4040
            phiU = H1Basis_Phi(i,u(l))
4041
            phiV = H1Basis_Phi(j,v(l))
4042
            phiW = H1Basis_Phi(k+1,w(l))
4043

4044
            grad(l,nbasis+k,1) = H1Basis_dPhi(i,u(l))*phiV*phiW
4045
            grad(l,nbasis+k,2) = phiU*H1Basis_dPhi(j,v(l))*phiW
4046
            grad(l,nbasis+k,3) = phiU*phiV*H1Basis_dPhi(k+1,w(l))
4047
          END DO
4048
        END DO
4049
        nbasis = nbasis+pmax-1
4050
      END DO
4051
    END DO
4052
  END SUBROUTINE H1Basis_dBrickBubbleP
4053

4054
  PURE FUNCTION H1Basis_Phi(k, x) RESULT(fval)
4055
    IMPLICIT NONE
4056

4057
    ! Parameters
4058
    INTEGER, INTENT(IN) :: k
4059
    REAL (KIND=dp), INTENT(IN) :: x
4060
    ! Return value
4061
    REAL (KIND=dp) :: fval
4062
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(k) _ELMER_LINEAR_REF(x) NOTINBRANCH
4063
    
4064
    ! Phi function values (autogenerated to Horner form)
4065
    SELECT CASE(k)
4066
    CASE(2)
4067
      fval = -SQRT(0.6D1) / 0.4D1 + SQRT(0.6D1) * x ** 2 / 0.4D1
4068
    CASE(3)
4069
      fval = (-SQRT(0.10D2) / 0.4D1 + SQRT(0.10D2) * x ** 2 / 0.4D1) * x
4070
    CASE(4)
4071
      fval = SQRT(0.14D2) / 0.16D2 + (-0.3D1 / 0.8D1 * SQRT(0.14D2) &
4072
              + 0.5D1 / 0.16D2 * SQRT(0.14D2) * x ** 2) * x ** 2
4073
    CASE(5)
4074
      fval = (0.9D1 / 0.16D2 * SQRT(0.2D1) + (-0.15D2 / 0.8D1 * SQRT(0.2D1)&
4075
              + 0.21D2 / 0.16D2 * SQRT(0.2D1) * x ** 2) * x ** 2) * x
4076
    CASE(6)
4077
      fval = -SQRT(0.22D2) / 0.32D2 + (0.15D2 / 0.32D2 * SQRT(0.22D2) + &
4078
              (-0.35D2 / 0.32D2 * SQRT(0.22D2) + 0.21D2 / 0.32D2 * &
4079
              SQRT(0.22D2) * x ** 2) * x ** 2) * x ** 2
4080
    CASE(7)
4081
      fval = (-0.5D1 / 0.32D2 * SQRT(0.26D2) + (0.35D2 / 0.32D2 * SQRT(0.26D2) &
4082
              + (-0.63D2 / 0.32D2 * SQRT(0.26D2) + 0.33D2 / 0.32D2 * &
4083
              SQRT(0.26D2) * x ** 2) * x ** 2) * x ** 2) * x
4084
    CASE(8)
4085
      fval = 0.5D1 / 0.256D3 * SQRT(0.30D2) + (-0.35D2 / 0.64D2 * SQRT(0.30D2) &
4086
              + (0.315D3 / 0.128D3 * SQRT(0.30D2) + (-0.231D3 / 0.64D2 * SQRT(0.30D2) &
4087
              + 0.429D3 / 0.256D3 * SQRT(0.30D2) * x ** 2) * x ** 2) * x ** 2) * x ** 2
4088
    CASE(9)
4089
      fval = (0.35D2 / 0.256D3 * SQRT(0.34D2) + (-0.105D3 / 0.64D2 * SQRT(0.34D2) +&
4090
              (0.693D3 / 0.128D3 * SQRT(0.34D2) + (-0.429D3 / 0.64D2 * SQRT(0.34D2) +&
4091
              0.715D3 / 0.256D3 * SQRT(0.34D2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x
4092
    CASE(10)
4093
      fval = -0.7D1 / 0.512D3 * SQRT(0.38D2) + (0.315D3 / 0.512D3 * SQRT(0.38D2) + &
4094
              (-0.1155D4 / 0.256D3 * SQRT(0.38D2) + (0.3003D4 / 0.256D3 * SQRT(0.38D2) +&
4095
              (-0.6435D4 / 0.512D3 * SQRT(0.38D2) + 0.2431D4 / 0.512D3 * SQRT(0.38D2) &
4096
              * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2
4097
    CASE(11)
4098
      fval = (-0.63D2 / 0.512D3 * SQRT(0.42D2) + (0.1155D4 / 0.512D3 * SQRT(0.42D2) + &
4099
              (-0.3003D4 / 0.256D3 * SQRT(0.42D2) + (0.6435D4 / 0.256D3 * SQRT(0.42D2) +&
4100
              (-0.12155D5 / 0.512D3 * SQRT(0.42D2) + 0.4199D4 / 0.512D3 * &
4101
              SQRT(0.42D2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x
4102
    CASE(12)
4103
      fval = 0.21D2 / 0.2048D4 * SQRT(0.46D2) + (-0.693D3 / 0.1024D4 * SQRT(0.46D2) + &
4104
              (0.15015D5 / 0.2048D4 * SQRT(0.46D2) + (-0.15015D5 / 0.512D3 * SQRT(0.46D2) +&
4105
              (0.109395D6 / 0.2048D4 * SQRT(0.46D2) + (-0.46189D5 / 0.1024D4 * SQRT(0.46D2) +&
4106
              0.29393D5 / 0.2048D4 * SQRT(0.46D2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) &
4107
              * x ** 2) * x ** 2
4108
    CASE(13)
4109
      fval = (0.1155D4 / 0.2048D4 * SQRT(0.2D1) + (-0.15015D5 / 0.1024D4 * SQRT(0.2D1) +&
4110
              (0.225225D6 / 0.2048D4 * SQRT(0.2D1) + (-0.182325D6 / 0.512D3 * SQRT(0.2D1) +&
4111
              (0.1154725D7 / 0.2048D4 * SQRT(0.2D1) + (-0.440895D6 / 0.1024D4 * SQRT(0.2D1) +&
4112
              0.260015D6 / 0.2048D4 * SQRT(0.2D1) * x ** 2) * x ** 2) * x ** 2) * x ** 2) &
4113
              * x ** 2) * x ** 2) * x
4114
    CASE(14)
4115
      fval = -0.99D2 / 0.4096D4 * SQRT(0.6D1) + (0.9009D4 / 0.4096D4 * SQRT(0.6D1) + &
4116
              (-0.135135D6 / 0.4096D4 * SQRT(0.6D1) + (0.765765D6 / 0.4096D4 * SQRT(0.6D1) + &
4117
              (-0.2078505D7 / 0.4096D4 * SQRT(0.6D1) + (0.2909907D7 / 0.4096D4 * SQRT(0.6D1) +&
4118
              (-0.2028117D7 / 0.4096D4 * SQRT(0.6D1) + 0.557175D6 / 0.4096D4 * SQRT(0.6D1) &
4119
              * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2
4120
    CASE DEFAULT
4121
      ! TODO: Handle error somehow
4122
    END SELECT
4123
  END FUNCTION H1Basis_Phi
4124

4125
  PURE FUNCTION H1Basis_dPhi(k, x) RESULT(fval)
4126
    IMPLICIT NONE
4127

4128
    ! Parameters
4129
    INTEGER, INTENT(IN) :: k
4130
    REAL (KIND=dp), INTENT(IN) :: x
4131
    ! Return value
4132
    REAL (KIND=dp) :: fval
4133
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(k) _ELMER_LINEAR_REF(x) NOTINBRANCH
4134
    
4135
    ! Phi function values (autogenerated to Horner form)
4136
    SELECT CASE(k)
4137
    CASE(2)
4138
      fval = x * SQRT(0.6D1) / 0.2D1
4139
    CASE(3)
4140
      fval = DBLE(3 * x ** 2 - 1) * SQRT(0.10D2) / 0.4D1
4141
    CASE(4)
4142
      fval = DBLE(5 * x ** 2 - 3) * DBLE(x) * SQRT(0.14D2) / 0.4D1
4143
    CASE(5)
4144
      fval = 0.3D1 / 0.16D2 * DBLE(3 + (35 * x ** 2 - 30) * x ** 2) * SQRT(0.2D1)
4145
    CASE(6)
4146
      fval = DBLE(15 + (63 * x ** 2 - 70) * x ** 2) * DBLE(x) * SQRT(0.22D2) / 0.16D2
4147
    CASE(7)
4148
      fval = DBLE(-5 + (105 + (231 * x ** 2 - 315) * x ** 2) * x ** 2) * &
4149
              SQRT(0.26D2) / 0.32D2
4150
    CASE(8)
4151
      fval = DBLE(-35 + (315 + (429 * x ** 2 - 693) * x ** 2) * x ** 2) * &
4152
              DBLE(x) * SQRT(0.30D2) / 0.32D2
4153
    CASE(9)
4154
      fval = DBLE(35 + (-1260 + (6930 + (6435 * x ** 2 - 12012) * x ** 2) &
4155
              * x ** 2) * x ** 2) * SQRT(0.34D2) / 0.256D3
4156
    CASE(10)
4157
      fval = DBLE(315 + (-4620 + (18018 + (12155 * x ** 2 - 25740) * x ** 2) * &
4158
              x ** 2) * x ** 2) * DBLE(x) * SQRT(0.38D2) / 0.256D3
4159
    CASE(11)
4160
      fval = DBLE(-63 + (3465 + (-30030 + (90090 + (46189 * x ** 2 - 109395) * &
4161
              x ** 2) * x ** 2) * x ** 2) * x ** 2) * SQRT(0.42D2) / 0.512D3
4162
    CASE(12)
4163
      fval = DBLE(-693 + (15015 + (-90090 + (218790 + (88179 * x ** 2 - 230945) *&
4164
              x ** 2) * x ** 2) * x ** 2) * x ** 2) * DBLE(x) * SQRT(0.46D2) / 0.512D3
4165
    CASE(13)
4166
      fval = 0.5D1 / 0.2048D4 * DBLE(231 + (-18018 + (225225 + (-1021020 + &
4167
              (2078505 + (676039 * x ** 2 - 1939938) * x ** 2) * x ** 2) * &
4168
              x ** 2) * x ** 2) * x ** 2) * SQRT(0.2D1)
4169
    CASE(14)
4170
      fval = 0.3D1 / 0.2048D4 * DBLE(3003 + (-90090 + (765765 + (-2771340 + &
4171
              (4849845 + (1300075 * x ** 2 - 4056234) * x ** 2) * x ** 2) * &
4172
              x ** 2) * x ** 2) * x ** 2) * DBLE(x) * SQRT(0.6D1)
4173
    CASE DEFAULT
4174
      ! TODO: Handle error somehow
4175
    END SELECT
4176
  END FUNCTION H1Basis_dPhi
4177

4178
  PURE FUNCTION H1Basis_VarPhi(k, x) RESULT(fval)
4179
    IMPLICIT NONE
4180

4181
    ! Parameters
4182
    INTEGER, INTENT(IN) :: k
4183
    REAL (KIND=dp), INTENT(IN) :: x
4184
    ! Return value
4185
    REAL (KIND=dp) :: fval
4186
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(k) _ELMER_LINEAR_REF(x) NOTINBRANCH
4187

4188
    SELECT CASE(k)
4189
    CASE(2)
4190
      fval = -SQRT(0.6D1)
4191
    CASE(3)
4192
      fval = -x * SQRT(0.10D2)
4193
    CASE(4)
4194
      fval = -DBLE(5 * x ** 2 - 1) * SQRT(0.14D2) / 0.4D1
4195
    CASE(5)
4196
      fval = -0.3D1 / 0.4D1 * DBLE(7 * x ** 2 - 3) * DBLE(x) * SQRT(0.2D1)
4197
    CASE(6)
4198
      fval = -DBLE(1 + (21 * x ** 2 - 14) * x ** 2) * SQRT(0.22D2) / 0.8D1
4199
    CASE(7)
4200
      fval = -DBLE(5 + (33 * x ** 2 - 30) * x ** 2) * DBLE(x) * &
4201
              SQRT(0.26D2) / 0.8D1
4202
    CASE(8)
4203
      fval = -DBLE(-5 + (135 + (429 * x ** 2 - 495) * x ** 2) * x ** 2) *&
4204
              SQRT(0.30D2) / 0.64D2
4205
    CASE(9)
4206
      fval = -DBLE(-35 + (385 + (715 * x ** 2 - 1001) * x ** 2) * x ** 2) * &
4207
              DBLE(x) * SQRT(0.34D2) / 0.64D2
4208
    CASE(10)
4209
      fval = -DBLE(7 + (-308 + (2002 + (2431 * x ** 2 - 4004) * x ** 2) * &
4210
              x ** 2) * x ** 2) * SQRT(0.38D2) / 0.128D3
4211
    CASE(11)
4212
      fval = -DBLE(63 + (-1092 + (4914 + (4199 * x ** 2 - 7956) * x ** 2) *&
4213
              x ** 2) * x ** 2) * DBLE(x) * SQRT(0.42D2) / 0.128D3
4214
    CASE(12)
4215
      fval = -DBLE(-21 + (1365 + (-13650 + (46410 + (29393 * x ** 2 - 62985) *&
4216
              x ** 2) * x ** 2) * x ** 2) * x ** 2) * SQRT(0.46D2) / 0.512D3
4217
    CASE(13)
4218
      fval = -0.5D1 / 0.512D3 * DBLE(-231 + (5775 + (-39270 + (106590 + &
4219
              (52003 * x ** 2 - 124355) * x ** 2) * x ** 2) * x ** 2) * &
4220
              x ** 2) * DBLE(x) * SQRT(0.2D1)
4221
    CASE(14)
4222
      fval = -0.3D1 / 0.1024D4 * DBLE(33 + (-2970 + (42075 + (-213180 + &
4223
              (479655 + (185725 * x ** 2 - 490314) * x ** 2) * x ** 2) *&
4224
              x ** 2) * x ** 2) * x ** 2) * SQRT(0.6D1)
4225
    CASE DEFAULT
4226
      ! TODO: Handle error somehow
4227
    END SELECT
4228
  END FUNCTION H1Basis_VarPhi
4229

4230
  PURE FUNCTION H1Basis_dVarPhi(k, x) RESULT(fval)
4231
    IMPLICIT NONE
4232

4233
    ! Parameters
4234
    INTEGER, INTENT(IN) :: k
4235
    REAL (KIND=dp), INTENT(IN) :: x
4236
    ! Return value
4237
    REAL (KIND=dp) :: fval
4238
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(k) _ELMER_LINEAR_REF(x) NOTINBRANCH
4239
    
4240
    SELECT CASE(k)
4241
    CASE(2)
4242
      fval = 0
4243
    CASE(3)
4244
      fval = -SQRT(0.10D2)
4245
    CASE(4)
4246
      fval = -0.5D1 / 0.2D1 * SQRT(0.14D2) * x
4247
    CASE(5)
4248
      fval = -0.9D1 / 0.4D1 * SQRT(0.2D1) * DBLE(7 * x ** 2 - 1)
4249
    CASE(6)
4250
      fval = -0.7D1 / 0.2D1 * SQRT(0.22D2) * DBLE(x) * DBLE(3 * x ** 2 - 1)
4251
    CASE(7)
4252
      fval = -0.5D1 / 0.8D1 * SQRT(0.26D2) * DBLE(1 + (33 * x ** 2 - 18) * &
4253
              x ** 2)
4254
    CASE(8)
4255
      fval = -0.9D1 / 0.32D2 * SQRT(0.30D2) * DBLE(x) * DBLE(15 + (143 * x ** 2 -&
4256
              110) * x ** 2)
4257
    CASE(9)
4258
      fval = -0.35D2 / 0.64D2 * SQRT(0.34D2) * DBLE(-1 + (33 + (143 * x ** 2 - 143) &
4259
              * x ** 2) * x ** 2)
4260
    CASE(10)
4261
      fval = -0.11D2 / 0.16D2 * SQRT(0.38D2) * DBLE(x) * DBLE(-7 + (91 + (221 * &
4262
              x ** 2 - 273) * x ** 2) * x ** 2)
4263
    CASE(11)
4264
      fval = -0.9D1 / 0.128D3 * SQRT(0.42D2) * DBLE(7 + (-364 + (2730 + (4199 * &
4265
              x ** 2 - 6188) * x ** 2) * x ** 2) * x ** 2)
4266
    CASE(12)
4267
      fval = -0.65D2 / 0.256D3 * SQRT(0.46D2) * DBLE(x) * DBLE(21 + (-420 + (2142 +&
4268
              (2261 * x ** 2 - 3876) * x ** 2) * x ** 2) * x ** 2)
4269
    CASE(13)
4270
      fval = -0.385D3 / 0.512D3 * SQRT(0.2D1) * DBLE(-3 + (225 + (-2550 + (9690 + &
4271
              (7429 * x ** 2 - 14535) * x ** 2) * x ** 2) * x ** 2) * x ** 2)
4272
    CASE(14)
4273
      fval = -0.45D2 / 0.256D3 * DBLE(x) * SQRT(0.6D1) * DBLE(-99 + (2805 + (-21318 + &
4274
              (63954 + (37145 * x ** 2 - 81719) * x ** 2) * x ** 2) * x ** 2) * x ** 2)
4275
    CASE DEFAULT
4276
      ! TODO: Handle error somehow
4277
    END SELECT
4278
  END FUNCTION H1Basis_dVarPhi
4279

4280
  PURE FUNCTION H1Basis_LegendreP(k, x) RESULT(fval)
4281
    IMPLICIT NONE
4282

4283
    ! Parameters
4284
    INTEGER, INTENT(IN) :: k
4285
    REAL (KIND=dp), INTENT(IN) :: x
4286
    ! Return value
4287
    REAL (KIND=dp) :: fval
4288
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(k) _ELMER_LINEAR_REF(x) NOTINBRANCH
4289
    
4290
    SELECT CASE(k)
4291
    CASE(0)
4292
      fval = 1
4293
    CASE(1)
4294

4295
      fval = x
4296
    CASE(2)
4297

4298
      fval = -0.1D1 / 0.2D1 + 0.3D1 / 0.2D1 * x ** 2
4299
    CASE(3)
4300

4301
      fval = (-0.3D1 / 0.2D1 + 0.5D1 / 0.2D1 * x ** 2) * x
4302
    CASE(4)
4303

4304
      fval = 0.3D1 / 0.8D1 + (-0.15D2 / 0.4D1 + 0.35D2 / 0.8D1 * x ** 2) * &
4305
              x ** 2
4306
    CASE(5)
4307

4308
      fval = (0.15D2 / 0.8D1 + (-0.35D2 / 0.4D1 + 0.63D2 / 0.8D1 * x ** 2) * &
4309
              x ** 2) * x
4310
    CASE(6)
4311

4312
      fval = -0.5D1 / 0.16D2 + (0.105D3 / 0.16D2 + (-0.315D3 / 0.16D2 + 0.231D3 / &
4313
              0.16D2 * x ** 2) * x ** 2) * x ** 2
4314
    CASE(7)
4315
      fval = (-0.35D2 / 0.16D2 + (0.315D3 / 0.16D2 + (-0.693D3 / 0.16D2 + 0.429D3 / &
4316
              0.16D2 * x ** 2) * x ** 2) * x ** 2) * x
4317
    CASE(8)                              
4318
      fval = 0.35D2 / 0.128D3 + (-0.315D3 / 0.32D2 + (0.3465D4 / 0.64D2 + &
4319
              (-0.3003D4 / 0.32D2 + 0.6435D4 / 0.128D3 * x ** 2) * x ** 2) * x ** 2) * x ** 2
4320
    CASE(9)
4321
      fval = (0.315D3 / 0.128D3 + (-0.1155D4 / 0.32D2 + (0.9009D4 / 0.64D2 + &
4322
              (-0.6435D4 / 0.32D2 + 0.12155D5 / 0.128D3 * x ** 2) * x ** 2) * &
4323
              x ** 2) * x ** 2) * x
4324
    CASE(10)
4325
      fval = -0.63D2 / 0.256D3 + (0.3465D4 / 0.256D3 + (-0.15015D5 / 0.128D3 + &
4326
              (0.45045D5 / 0.128D3 + (-0.109395D6 / 0.256D3 + 0.46189D5 / 0.256D3 *&
4327
              x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2
4328
    CASE(11)
4329
      fval = (-0.693D3 / 0.256D3 + (0.15015D5 / 0.256D3 + (-0.45045D5 / 0.128D3 +&
4330
              (0.109395D6 / 0.128D3 + (-0.230945D6 / 0.256D3 + 0.88179D5 / 0.256D3 * &
4331
              x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x
4332
    CASE(12)
4333
      fval = 0.231D3 / 0.1024D4 + (-0.9009D4 / 0.512D3 + (0.225225D6 / 0.1024D4 + &
4334
              (-0.255255D6 / 0.256D3 + (0.2078505D7 / 0.1024D4 + (-0.969969D6 / 0.512D3 +&
4335
              0.676039D6 / 0.1024D4 * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2
4336
    CASE(13)
4337
      fval = (0.3003D4 / 0.1024D4 + (-0.45045D5 / 0.512D3 + (0.765765D6 / 0.1024D4 + &
4338
              (-0.692835D6 / 0.256D3 + (0.4849845D7 / 0.1024D4 + (-0.2028117D7 / 0.512D3 + &
4339
              0.1300075D7 / 0.1024D4 * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) *&
4340
              x ** 2) * x
4341
    CASE(14)
4342
      fval = -0.429D3 / 0.2048D4 + (0.45045D5 / 0.2048D4 + (-0.765765D6 / 0.2048D4 +&
4343
              (0.4849845D7 / 0.2048D4 + (-0.14549535D8 / 0.2048D4 + (0.22309287D8 / &
4344
              0.2048D4 + (-0.16900975D8 / 0.2048D4 + 0.5014575D7 / 0.2048D4 * x ** 2) *&
4345
              x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2
4346
    CASE DEFAULT
4347
      ! TODO: Handle error somehow
4348
    END SELECT
4349
  END FUNCTION H1Basis_LegendreP
4350

4351
  PURE FUNCTION H1Basis_dLegendreP(k, x) RESULT(fval)
4352
    IMPLICIT NONE
4353

4354
    ! Parameters
4355
    INTEGER, INTENT(IN) :: k
4356
    REAL (KIND=dp), INTENT(IN) :: x
4357
    ! Return value
4358
    REAL (KIND=dp) :: fval
4359
    !_ELMER_OMP_DECLARE_SIMD UNIFORM(k) _ELMER_LINEAR_REF(x) NOTINBRANCH
4360
    
4361
    SELECT CASE(k)
4362
    CASE(0)
4363
      fval = 0
4364
    CASE(1)
4365
      fval = 1
4366
    CASE(2)
4367
      fval = 3 * x
4368
    CASE(3)
4369
      fval = 0.15D2 / 0.2D1 * x ** 2 - 0.3D1 / 0.2D1
4370
    CASE(4)
4371
      fval = 0.5D1 / 0.2D1 * DBLE(7 * x ** 2 - 3) * DBLE(x)
4372
    CASE(5)
4373
      fval = 0.15D2 / 0.8D1 + (-0.105D3 / 0.4D1 + 0.315D3 / 0.8D1 * x ** 2) *&
4374
              x ** 2
4375
    CASE(6)
4376
      fval = 0.21D2 / 0.8D1 * DBLE(5 + (33 * x ** 2 - 30) * x ** 2) * DBLE(x)
4377
    CASE(7)
4378
      fval = -0.35D2 / 0.16D2 + (0.945D3 / 0.16D2 + (-0.3465D4 / 0.16D2 + 0.3003D4 / &
4379
              0.16D2 * x ** 2) * x ** 2) * x ** 2
4380
    CASE(8)
4381
      fval = 0.9D1 / 0.16D2 * DBLE(-35 + (385 + (715 * x ** 2 - 1001) * x ** 2) * x ** 2) *&
4382
              DBLE(x)
4383
    CASE(9)
4384
      fval = 0.315D3 / 0.128D3 + (-0.3465D4 / 0.32D2 + (0.45045D5 / 0.64D2 + (-0.45045D5 / &
4385
              0.32D2 + 0.109395D6 / 0.128D3 * x ** 2) * x ** 2) * x ** 2) * x ** 2
4386
    CASE(10)
4387
      fval = 0.55D2 / 0.128D3 * DBLE(63 + (-1092 + (4914 + (4199 * x ** 2 - 7956) * x ** 2) * &
4388
              x ** 2) * x ** 2) * DBLE(x)
4389
    CASE(11)
4390
      fval = -0.693D3 / 0.256D3 + (0.45045D5 / 0.256D3 + (-0.225225D6 / 0.128D3 + &
4391
              (0.765765D6 / 0.128D3 + (-0.2078505D7 / 0.256D3 + 0.969969D6 / 0.256D3 *&
4392
              x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2
4393
    CASE(12)
4394
      fval = 0.39D2 / 0.256D3 * DBLE(-231 + (5775 + (-39270 + (106590 + (52003 * x ** 2 - &
4395
              124355) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * DBLE(x)
4396
    CASE(13)
4397
      fval = 0.3003D4 / 0.1024D4 + (-0.135135D6 / 0.512D3 + (0.3828825D7 / 0.1024D4 + &
4398
              (-0.4849845D7 / 0.256D3 + (0.43648605D8 / 0.1024D4 + (-0.22309287D8 / 0.512D3 + &
4399
              0.16900975D8 / 0.1024D4 * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2
4400
    CASE(14)
4401
      fval = 0.105D3 / 0.1024D4 * DBLE(429 + (-14586 + (138567 + (-554268 + (1062347 + (334305 *&
4402
              x ** 2 - 965770) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * x ** 2) * DBLE(x)
4403
    CASE DEFAULT
4404
      ! TODO: Handle error somehow
4405
    END SELECT
4406
  END FUNCTION H1Basis_dLegendreP
4407

4408
  ! To be deprecated
4409
  
4410
  ! WARNING: this is not a barycentric triangle
4411
  SUBROUTINE H1Basis_TriangleNodal(nvec, u, v, nbasismax, fval)
4412
    IMPLICIT NONE
4413

4414
    ! Parameters
4415
    INTEGER, INTENT(IN) :: nvec
4416
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v
4417
    INTEGER, INTENT(IN) :: nbasismax
4418
    ! Variables
4419
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
4420
    INTEGER :: j
4421

4422
    !_ELMER_OMP_SIMD
4423
    DO j=1,nvec
4424
      ! Node 1
4425
      fval(j,1) = 1-u(j)-v(j)
4426
      ! Node 2
4427
      fval(j,2) = u(j)
4428
      ! Node 3
4429
      fval(j,3) = v(j)
4430
    END DO
4431
  END SUBROUTINE H1Basis_TriangleNodal
4432

4433
  ! WARNING: this is not a barycentric triangle
4434
  SUBROUTINE H1Basis_dTriangleNodal(nvec, u, v, nbasismax, grad)
4435
    IMPLICIT NONE
4436

4437
    ! Parameters
4438
    INTEGER, INTENT(IN) :: nvec
4439
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v
4440
    INTEGER, INTENT(IN) :: nbasismax
4441
    ! Variables
4442
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
4443
    INTEGER :: j
4444

4445
    ! First coordinate (xi)
4446
    !_ELMER_OMP_SIMD
4447
    DO j=1,nvec
4448
      grad(j,1,1) = -1
4449
    END DO
4450
    !_ELMER_OMP_SIMD
4451
    DO j=1,nvec
4452
      grad(j,2,1) =  1
4453
    END DO
4454
    !_ELMER_OMP_SIMD
4455
    DO j=1,nvec
4456
      grad(j,3,1) =  0
4457
    END DO
4458
    
4459
    ! Second coordinate (eta)
4460
    !_ELMER_OMP_SIMD
4461
    DO j=1,nvec
4462
      grad(j,1,2) = -1
4463
    END DO
4464
    !_ELMER_OMP_SIMD
4465
    DO j=1,nvec
4466
      grad(j,2,2) =  0
4467
    END DO
4468
    !_ELMER_OMP_SIMD
4469
    DO j=1,nvec
4470
      grad(j,3,2) =  1
4471
    END DO
4472
  END SUBROUTINE H1Basis_dTriangleNodal
4473
  
4474
  ! WARNING: this is not a barycentric tetra
4475
  SUBROUTINE H1Basis_TetraNodal(nvec, u, v, w, nbasismax, fval)
4476
    IMPLICIT NONE
4477

4478
    ! Parameters
4479
    INTEGER, INTENT(IN) :: nvec
4480
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
4481
    INTEGER, INTENT(IN) :: nbasismax
4482
    ! Variables
4483
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
4484
    INTEGER :: j
4485

4486
    !_ELMER_OMP_SIMD
4487
    DO j=1,nvec
4488
      ! Node 1
4489
      fval(j,1) = 1-u(j)-v(j)-w(j)
4490
      ! Node 2
4491
      fval(j,2) = u(j)
4492
      ! Node 3
4493
      fval(j,3) = v(j)
4494
      ! Node 4
4495
      fval(j,4) = w(j)
4496
    END DO
4497
  END SUBROUTINE H1Basis_TetraNodal
4498

4499
  ! WARNING: this is not a barycentric tetra
4500
  SUBROUTINE H1Basis_dTetraNodal(nvec, u, v, w, nbasismax, grad)
4501
    IMPLICIT NONE
4502

4503
    ! Parameters
4504
    INTEGER, INTENT(IN) :: nvec
4505
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
4506
    INTEGER, INTENT(IN) :: nbasismax
4507
    ! Variables
4508
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
4509
    INTEGER :: j
4510

4511
    ! First coordinate (xi)
4512
    !_ELMER_OMP_SIMD
4513
    DO j=1,nvec
4514
      grad(j,1,1) = -1
4515
    END DO
4516
    !_ELMER_OMP_SIMD
4517
    DO j=1,nvec
4518
      grad(j,2,1) =  1
4519
    END DO
4520
    !_ELMER_OMP_SIMD
4521
    DO j=1,nvec
4522
      grad(j,3,1) =  0
4523
    END DO
4524
    !_ELMER_OMP_SIMD
4525
    DO j=1,nvec
4526
      grad(j,4,1) =  0
4527
    END DO
4528
    
4529
    ! Second coordinate (eta)
4530
    !_ELMER_OMP_SIMD
4531
    DO j=1,nvec
4532
      grad(j,1,2) = -1
4533
    END DO
4534
    !_ELMER_OMP_SIMD
4535
    DO j=1,nvec
4536
      grad(j,2,2) =  0
4537
    END DO
4538
    !_ELMER_OMP_SIMD
4539
    DO j=1,nvec
4540
      grad(j,3,2) =  1
4541
    END DO
4542
    !_ELMER_OMP_SIMD
4543
    DO j=1,nvec
4544
      grad(j,4,2) =  0
4545
    END DO
4546
    
4547
    ! Third coordinate (zeta)
4548
    !_ELMER_OMP_SIMD
4549
    DO j=1,nvec
4550
      grad(j,1,3) = -1
4551
    END DO
4552
    !_ELMER_OMP_SIMD
4553
    DO j=1,nvec
4554
      grad(j,2,3) =  0
4555
    END DO
4556
    !_ELMER_OMP_SIMD
4557
    DO j=1,nvec
4558
      grad(j,3,3) =  0
4559
    END DO
4560
    !_ELMER_OMP_SIMD
4561
    DO j=1,nvec
4562
      grad(j,4,3) =  1
4563
    END DO
4564
  END SUBROUTINE H1Basis_dTetraNodal
4565

4566
  ! WARNING: this is not a barycentric wedge
4567
  SUBROUTINE H1Basis_WedgeNodal(nvec, u, v, w, nbasismax, fval)
4568
    IMPLICIT NONE
4569

4570
    ! Parameters
4571
    INTEGER, INTENT(IN) :: nvec
4572
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
4573
    INTEGER, INTENT(IN) :: nbasismax
4574
    ! Variables
4575
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax), INTENT(INOUT) :: fval
4576
    REAL(Kind=dp), PARAMETER :: c = 1D0/2
4577
    INTEGER :: j
4578

4579
    !_ELMER_OMP_SIMD
4580
    DO j=1,nvec
4581
      ! Node 1
4582
      fval(j,1) = c*(1-u(j)-v(j))*(1-w(j))
4583
      ! Node 2
4584
      fval(j,2) = c*u(j)*(1-w(j))
4585
      ! Node 3
4586
      fval(j,3) = c*v(j)*(1-w(j))
4587
      ! Node 4
4588
      fval(j,4) = c*(1-u(j)-v(j))*(1+w(j))
4589
      ! Node 5
4590
      fval(j,5) = c*u(j)*(1+w(j))
4591
      ! Node 6
4592
      fval(j,6) = c*v(j)*(1+w(j))
4593
    END DO
4594
  END SUBROUTINE H1Basis_WedgeNodal
4595

4596
  ! WARNING: this is not a barycentric wedge
4597
  SUBROUTINE H1Basis_dWedgeNodal(nvec, u, v, w, nbasismax, grad)
4598
    IMPLICIT NONE
4599

4600
    ! Parameters
4601
    INTEGER, INTENT(IN) :: nvec
4602
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH), INTENT(IN) :: u,v,w
4603
    INTEGER, INTENT(IN) :: nbasismax
4604
    ! Variables
4605
    REAL(Kind=dp), DIMENSION(VECTOR_BLOCK_LENGTH,nbasismax,3), INTENT(INOUT) :: grad
4606
    REAL(Kind=dp), PARAMETER :: c = 1D0/2
4607
    INTEGER :: j
4608

4609
    ! First coordinate (xi)
4610
    !_ELMER_OMP_SIMD
4611
    DO j=1,nvec
4612
      grad(j,1,1) = -c*(1-w(j))
4613
      grad(j,2,1) =  c*(1-w(j))
4614
      grad(j,3,1) =  0
4615
      grad(j,4,1) = -c*(1+w(j))
4616
      grad(j,5,1) =  c*(1+w(j))
4617
      grad(j,6,1) =  0
4618
    END DO
4619
    
4620
    ! Second coordinate (eta)
4621
    !_ELMER_OMP_SIMD
4622
    DO j=1,nvec
4623
      grad(j,1,2) = -c*(1-w(j))
4624
      grad(j,2,2) =  0
4625
      grad(j,3,2) =  c*(1-w(j))
4626
      grad(j,4,2) = -c*(1+w(j))
4627
      grad(j,5,2) =  0
4628
      grad(j,6,2) =  c*(1+w(j))
4629
    END DO
4630
    
4631
    ! Third coordinate (zeta)
4632
    !_ELMER_OMP_SIMD
4633
    DO j=1,nvec
4634
      grad(j,1,3) = -c*(1-u(j)-v(j))
4635
      grad(j,2,3) = -c*u(j)
4636
      grad(j,3,3) = -c*v(j)
4637
      grad(j,4,3) =  c*(1-u(j)-v(j))
4638
      grad(j,5,3) =  c*u(j)
4639
      grad(j,6,3) =  c*v(j)
4640
    END DO
4641
  END SUBROUTINE H1Basis_dWedgeNodal
4642

4643
END MODULE H1Basis
4644

4645
Product

Resources

Company
